2- approx-exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 1: Linear algebra refresher\n",
    "<div id=\"fem:approx:exer:linalg1\"></div>\n",
    "\n",
    "Look up the topic of *vector space* in your favorite linear algebra\n",
    "book or search for the term at Wikipedia.\n",
    "\n",
    "\n",
    "**a)**\n",
    "Prove that vectors in the plane spanned by the vector $(a,b)$  form a vector space\n",
    "by showing that all the axioms of a vector space are satisfied.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The axioms of a vector space go as follows (see\n",
    "[Wikipedia](https://en.wikipedia.org/wiki/Vector_space), but\n",
    "we use slightly different notation below):\n",
    "\n",
    "1. The sum of $u$ and $v$, denoted by $u + v$, is in V.\n",
    "\n",
    "2. $u + v = v + u$.\n",
    "\n",
    "3. $(u + v) + w = u + (v + w)$.\n",
    "\n",
    "4. There is a *zero* vector $0$ in V such that $u + 0 = u$.\n",
    "\n",
    "5. For each $u$ in V, there is a vector $-u$ in V such that $u + (-u) = 0$.\n",
    "\n",
    "6. The scalar multiple of $u$ by $\\gamma$, denoted by $\\gamma u$, is in V.\n",
    "\n",
    "7. $\\gamma (u + v) = \\gamma u + \\gamma v$.\n",
    "\n",
    "8. $(\\gamma + \\delta)u = \\gamma u + \\delta u$ for scalar $\\gamma$ and $\\delta$.\n",
    "\n",
    "9. $\\gamma(\\delta u) = \\gamma\\delta u$.\n",
    "\n",
    "10. $1u = u$.\n",
    "\n",
    "We must show that each axiom is fulfilled by planar vectors and their\n",
    "mathematical rules. Let $u=(a,b)$ and $v=(c,d)$.\n",
    "\n",
    "Axiom 1:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u+v = (a,b)+(c,d)=(a+c, b+d),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is a planar vector and therefore in $V$.\n",
    "\n",
    "Axiom 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u+v = (a,b)+(c,d)=(a+c, b+d) = (c+a, d+b) = v + u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 3:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u+v) + w = ((a,b) + (c,d)) + (e,f) = (a,b) + ((c,d) + (e,f)) = u + (v+w){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 4: The $(0,0)$ vector is the 0 element,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u + 0 = (a,b) + (0,0) = (a+0, b+0) = (a,b) = u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 5: Let $-u$ element is $(-a,-b)$, so"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u + (-u) = (a,b) + (-a,-b) = (0,0) = 0{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 6:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\gamma u = \\gamma\\cdot(a,b)= (\\gamma a, \\gamma b),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is also a planar vector and therefore in $V$.\n",
    "\n",
    "Axiom 7:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\gamma (u + v) = \\gamma \\cdot((a,b) + (c,d) = \\gamma (a,b) + \\gamma (c,d){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 8:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(\\gamma + \\delta)u = (\\gamma + \\delta)(a,b) = \\gamma (a,b) + \\delta (a,b) = \\gamma u + \\delta u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 9:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\gamma (\\delta u) = \\gamma (\\delta\\cdot (a,b)) = \\gamma (\\delta a, \\delta b) = (\\gamma\\delta a, \\gamma\\delta b) = \\gamma\\delta u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 10:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "1u = 1(a,b) = (1\\cdot a, 1\\cdot b) = (a,b) = u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**b)**\n",
    "Prove that all linear functions of the form $ax+b$ constitute a vector space,\n",
    "$a,b\\in\\mathbb{R}$.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "Let $u=ax+b$ and $v=cx+d$. We verify each axiom.\n",
    "\n",
    "Axiom 1:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u+v = ax + b + cx + d = (a+c)x + (b + d),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is also a linear function and therefore in $V$.\n",
    "\n",
    "Axiom 2:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u+v = ax+b + cx+ d = cx+d + ax+ b = v + u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 3:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u+v) + w = (ax+b + cx+ d) + ex + f = ax+b + cx+ d + ex + f\n",
    "ax+b + (cx+ d + ex + f) = u + (v+w){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 4: The $0x+0=0$ function is the 0 element,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u + 0 = ax + b + 0 = ax + b = u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 5: Let $-u$ element is $-ax-b$, so"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u + (-u) = ax + b + (-ax - b) = 0{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 6:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\gamma u = \\gamma(ax +b)= \\gamma ax + \\gamma b,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is also a linear function and therefore in $V$.\n",
    "\n",
    "Axiom 7:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\gamma (u + v) = \\gamma (ax+b + cx+ d) = \\gamma ax + \\gamma b + \\gamma cx + \\gamma d = \\gamma(ax + b) + \\gamma(cd + d) = \\gamma u + \\gamma v{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 8:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(\\gamma + \\delta)u = (\\gamma + \\delta)(ax + b) = \\gamma (ax+b) + \\delta (ax+b) = \\gamma u + \\delta u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 9:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\gamma (\\delta u) = \\gamma (\\delta(ax+b))) = \\gamma\\delta u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Axiom 10:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "1u = 1(ax+b) = ax+b = u{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**c)**\n",
    "Show that all quadratic functions of the form $1 + ax^2 + bx$ *do not*\n",
    "constitute a vector space.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "Let $u=ax^2+bx + 1$ and $v=cx^2 + dx + 1$. We try to verify each axiom.\n",
    "\n",
    "Axiom 1:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u+v = ax^2 + bx + 1 + cx^2 + dx + 1 = (a+c)x^2 + (b+d)x + 2,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "but this quadratic function is not in $V$ because the constant term is 2\n",
    "and not 1. Consequently, quadratic functions of the particular\n",
    "form $1 + ax^2 + bx$ do not constitute a vector space.\n",
    "The more general form  $ax^2 + bx +c$ for arbitrary constants $a$, $b$, and\n",
    "$c$ makes functions that span a function space.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**d)**\n",
    "Check out\n",
    "the topic of *inner product spaces*. Suggest a possible inner product\n",
    "for the space of all linear functions of the form $ax+b$, $a,b\\in\\mathbb{R}$,\n",
    "defined on some interval $\\Omega=[A,B]$.\n",
    "Show that this particular inner product satisfies the\n",
    "general requirements of an inner product in a vector space.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "According to [Wikipedia](https://en.wikipedia.org/wiki/Inner_product_space),\n",
    "an inner product space, with an inner product $(u,v)$, has three axioms\n",
    "(we drop the possibility of complex numbers and assume everything is real):\n",
    "\n",
    "1. Symmetry: $(u,v)$ = $(v,u)$\n",
    "\n",
    "2. Linearity in the first argument: $\\gamma u, v) = \\gamma (u,v)$\n",
    "\n",
    "3. Positive-definiteness: $(u,u)\\geq 0$ and $(u,u)=0$ implies $u=0$\n",
    "\n",
    "A possible inner product for linear functions on a domain $\\Omega$ is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u,v) = \\int_A^B uv{\\, \\mathrm{d}x} = \\int_A^B (ax+b)(cx+d){\\, \\mathrm{d}x}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Symmetry is obvious since $uv=vu$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u,v) = \\int_A^B uv{\\, \\mathrm{d}x} = \\int_A^B vu d = (v,u){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Linearity in the first argument:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(\\gamma u,v) = \\gamma\\int_A^B(\\gamma u)v{\\, \\mathrm{d}x} = \\gamma\\int_A^B uv{\\, \\mathrm{d}x} = \\gamma (u,v){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Positive-definiteness:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u,u) = \\int_A^B(ax+b)^2{\\, \\mathrm{d}x}\\geq 0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "since the integral of a function $f(x)\\geq$ must be greater than or equal to\n",
    "zero, and in particular,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u,u)=0\\quad\\Rightarrow\\quad \\int_A^B(ax+b)^2{\\, \\mathrm{d}x} = 0,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "implies $ax+b=0$.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "\n",
    "Filename: `linalg1`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 2: Approximate a three-dimensional vector in a plane\n",
    "<div id=\"fem:approx:exer:vec:3Dby2D\"></div>\n",
    "\n",
    "Given $\\boldsymbol{f} = (1,1,1)$ in $\\mathbb{R}^3$, find the best approximation vector\n",
    "$\\boldsymbol{u}$ in the plane spanned by the unit vectors $(1,0)$ and $(0,1)$.\n",
    "Repeat the calculations using the vectors $(2,1)$ and $(1,2)$.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "We have the vector $\\pmb{f} = (1,1,1)$.\n",
    "Our aim is to approximate this with a vector in the vector-space spanned by the unit vectors $\\pmb{\\phi_0}$ and $\\pmb{\\phi_1}$. We seek a solution $u = c_0 \\phi_0 + c_1 \\phi_1$. To do this we use the least-square-method and solve the equation $\\pmb{A c} = \\pmb{b}$. Or written out:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left[\\begin{array}{cc}\n",
    "A_{0,0} & A_{0,1} \\\\A_{1,0} & A_{1,1}\n",
    "\\end{array}\\right]\n",
    "\\left[ \\begin{array}{c} c_0 \\\\ c_1 \\end{array} \\right]\n",
    "=\n",
    "\\left[ \\begin{array}{c} b_0 \\\\ b_1 \\end{array} \\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $A_{i,j} = ( \\pmb{\\phi_i}, \\pmb{\\phi_j} )$ and $b_i = ( \\pmb{\\phi_i}, \\pmb{f} )$.\n",
    "\n",
    "We start with $\\pmb{\\phi_0} = (0,1)$, $\\pmb{\\phi_1} = (1,0)$.\n",
    "Calculations give"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A_{0,0} &= ([1,0],[1,0]) = 1,\\\\\n",
    "A_{0,1} &= ([1,0],[0,1]) = 0,\n",
    "A_{1,0} &= ([0,1],[1,0]) = 0,\n",
    "A_{1,1} &= ([0,1],[0,1]) = 1,\n",
    "b_0 &= ([1,0],[1,1,1]) = 1,\n",
    "b_0 &= ([0,1],[1,1,1]) = 1.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The result becomes $ c_0 = c_1 =1$, and hence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u = 1\\cdot\\pmb{\\phi_0}  + 1\\cdot\\pmb{\\phi_1}.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we proceed with\n",
    "$\\pmb{\\phi_0} = (2,1)$, $\\pmb{\\phi_1} = (1,2)$.\n",
    "Calculations give"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A_{0,0} &= ([2,1],[2,1]) = 2\\cdot 2 + 1\\cdot 1 = 5,\n",
    "A_{0,1} &= ([2,1],[1,2]) = 2\\cdot 1 + 1\\cdot 2 = 4,\n",
    "A_{1,0} &= ([1,2],[2,1]) = 1\\cdot 2 + 2\\cdot 1 = 4,\n",
    "A_{1,1} &= ([1,2],[1,2]) = 1\\cdot 1 + 2\\cdot 2 = 5,\n",
    "b_0 &= ([2,1],[1,1,1]) = 2\\cdot 1 + 1\\cdot 1 = 3,\n",
    "b_0 &= ([1,2],[1,1,1]) = 1\\cdot 1 + 2\\cdot 1 = 3.\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Solving for the $c_i$ values results in $ c_0 = c_1 =\\frac{1}{3}$ and\n",
    "hence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u = \\frac{1}{3}\\cdot \\pmb{\\phi_0}  + \\frac{1}{3}\\cdot \\pmb{\\phi_1} .\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `vec111_approx`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 3: Approximate a parabola by a sine\n",
    "<div id=\"fem:approx:exer:parabola_sine\"></div>\n",
    "\n",
    "Given the function $f(x)=1 + 2x(1-x)$ on $\\Omega=[0,1]$, we want to\n",
    "find an approximation in the function space"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "V = \\hbox{span}\\{1, \\sin(\\pi x)\\}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)**\n",
    "Sketch or plot $f(x)$. Think intuitively how an expansion in terms\n",
    "of the basis functions of $V$, ${\\psi}_0(x)=1$, ${\\psi}_1(x)=\\sin(\\pi x)$,\n",
    "can be constructed to yield a best approximation to $f$. Or\n",
    "phrased differently, see if you can guess the coefficients $c_0$ and $c_1$\n",
    "in the expansion"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = c_0{\\psi}_0 + c_1{\\psi}_1 = c_0 + c_1\\sin (\\pi x){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compute the $L^2$ error $||f-u||_{L^2}=(\\int_0^1(f-u)^2{\\, \\mathrm{d}x})^{1/2}$.\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint.**\n",
    "If you make a mesh function `e` of the error\n",
    "on some mesh with uniformly spaced coordinates in\n",
    "the array `xc`, the integral can be approximated as `np.sqrt(dx*np.sum(e**2))`,\n",
    "where `dx=xc[0]-xc[1]` is the mesh spacing and `np` is an alias for the\n",
    "`numpy` module in Python.\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The function ${\\psi}_0=1$ can be used to \"take care of\" of the\n",
    "constant 1 in $f$, while ${\\psi}_1=\\sin(\\pi x)$ can approximate\n",
    "the parabola. The maximum of $f$ is 3/2, so we should use\n",
    "$u(x) = 1\\cdot{\\psi}_0 + \\frac{1}{2}{\\psi}_1$.\n",
    "\n",
    "Plotting and the computation of the $L^2$ error can be done by"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "xc = np.linspace(0, 1, 101)  # x coordinates for plotting\n",
    "\n",
    "def f(x):\n",
    "    return 1 + 2*x*(1-x)\n",
    "\n",
    "import sympy as sym\n",
    "x = sym.symbols('x')\n",
    "psi_0 = 1\n",
    "psi_1 = sym.sin(sym.pi*x)\n",
    "\n",
    "half = sym.Rational(1,2)\n",
    "u = 1*psi_0 + half*psi_1\n",
    "\n",
    "# How to combine c_0*psi_0 + c_1*psi_1 to match f?\n",
    "# Intuitively, c_0=c_1=1...\n",
    "# Turn u to function so we can plot and compute with it\n",
    "u = sym.lambdify([x], u, modules='numpy')\n",
    "\n",
    "print('L2 error of intuitive approximation:', end=' ')\n",
    "e = f(xc) - u(xc)\n",
    "dx = xc[1] - xc[0]\n",
    "print(np.sqrt(dx*np.sum(e**2)))\n",
    "\n",
    "plt.plot(xc, f(xc), 'r--')\n",
    "plt.plot(xc, u(xc), 'b-')\n",
    "plt.legend(['exact', 'intuitive approximation'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $L^2$ error becomes $0.0179$.\n",
    "\n",
    "<!-- dom:FIGURE: [fig/parabola_sin_a.png, width=500 frac=0.7] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/parabola_sin_a.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**b)**\n",
    "Perform the hand calculations for a least squares approximation.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The least squares method ends up with a linear system where the\n",
    "coefficient matrix has entries $A_{i,j}=\\int_0^1{\\psi}_i{\\psi}_j{\\, \\mathrm{d}x}$\n",
    "and the right-hand side has entries $b_i=\\int_0^1{\\psi}_i{\\, \\mathrm{d}x}$.\n",
    "We can use `sympy` do carry out the integrals:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Do the calculations in the least squares or project method\n",
    "A = sym.zeros(2, 2)\n",
    "b = sym.zeros(2, 1)\n",
    "A[0,0] = sym.integrate(psi_0*psi_0, (x, 0, 1))\n",
    "A[0,1] = sym.integrate(psi_0*psi_1, (x, 0, 1))\n",
    "A[1,0] = A[0,1]\n",
    "A[1,1] = sym.integrate(psi_1*psi_1, (x, 0, 1))\n",
    "b[0] = sym.integrate(f(x)*psi_0, (x, 0, 1))\n",
    "b[1] = sym.integrate(f(x)*psi_1, (x, 0, 1))\n",
    "print(('A:', A))\n",
    "print(('b:', b))\n",
    "c = A.LUsolve(b)\n",
    "c = [sym.simplify(c[i,0]) for i in range(c.shape[0])]\n",
    "print(('c:', c, [c_.evalf() for c_ in c]))\n",
    "u = c[0]*psi_0 + c[1]*psi_1\n",
    "print(('u:', u))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the results, we get the linear system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\\begin{matrix}1 & \\frac{2}{\\pi}\\\\\\frac{2}{\\pi} & \\frac{1}{2}\\end{matrix}\\right)\n",
    "\\left(\\begin{bmatrix}\n",
    "\\frac{- 24 \\pi^{2} - 96 + 4 \\pi^{4}}{3 \\pi^{2} \\left(-8 + \\pi^{2}\\right)}\\\\\n",
    "\\frac{- 4 \\pi^{2} + 48}{3 \\pi \\left(-8 + \\pi^{2}\\right)}\n",
    "\\end{bmatrix}\\right) =\n",
    "\\left(\\begin{matrix}\n",
    "\\frac{4}{3}\\\\\n",
    "\\frac{8}{\\pi^{3}} + \\frac{2}{\\pi}\n",
    "\\end{matrix}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting best approximation reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = \\frac{4 \\left(- \\pi^{2} + 12\\right) \\sin{\\left (\\pi x \\right )}}{3 \\pi \\left(-8 + \\pi^{2}\\right)} + \\frac{- 24 \\pi^{2} - 96 + 4 \\pi^{4}}{3 \\pi^{2} \\left(-8 + \\pi^{2}\\right)}\\approx 1.025 + 0.484\\sin(\\pi x)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The $L^2$ error turns out to be $0.00876$.\n",
    "To conclude, the least squares method is slightly better than the intuition\n",
    "in this case.\n",
    "\n",
    "<!-- dom:FIGURE: [fig/parabola_sin_b.png, width=500 frac=0.7] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/parabola_sin_b.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "Filename: `parabola_sin`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 4: Approximate the exponential function by power functions\n",
    "<div id=\"fem:approx:exer:exp:powers\"></div>\n",
    "\n",
    "Let $V$ be a function space with basis functions $x^i$,\n",
    "$i=0,1,\\ldots,N$.  Find the best approximation to $f(x)=\\exp(-x)$ on\n",
    "$\\Omega =[0,8]$ among all functions in $V$ for $N=2,4,6$. Illustrate\n",
    "the three approximations in three separate plots.\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint.**\n",
    "Apply the `lest_squares` and\n",
    "`comparison_plot` functions in the `approx1D.py` module as these\n",
    "make the exercise easier to solve.\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "A suitable code is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "from approx1D import least_squares, comparison_plot\n",
    "import matplotlib.pyplot as plt\n",
    "import sympy as sym\n",
    "from math import factorial\n",
    "import numpy as np\n",
    "\n",
    "x = sym.Symbol('x')\n",
    "f = sym.exp(-x)\n",
    "\n",
    "Omega = [0, 8]\n",
    "\n",
    "for N in 2,4,6:\n",
    "    psi = [x**i for i in range(N+1)]\n",
    "    u, c = least_squares(f,psi,Omega)\n",
    "    print((N, u))\n",
    "    comparison_plot(f, u, Omega, filename='tmp_exp_%d' % N,\n",
    "                    plot_title='N=%d' % N)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the case $N=2$ the program prints the following, here edited for clearer\n",
    "reading:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        A:\n",
    "        Matrix([[8, 32, 512/3], [32, 512/3, 1024], [512/3, 1024, 32768/5]])\n",
    "        b:\n",
    "        Matrix([[-exp(-8) + 1], [-9*exp(-8) + 1], [-82*exp(-8) + 2]])\n",
    "        \n",
    "        f: exp(-x)\n",
    "        u: x**2*(-465*exp(-8)/4096 + 105/4096) + x*(-141/512 +\n",
    "           405*exp(-8)/512) - 111*exp(-8)/128 + 87/128\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:FIGURE: [fig/exp_powers_N2.png, width=500 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/exp_powers_N2.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `exp_powers`.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- Taylor: these polynomials go so far off on [0,8] that it is not a -->\n",
    "<!-- good idea to add them. -->\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 5: Approximate the sine function by power functions\n",
    "<div id=\"fem:approx:exer:sin:powers\"></div>\n",
    "\n",
    "In this exercise we want to approximate the sine function by polynomials\n",
    "of order $N+1$. Consider two bases:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "V_1 &= \\{x, x^3, x^5, \\ldots, x^{N-2}, x^N \\}, \\\\\n",
    "V_2 &= \\{1,x,x^2,x^3,\\ldots, x^N\\}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The basis $V_1$ is motivated by the fact that the Taylor polynomial\n",
    "approximation to the sine function has only odd powers, while $V_2$\n",
    "is motivated by the assumption that including the even powers could\n",
    "improve the approximation in a least-squares setting.\n",
    "\n",
    "Compute the best approximation to $f(x)=\\sin(x)$ among all functions\n",
    "in $V_1$ and $V_2$ on two domains of increasing sizes: $\\Omega_{1,k} =\n",
    "[0, k\\pi]$, $k=2,3\\ldots,6$ and $\\Omega_{2,k} = [-k\\pi /2, k\\pi/2]$,\n",
    "$k=2,3,4,5$.  Make plots for all combinations of $V_1$, $V_2$,\n",
    "$\\Omega_1$, $\\Omega_2$, $k=2,3,\\ldots,6$.\n",
    "\n",
    "Add a plot of the $N$-th degree Taylor polynomial approximation of\n",
    "$\\sin(x)$ around $x=0$.\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint.**\n",
    "You can make a loop over $V_1$ and $V_2$, a loop over\n",
    "$\\Omega_1$ and $\\Omega_2$, and a loop over $k$. Inside the loops,\n",
    "call the functions `least_squares` and\n",
    "`comparison_plot` from the `approx1D` module.\n",
    "$N=7$ is a suggested value.\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "Suitable code is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "from approx1D import least_squares, comparison_plot\n",
    "from math import pi\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = sym.Symbol('x')\n",
    "f = sym.sin(x)\n",
    "N = 7\n",
    "psi_bases = [[x**i for i in range(1, N+1, 2)],  # V_1\n",
    "             [x**i for i in range(0, N+1)]]     # V_2\n",
    "symbolic = False\n",
    "\n",
    "for V, psi in enumerate(psi_bases):\n",
    "    for domain_no in range(1, 3):\n",
    "        for k in range(2, 6):\n",
    "            if symbolic:\n",
    "                Omega = [0, k*sym.pi] if domain_no == 1 else \\\n",
    "                        [-k*sym.pi/2, k*sym.pi/2]\n",
    "            else:\n",
    "                # cannot use sym.pi with numerical sympy computing\n",
    "                Omega = [0, k*pi] if domain_no == 1 else \\\n",
    "                        [-k*pi/2, k*pi/2]\n",
    "\n",
    "            u, c = least_squares(f, psi, Omega, symbolic=symbolic)\n",
    "\n",
    "            comparison_plot(\n",
    "                f, u, Omega,\n",
    "                ymin=-2, ymax=2,\n",
    "                filename='tmp_N%d_V%dOmega%dk%d' %\n",
    "                (N, V, k, domain_no),\n",
    "                plot_title='sin(x) on [0,%d*pi/2] by %s' %\n",
    "                (k, ','.join([str(p) for p in psi])))\n",
    "            # Need to kill the plot to proceed!\n",
    "        for ext in 'png', 'pdf':\n",
    "            cmd = 'doconce combine_images -2 ' + \\\n",
    "                  ' '.join(['tmp_N%d_V%dOmega%dk%d.' %\n",
    "                            (N, V, k, domain_no) + ext\n",
    "                            for k in range(2, 6)]) + \\\n",
    "                  ' sin_powers_N%d_V%d_Omega%d.' % (N, V, domain_no) + ext\n",
    "            print(cmd)\n",
    "            os.system(cmd)\n",
    "\n",
    "# Show the standard Taylor series approximation\n",
    "from math import factorial, pi\n",
    "import time\n",
    "Omega = [0, 12*pi/2.]\n",
    "u = 0\n",
    "for k in range(0,N+1):\n",
    "    u = u + ((-1)**k*x**(1+2*k))/float(factorial(1+2*k))\n",
    "# Shorter: u = sum(((-1)**k*x**(1+2*k))/float(factorial(1+2*k))\n",
    "# for k in range(0,10))\n",
    "comparison_plot(f, u, Omega, 'sin_taylor%d' % k,\n",
    "                ymin=-1.5, ymax=1.5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The odd powers ($V_1$ space) behave not so good on $\\Omega_{1,k}$,\n",
    "but better on $\\Omega_{2,k}$:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/sin_powers_N7_V1_Omega1.png, width=800 frac=1] $V_1$ space, $\\Omega_{1,k}$ domain. -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p>$V_1$ space, $\\Omega_{1,k}$ domain.</p>\n",
    "<img src=\"fig/sin_powers_N7_V1_Omega1.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [fig/sin_powers_N7_V1_Omega2.png, width=800 frac=1] $V_1$ space, $\\Omega_{2,k}$ domain. -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p>$V_1$ space, $\\Omega_{2,k}$ domain.</p>\n",
    "<img src=\"fig/sin_powers_N7_V1_Omega2.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "Including also even powers ($V_2$ space) is clearly much better:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/sin_powers_N7_V2_Omega1.png, width=800 frac=1] $V_2$ space, $\\Omega_{1,k}$ domain. -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p>$V_2$ space, $\\Omega_{1,k}$ domain.</p>\n",
    "<img src=\"fig/sin_powers_N7_V2_Omega1.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [fig/sin_powers_N7_V2_Omega2.png, width=800 frac=1] $V_2$ space, $\\Omega_{2,k}$ domain. -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p>$V_2$ space, $\\Omega_{2,k}$ domain.</p>\n",
    "<img src=\"fig/sin_powers_N7_V2_Omega2.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "Comparison with a standard Taylor series shows that it is very inferior\n",
    "as an approximation over the entire domain, but much more accurate\n",
    "close to the origin (as expected, since the Taylor series is constructed\n",
    "with this property, while the least squares method tries to find a good\n",
    "approximation over the entire domain).\n",
    "\n",
    "<!-- dom:FIGURE: [fig/sin_taylor7.png, width=500 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/sin_taylor7.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `sin_powers`.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- Solveig explanations based on f-B approx -->\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 6: Approximate a steep function by sines\n",
    "<div id=\"fem:approx:exer:tanh:sine1\"></div>\n",
    "\n",
    "Find the best approximation of $f(x) = \\tanh (s(x-\\pi))$ on\n",
    "$[0, 2\\pi]$ in the space $V$ with basis\n",
    "${\\psi}_i(x) = \\sin((2i+1)x)$, $i\\in{\\mathcal{I}_s} = \\{0,\\ldots,N\\}$.\n",
    "Make a movie showing how $u=\\sum_{j\\in{\\mathcal{I}_s}}c_j{\\psi}_j(x)$\n",
    "approximates $f(x)$ as $N$ grows. Choose $s$ such that $f$ is\n",
    "steep ($s=20$ is appropriate).\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint 1.**\n",
    "One may naively call the `least_squares_orth` and `comparison_plot`\n",
    "from the `approx1D` module in a loop and extend the basis with\n",
    "one new element in each pass. This approach\n",
    "implies a lot of recomputations.\n",
    "A more efficient strategy is to let `least_squares_orth`\n",
    "compute with only one basis function at a time and accumulate\n",
    "the corresponding `u` in the total solution.\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint 2.**\n",
    "`ffmpeg` or `avconv` may skip frames when plot files are combined to\n",
    "a movie. Since there are few files and we want to see each of them,\n",
    "use `convert` to make an animated GIF file (`-delay 200` is\n",
    "suitable).\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The code may read"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "from approx1D import least_squares_orth, comparison_plot\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "x = sym.Symbol('x')\n",
    "\n",
    "# Naive approach: (not utilizing the fact that i+1 computations can\n",
    "# make use of i computations)\n",
    "def naive(f, s, Omega, N=10):\n",
    "    psi = []\n",
    "    for i in range(N+1):\n",
    "        psi.append(sym.sin((2*i+1)*x))\n",
    "        u, c = least_squares_orth(f, psi, Omega, symbolic=False)\n",
    "        comparison_plot(f, u, Omega, 'tmp_sin%02dx' % i,\n",
    "                        legend_loc='upper left', show=True)\n",
    "\n",
    "# Efficient approach: compute just the matrix diagonal\n",
    "def efficient(f, s, Omega, N=10):\n",
    "    u = 0\n",
    "    for i in range(N+1):\n",
    "        psi = [sym.sin((2*i+1)*x)]\n",
    "        next_term, c = least_squares_orth(f, psi, Omega, False)\n",
    "        u = u + next_term\n",
    "        comparison_plot(f, u, Omega, 'tmp_sin%02dx' % i,\n",
    "                        legend_loc='upper left', show=False,\n",
    "                        plot_title='s=%g, i=%d' % (s, i))\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    s = 20  # steepness\n",
    "    f = sym.tanh(s*(x-sym.pi))\n",
    "    from math import pi\n",
    "    Omega = [0, 2*pi]  # sym.pi did not work here\n",
    "    efficient(f, s, Omega, N=10)\n",
    "    # Make movie\n",
    "    # avconv/ffmpeg skips frames, use convert instead (few files)\n",
    "    cmd = 'convert -delay 200 tmp_sin*.png tanh_sines_approx.gif'\n",
    "    os.system(cmd)\n",
    "    # Make static plots, 3 figures on 2 lines\n",
    "    for ext in 'pdf', 'png':\n",
    "        cmd = 'doconce combine_images %s -3 ' % ext\n",
    "        cmd += 'tmp_sin00x tmp_sin01x tmp_sin02x tmp_sin04x '\n",
    "        cmd += 'tmp_sin07x tmp_sin10x tanh_sines_approx'\n",
    "        os.system(cmd)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:FIGURE: [fig/tanh_sines.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_sines.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `tanh_sines`.\n",
    "\n",
    "<!-- Closing remarks for this Problem -->\n",
    "\n",
    "### Remarks\n",
    "\n",
    "Approximation of a discontinuous (or steep) $f(x)$ by sines,\n",
    "results in slow convergence and oscillatory behavior of the\n",
    "approximation close to the abrupt changes in $f$.\n",
    "This is known as the [Gibb's phenomenon](http://en.wikipedia.org/wiki/Gibbs_phenomenon).\n",
    "\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 7: Approximate a steep function by sines with boundary adjustment\n",
    "<div id=\"fem:approx:exer:tanh:sine3\"></div>\n",
    "\n",
    "We study the same approximation problem as in\n",
    "[Problem 6: Approximate a steep function by sines](#fem:approx:exer:tanh:sine1). Since ${\\psi}_i(0)={\\psi}_i(2\\pi)=0$\n",
    "for all $i$, $u=0$ at the boundary points $x=0$ and $x=2\\pi$, while\n",
    "$f(0)=-1$ and $f(2\\pi)=1$. This discrepancy at the boundary can be\n",
    "removed by adding a boundary function $B(x)$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = B(x) + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(x),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $B(x)$ has the right boundary values: $B(x_L)=f(x_L)$ and\n",
    "$B(x_R)=f(x_R)$, with $x_L=0$ and $x_R=2\\pi$ as the boundary points.\n",
    "A linear choice of $B(x)$ is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "B(x) = \\frac{(x_R-x)f(x_L) + (x-x_L)f(x_R)}{x_R-x_L}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)**\n",
    "Use the basis\n",
    "${\\psi}_i(x) = \\sin((i+1)x)$, $i\\in{\\mathcal{I}_s} = \\{0,\\ldots,N\\}$\n",
    "and plot $u$ and $f$ for $N=16$. (It suffices to make plots for even $i$.)\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "With a boundary term $B(x)$ we call `least_squares_orth` with\n",
    "`f-B` as right-hand side function, and we must remember to add $B$ to $u$.\n",
    "\n",
    "We can extend the code from [Problem 6: Approximate a steep function by sines](#fem:approx:exer:tanh:sine1) and\n",
    "let the function `efficient` handle different choices of basis.\n",
    "Appropriate code for all three subexercises is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "from approx1D import least_squares_orth, comparison_plot\n",
    "import matplotlib.pyplot as plt\n",
    "x = sym.Symbol('x')\n",
    "\n",
    "def efficient(f, B, s, Omega, N=10, basis='a'):\n",
    "    u = B\n",
    "    for i in range(N+1):\n",
    "        if basis == 'a':\n",
    "            psi = [sym.sin((i+1)*x)]\n",
    "        elif basis == 'b':\n",
    "            psi = [sym.sin((2*i+1)*x)]\n",
    "        elif basis == 'c':\n",
    "            psi = [sym.sin(2*(i+1)*x)]\n",
    "        next_term, c = least_squares_orth(f-B, psi, Omega, False)\n",
    "        u = u + next_term\n",
    "        # Make only plot for i even\n",
    "        if i % 2 == 0:\n",
    "            comparison_plot(f, u, Omega, 'tmp_sin%02dx' % i,\n",
    "                            legend_loc='upper left', show=False,\n",
    "                            plot_title='s=%g, i=%d' % (s, i))\n",
    "\n",
    "\n",
    "if __name__ == '__main__':\n",
    "    s = 20  # steepness\n",
    "    f = sym.tanh(s*(x-sym.pi))\n",
    "    from math import pi\n",
    "    Omega = [0, 2*pi]  # sym.pi did not work here\n",
    "\n",
    "    # sin((i+1)*x) basis\n",
    "    xL = Omega[0]\n",
    "    xR = Omega[1]\n",
    "    B = ((xR-x)*f.subs(x, xL) + (x-xL)*f.subs(x, xR))/(xR-xL)\n",
    "    for exercise in 'a', 'b', 'c':\n",
    "        efficient(f, B, s, Omega, N=16, basis=exercise)\n",
    "        # Make movie\n",
    "        cmd = 'convert -delay 200 tmp_sin*.png '\n",
    "        cmd += 'tanh_sines_boundary_term_%s.gif' % exercise\n",
    "        os.system(cmd)\n",
    "        # Make static plots, 3 figures on 2 lines\n",
    "        for ext in 'pdf', 'png':\n",
    "            cmd = 'doconce combine_images %s -3 ' % ext\n",
    "            cmd += 'tmp_sin00x tmp_sin02x tmp_sin04x tmp_sin08x '\n",
    "            cmd += 'tmp_sin12x tmp_sin16x '\n",
    "            cmd += 'tanh_sines_boundary_term_%s' % exercise\n",
    "            os.system(cmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:FIGURE: [fig/tanh_sines_boundary_term_a.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_sines_boundary_term_a.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**b)**\n",
    "Use the basis from [Problem 6: Approximate a steep function by sines](#fem:approx:exer:tanh:sine1),\n",
    "${\\psi}_i(x) = \\sin((2i+1)x)$, $i\\in{\\mathcal{I}_s} = \\{0,\\ldots,N\\}$.\n",
    "(It suffices to make plots for even $i$.)\n",
    "Observe that the approximation converges to a piecewise\n",
    "linear function!\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "<!-- dom:FIGURE: [fig/tanh_sines_boundary_term_b.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_sines_boundary_term_b.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**c)**\n",
    "Use the basis\n",
    "${\\psi}_i(x) = \\sin(2(i+1)x)$, $i\\in{\\mathcal{I}_s} = \\{0,\\ldots,N\\}$,\n",
    "and observe that the approximation converges to a piecewise\n",
    "constant function.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "<!-- dom:FIGURE: [fig/tanh_sines_boundary_term_c.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_sines_boundary_term_c.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "Filename: `tanh_sines_boundary_term`.\n",
    "\n",
    "<!-- Closing remarks for this Problem -->\n",
    "\n",
    "### Remarks\n",
    "\n",
    "The strange results in b) and c) are due to the choice of\n",
    "basis. In b), ${\\varphi}_i(x)$ is an odd function around\n",
    "$x=\\pi/2$ and $x=3\\pi/2$. No combination of basis functions\n",
    "is able to approximate the flat regions of $f$.\n",
    "All basis functions in c) are even around $x=\\pi/2$ and $x=3\\pi/2$,\n",
    "but odd at $x=0,\\pi,2\\pi$. With all the sines represented, as in a),\n",
    "the approximation is not constrained with a particular symmetry\n",
    "behavior.\n",
    "\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 8: Fourier series as a least squares approximation\n",
    "<div id=\"fem:approx:exer:Fourier\"></div>\n",
    "\n",
    "\n",
    "**a)**\n",
    "Given a function $f(x)$ on an interval $[0,L]$, look up the formula\n",
    "for the coefficients $a_j$ and $b_j$ in the Fourier series of $f$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "f(x) = \\frac{1}{2}a_0 +\n",
    "\\sum_{j=1}^\\infty a_j\\cos \\left(j\\frac{2\\pi x}{L}\\right)\n",
    "+ \\sum_{j=1}^\\infty b_j\\sin \\left(j\\frac{2\\pi x}{L}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "From [Wikipedia](https://en.wikipedia.org/wiki/Fourier_series) we\n",
    "have"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "a_j &= \\frac{2}{L}\\int_0^P f(x)\\cos\\left(j\\frac{2\\pi x}{L}\\right) {\\, \\mathrm{d}x},\\\\\n",
    "b_j &= \\frac{2}{L}\\int_0^P f(x)\\sin\\left(j\\frac{2\\pi x}{L}\\right) {\\, \\mathrm{d}x}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**b)**\n",
    "Let an infinite-dimensional vector space $V$ have the basis functions\n",
    "$\\cos j\\frac{2\\pi x}{L}$ for $j=0,1,\\dots,\\infty$ and\n",
    "$\\sin j\\frac{2\\pi x}{L}$ for $j=1,\\dots,\\infty$.  Show that the least squares\n",
    "approximation method from the section [Approximation principles](#fem:approx:global) leads to a\n",
    "linear system whose solution coincides with the standard formulas for\n",
    "the coefficients in a Fourier series of $f(x)$ (see also\n",
    "the section [Fourier series](#fem:approx:global:Fourier)).\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint.**\n",
    "You may choose"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:approx:exer:Fourier:basis\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "{\\psi}_{2i} = \\cos\\left( i\\frac{2\\pi}{L}x\\right),\\quad\n",
    "{\\psi}_{2i+1} = \\sin\\left( i\\frac{2\\pi}{L}x\\right),\n",
    "\\label{fem:approx:exer:Fourier:basis} \\tag{67}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $i=0,1,\\ldots,N\\rightarrow\\infty$.\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The entries in the linear system arising from the least squares method\n",
    "are $A_{i,j}=\\int_0^L{\\psi}_i{\\psi}_j{\\, \\mathrm{d}x}$ and $b_i=\\int_0^L\n",
    "f(x){\\psi}_i {\\, \\mathrm{d}x}$. To avoid name clash between the right-hand side\n",
    "components of the linear system and the $b_i$ coefficients in the\n",
    "Fourier series, we use the symbol $q_i$ for the former.\n",
    "With the basis functions in\n",
    "([67](#fem:approx:exer:Fourier:basis)) we get four different types\n",
    "of integrals:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A_{2i,2j} &= \\int_0^L \\cos\\left( i\\frac{2\\pi}{L}x\\right)\n",
    "\\cos\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x} = A_{2j,2i},\\\\\n",
    "A_{2i, 2j+1} &= \\int_0^L \\cos\\left( i\\frac{2\\pi}{L}x\\right)\n",
    "\\sin\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x},\\\\\n",
    "A_{2i+1,2j} &= \\int_0^L \\sin\\left( i\\frac{2\\pi}{L}x\\right)\n",
    "\\cos\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x},\\\\\n",
    "A_{2i+1, 2j+1} &= \\int_0^L \\sin\\left( i\\frac{2\\pi}{L}x\\right)\n",
    "\\sin\\left( j\\frac{2\\pi}{L}x\\right)dx,\\\\\n",
    "q_{2i} &= \\int_0^L f(x)\\cos\\left( i\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x},\\\\\n",
    "q_{2i+1} &= \\int_0^L f(x)\\sin\\left( i\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, the sine and cosine basis functions are orthogonal on $[0,L]$.\n",
    "We have in general"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "\\int_0^L\n",
    "\\cos\\left( i\\frac{2\\pi}{L}x\\right) \\cos\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x}\n",
    "&=0,\\quad i\\neq j,\\\\\n",
    "\\int_0^L\n",
    "\\cos\\left( i\\frac{2\\pi}{L}x\\right) \\cos\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x}\n",
    "&=\\frac{L}{2},\\quad i= j\\neq 0,\\\\\n",
    "\\int_0^L\n",
    "\\cos\\left( i\\frac{2\\pi}{L}x\\right) \\cos\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x}\n",
    "&=L,\\quad i= j= 0,\\\\\n",
    "\\int_0^L\n",
    "\\sin\\left( i\\frac{2\\pi}{L}x\\right) \\sin\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x}\n",
    "&=0,\\quad i\\neq j,\\\\\n",
    "\\int_0^L\n",
    "\\sin\\left( i\\frac{2\\pi}{L}x\\right) \\sin\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x}\n",
    "&=\\frac{L}{2},\\quad i= j,\\\\\n",
    "\\int_0^L\n",
    "\\cos\\left( i\\frac{2\\pi}{L}x\\right) \\sin\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x}\n",
    "&=0{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These results imply that only diagonal terms in the coefficient matrix\n",
    "are different from zero. We have"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "A_{0,0} &= L,\\\\\n",
    "A_{2i,2i} &= \\frac{L}{2},\\quad i>0,\\\\\n",
    "A_{2i+1,2i+1} &= \\frac{L}{2}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The unknown vector with components $c_i$ must be arranged as"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(a_0, b_1, a_1, b_2, a_2, b_3, \\ldots){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then get"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{0,0}a_0=q_0,\\quad A_{1,1}b_1=q_1,\\quad A_{2,2}a_1=q_2,\\quad A_{3,3}b_2=q_3,\\ldots\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These equations lead to the formulas"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "a_0 &= \\frac{1}{L}\\int_0^P f(x){\\, \\mathrm{d}x},\\\\\n",
    "b_1 &= \\frac{2}{L}\\int_0^P f(x)\\sin\\left( \\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x},\\\\\n",
    "a_1 &= \\frac{2}{L}\\int_0^P f(x)\\cos\\left( \\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x},\\\\\n",
    "b_2 &= \\frac{2}{L}\\int_0^P f(x)\\sin\\left( 2\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x},\\\\\n",
    "a_2 &= \\frac{2}{L}\\int_0^P f(x)\\cos\\left( 2\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which can be generalized to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "a_0 &= \\frac{1}{L}\\int_0^P f(x){\\, \\mathrm{d}x},\\\\\n",
    "a_j &= \\frac{2}{L}\\int_0^P f(x)\\cos\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x},\\ j>0,\\\\\n",
    "b_j &= \\frac{2}{L}\\int_0^P f(x)\\sin\\left( j\\frac{2\\pi}{L}x\\right){\\, \\mathrm{d}x},\\ j>0,\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and these are the standard formulas for the Fourier coefficients in a)\n",
    "if we recognize that the $a_0$ above is twice the $a_0$ in the\n",
    "expressions in a).\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**c)**\n",
    "Choose $f(x) = H(x-\\frac{1}{2})$ on $\\Omega=[0,1]$, where $H$ is the\n",
    "Heaviside function: $H(x)=0$ for $x<0$, $H(x)=1$ for $x>0$ and\n",
    "$H(0)=\\frac{1}{2}$. Find the coefficients $a_j$ and $b_j$ in the\n",
    "Fourier series for $f(x)$. Plot the sum for $j=2N+1$, where $N=5$ and\n",
    "$N=100$.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The formulas give"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "a_0 &= 2\\int_0^1 f(x){\\, \\mathrm{d}x} = 2\\int_{\\frac{1}{2}}^1 {\\, \\mathrm{d}x},\\\\\n",
    "a_j &= 2\\int_0^1 f(x)\\cos\\left( 2j\\pi x\\right){\\, \\mathrm{d}x}\n",
    "= 2\\int_{\\frac{1}{2}}^1 \\cos\\left( 2j\\pi x\\right){\\, \\mathrm{d}x},\\\\\n",
    "b_j &= 2\\int_{\\frac{1}{2}}^1 \\sin\\left( 2j\\pi x\\right){\\, \\mathrm{d}x}{\\thinspace .}\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The integrals are readily computed by `sympy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "j = sym.symbols('k', integer=True)\n",
    "x = sym.symbols('x', real=True)\n",
    "I = integrate(cos(2*j*pi*x), (x,Rational(1,2),1))\n",
    "I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = integrate(cos(2*0*pi*x), (x,Rational(1,2),1))\n",
    "I"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = integrate(sin(2*j*pi*x), (x,Rational(1,2),1))\n",
    "I"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that we have the series"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = \\frac{1}{2} + 2\\sum_{j=1}^\\infty \\frac{(-1)^j - 1}{2\\pi j}\n",
    "\\sin\\left( 2j\\pi x\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We only get a nonzero coefficient for $j$ odd:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = \\frac{1}{2} -2\\sum_{k=1}^\\infty \\frac{1}{(2k+1)\\pi}\n",
    "\\sin\\left( 2(2k+1)\\pi x\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Appropriate computer code for visualizing the series goes like"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import pi\n",
    "from numpy import sin\n",
    "\n",
    "def Heaviside_series(x, N):\n",
    "    s = 0.5\n",
    "    for k in range(N):\n",
    "        s += -2.0/((2*k+1)*pi)*sin(2*(2*k+1)*pi*x)\n",
    "    return s\n",
    "\n",
    "x = np.linspace(0, 1, 1001)\n",
    "for N in 5, 100:\n",
    "    H = Heaviside_series(x, N)\n",
    "    plt.figure()\n",
    "    plt.plot(x, H)\n",
    "    plt.legend(['$N=%d$' % N], loc='upper left')\n",
    "    plt.savefig('tmp_%d.png' % N)\n",
    "    plt.savefig('tmp_%d.pdf' % N)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:FIGURE: [fig/Fourier_Heaviside.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/Fourier_Heaviside.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "We clearly see the Gibbs' phenomenon: oscillations and overshoot around\n",
    "the point of discontinuity in the function we try to approximate.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "\n",
    "Filename: `Fourier_ls`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 9: Approximate a steep function by Lagrange polynomials\n",
    "<div id=\"fem:approx:exer:tanh:Lagrange\"></div>\n",
    "\n",
    "Use interpolation with uniformly distributed\n",
    "points and Chebychev nodes to approximate"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "f(x) = -\\tanh(s(x-\\frac{1}{2})),\\quad x\\in [0,1],\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "by Lagrange polynomials for $s=5$ and $s=20$, and $N=3,7,11,15$.\n",
    "Combine $2\\times 2$ plots of the approximation for the four\n",
    "$N$ values, and create such figures for the four combinations of\n",
    "$s$ values and point types.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The following code does the work (symbolically):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "sys.path.insert(0, os.path.join(os.pardir, 'src'))\n",
    "from approx1D import interpolation, comparison_plot\n",
    "from Lagrange import Lagrange_polynomials\n",
    "import sympy as sym\n",
    "\n",
    "x = sym.Symbol('x')\n",
    "Omega = [0,1]\n",
    "N_values = 3, 7, 11, 15\n",
    "\n",
    "for s in 5, 20:\n",
    "    f = -sym.tanh(s*(x-0.5))  # sympy expression\n",
    "    for distribution in 'uniform', 'Chebyshev':\n",
    "        for N in N_values:\n",
    "            phi, points = Lagrange_polynomials(\n",
    "                x, N, Omega,\n",
    "                point_distribution=distribution)\n",
    "\n",
    "            u, c = interpolation(f, phi, points)\n",
    "            filename = 'tmp_tanh_%d_%d_%s' % (N, s, distribution)\n",
    "            comparison_plot(f, u, Omega, filename,\n",
    "                            plot_title='s=%g, N=%d, %s points' %\n",
    "                            (s, N, distribution))\n",
    "        # Combine plot files (2x2)\n",
    "        for ext in 'png', 'pdf':\n",
    "            cmd = 'doconce combine_images ' + ext + ' '\n",
    "            cmd += ' '.join([\n",
    "                'tmp_tanh_%d_%d_%s' % (N, s, distribution)\n",
    "                for N in N_values])\n",
    "            cmd += ' tanh_Lagrange_%s_s%s' % (distribution, s)\n",
    "            os.system(cmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For a smooth function ($s=5$), the difference between uniform points and\n",
    "Chebyshev nodes is not substantial:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/tanh_Lagrange_uniform_s5.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_Lagrange_uniform_s5.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [fig/tanh_Lagrange_Chebyshev_s5.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_Lagrange_Chebyshev_s5.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "However, for a steep function ($s=20$) the overshoot and oscillations\n",
    "associated with uniform points must be considered unacceptable for\n",
    "larger $N$ values:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/tanh_Lagrange_uniform_s20.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_Lagrange_uniform_s20.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "Switching to Chebyshev points does give a great improvement, but\n",
    "we still have oscillatory approximations:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/tanh_Lagrange_Chebyshev_s20.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_Lagrange_Chebyshev_s20.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `tanh_Lagrange`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 10: Approximate a steep function by Lagrange polynomials and regression\n",
    "<div id=\"fem:approx:exer:tanh:Lagrange:regression\"></div>\n",
    "\n",
    "Redo [Problem 9: Approximate a steep function by Lagrange polynomials](#fem:approx:exer:tanh:Lagrange), but apply a regression\n",
    "method with $N$-degree Lagrange polynomials and $2N+1$\n",
    "data points. Recall that\n",
    "[Problem 9: Approximate a steep function by Lagrange polynomials](#fem:approx:exer:tanh:Lagrange) applies $N+1$ points and\n",
    "the resulting approximation interpolates $f$ at these points, while\n",
    "a regression method with more points does not interpolate $f$ at the\n",
    "data points.\n",
    "Do more points and a regression method help reduce\n",
    "the oscillatory behavior of Lagrange polynomial approximations?\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "We start out with the program from\n",
    "[Problem 9: Approximate a steep function by Lagrange polynomials](#fem:approx:exer:tanh:Lagrange). This time we need\n",
    "to call `Lagrange_polynomails` twice: first to compute the ${\\psi}(x)$\n",
    "functions (of degree $N$) and then to compute the data points corresponding\n",
    "to a uniform or Chebyshev distribution of $2N+1$ nodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys, os\n",
    "sys.path.insert(0, os.path.join(os.pardir, 'src'))\n",
    "from approx1D import regression, comparison_plot\n",
    "from Lagrange import Lagrange_polynomials\n",
    "import sympy as sym\n",
    "import numpy as np\n",
    "\n",
    "x = sym.Symbol('x')\n",
    "Omega = [0,1]\n",
    "N_values = 3, 7, 11, 15\n",
    "\n",
    "for s in 5, 20:\n",
    "    f = -sym.tanh(s*(x-0.5))  # sympy expression\n",
    "    for distribution in 'uniform', 'Chebyshev':\n",
    "        for N in N_values:\n",
    "            # Compute the points from a 2*N Lagrange polynomial\n",
    "            dummy, points = Lagrange_polynomials(\n",
    "                x, 2*N, Omega,\n",
    "                point_distribution=distribution)\n",
    "            # Compute phi from N points Lagrange polynomial\n",
    "            phi, dummy = Lagrange_polynomials(\n",
    "                x, N, Omega,\n",
    "                point_distribution=distribution)\n",
    "            points = np.array(points, dtype=float)\n",
    "            point_values = -np.tanh(s*(points-0.5))\n",
    "\n",
    "            u, c = regression(f, phi, points)\n",
    "            filename = 'tmp_tanh_%d_%d_%s' % (N, s, distribution)\n",
    "            comparison_plot(f, u, Omega, filename,\n",
    "                            plot_title='s=%g, N=%d, %s points' %\n",
    "                            (s, N, distribution),\n",
    "                            points=points, point_values=point_values,\n",
    "                            points_legend='%s points' % (2*N))\n",
    "        # Combine plot files (2x2)\n",
    "        for ext in 'png', 'pdf':\n",
    "            cmd = 'doconce combine_images ' + ext + ' '\n",
    "            cmd += ' '.join([\n",
    "                'tmp_tanh_%d_%d_%s' % (N, s, distribution)\n",
    "                for N in N_values])\n",
    "            cmd += ' tanh_Lagrange_regr_%s_s%s' % (distribution, s)\n",
    "            os.system(cmd)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An important point is to convert `points` to a `numpy` array using\n",
    "`dtype=float`. Leaving out this second argument makes an array\n",
    "of objects of symbolic expressions, and we cannot apply `tanh` to it.\n",
    "\n",
    "The oscillatory behavior is much reduced using more points and a\n",
    "regression method, and the difference between uniform and Chebyshev\n",
    "points is minor, even in the steep case $s=20$:\n",
    "\n",
    "<!-- dom:FIGURE: [fig/tanh_Lagrange_regr_uniform_s20.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_Lagrange_regr_uniform_s20.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- dom:FIGURE: [fig/tanh_Lagrange_regr_Chebyshev_s20.png, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/tanh_Lagrange_regr_Chebyshev_s20.png\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `tanh_Lagrange_regression`.\n",
    "\n",
    "<!-- --- end exercise --- -->"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}