5- varform-fe-exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 5: Compute the deflection of a cable with 2 P1 elements\n",
    "<div id=\"fem:deq:exer:cable:2P1\"></div>\n",
    "\n",
    "Solve the problem for $u$ in [Exercise 2: Compute the deflection of a cable with sine functions](#fem:deq:exer:tension:cable)\n",
    "using two P1 linear elements. Incorporate the condition $u(0)=0$\n",
    "by two methods: 1) excluding the unknown at $x=0$, 2) keeping the unknown\n",
    "at $x=0$, but modifying the linear system.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "From [Exercise 2: Compute the deflection of a cable with sine functions](#fem:deq:exer:tension:cable),\n",
    "the Galerkin method, after integration by parts, reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u',v')=-(1,v)\\quad\\forall v\\in V{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have two elements,\n",
    "$\\Omega^{(0)}=[0,\\frac{1}{2}]$ and $\\Omega^{(1)}=[\\frac{1}{2},1]$.\n",
    "\n",
    "**Method 1: Excluding the unknown at $x=0$.**\n",
    "Since $u(0)=0$, we exclude the value at $x=0$ as degree of freedom in\n",
    "the linear system. (There is no need for any boundary function.)\n",
    "The expansion reads $u=c_0{\\varphi}_1(x) + c_1{\\varphi}_2(x)$.\n",
    "The element matrix has then only one entry in the first element,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde A^{(0)} = \\frac{1}{h}(1){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From element 1 we get the usual element matrix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde A^{(1)} = \\frac{1}{h}\n",
    "\\left(\\begin{array}{rr}\n",
    "1 & -1\\\\ \n",
    "-1 & 1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The element vector in element 0 becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde b^{(0)} = \\frac{h}{2}(-1),\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "while the second element gives a contribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde b^{(1)} = \\frac{h}{2}\n",
    "\\left(\\begin{array}{c}\n",
    "-1 \\\\ \n",
    "-1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assembling the contributions gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{1}{h}\n",
    "\\left(\\begin{array}{cc}\n",
    "2 & -1\\\\ \n",
    "-1 & 1\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "c_1\n",
    "\\end{array}\\right)\n",
    "=\n",
    "- \\frac{h}{2}\n",
    "\\left(\\begin{array}{c}\n",
    "2 \\\\ \n",
    "1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that $h=\\frac{1}{2}$.\n",
    "Solving this system yields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "c_0 = -\\frac{3}{8},\\quad c_1=-\\frac{1}{2}\\quad\\Rightarrow\\quad u=-\\frac{3}{8}{\\varphi}_1(x)-\\frac{1}{2}{\\psi}_2(x){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Evaluating the exact solution for $x=\\frac{1}{2}$ and $x=1$, we get $3/8$\n",
    "and $1/2$, respectively, a result which shows\n",
    "that the numerical solution with P1 is exact at the three node points.\n",
    "The difference between the numerical and exact solution is that the\n",
    "numerical solution varies linearly over the two elements while the\n",
    "exact solution is quadratic.\n",
    "\n",
    "**Method 2: Modifying the linear system.**\n",
    "Now we let $c_i$ correspond to the value at node $x_{i}$, i.e.,\n",
    "all known Dirichlet values become part of the linear system.\n",
    "The expansion is now simply $u=\\sum_{i=0}^2c_i{\\varphi}_i(x)$, with\n",
    "three unknowns $c_0$, $c_1$, and $c_2$.\n",
    "Now the element matrix in the first and second element are equal.\n",
    "The same is true for the element vectors.\n",
    "Assembling yields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{1}{h}\n",
    "\\left(\\begin{array}{ccc}\n",
    "1 & -1 & 0\\\\ \n",
    "-1 & 2 & -1\\\\ \n",
    "0 & -1 & 1\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "c_1\\\\ \n",
    "c_2\n",
    "\\end{array}\\right)\n",
    "=\n",
    "- \\frac{h}{2}\n",
    "\\left(\\begin{array}{c}\n",
    "1\\\\ \n",
    "2\\\\ \n",
    "1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next step is to modify the linear system to implement the\n",
    "Dirichlet condition $c_0=0$. We first multiply by $h=\\frac{1}{2}$ and replace the\n",
    "first equation by $c_0=0$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\\begin{array}{ccc}\n",
    "1 & 0 & 0\\\\ \n",
    "-1 & 2 & -1\\\\ \n",
    "0 & -1 & 1\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "c_0 \\\\ \n",
    "c_1\\\\ \n",
    "c_2\n",
    "\\end{array}\\right)\n",
    "=\n",
    "- \\left(\\begin{array}{c}\n",
    "0\\\\ \n",
    "\\frac{1}{4}\\\\ \n",
    "\\frac{1}{8}\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the remaining $2\\times 2$ system is identical to the one\n",
    "previously solved, and the solution is the same."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u = 0{\\varphi}_0(x) - \\frac{3}{8}{\\varphi}_1(x) - \\frac{1}{2}{\\psi}_2(x){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `cable_2P1`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 6: Compute the deflection of a cable with 1 P2 element\n",
    "<div id=\"fem:deq:exer:cable:1P2\"></div>\n",
    "\n",
    "Solve the problem for $u$ in [Exercise 2: Compute the deflection of a cable with sine functions](#fem:deq:exer:tension:cable)\n",
    "using one P2 element with quadratic basis functions.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The P2 basis functions on a reference element $[-1,1]$ are"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "{\\tilde{\\varphi}}_0(X) &= \\frac{1}{2} (X-1)X\n",
    "\\\\ \n",
    "{\\tilde{\\varphi}}_1(X) &= 1 - X^2\n",
    "\\\\ \n",
    "{\\tilde{\\varphi}}_2(X) &= \\frac{1}{2} (X+1)X\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The element matrix and vector are easily calculated by some lines with\n",
    "`sympy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "X, h = sym.symbols('X h')\n",
    "half = sym.Rational(1, 2)\n",
    "psi = [half*(X-1)*X, 1-X**2, half*(X+1)*X]\n",
    "dpsi_dX = [sym.diff(psi[r], X) for r in range(len(psi))]\n",
    "\n",
    "# Element matrix\n",
    "# (2/h)*dpsi_dX[r]*(2/h)*dpsi_dX[s]*h/2\n",
    "import numpy as np\n",
    "d = 2\n",
    "# Use a numpy matrix with general objects to hold A\n",
    "A = np.empty((d+1, d+1), dtype=object)\n",
    "for r in range(d+1):\n",
    "    for s in range(d+1):\n",
    "        integrand = dpsi_dX[r]*dpsi_dX[s]*2/h\n",
    "        A[r,s] = sym.integrate(integrand, (X, -1, 1))\n",
    "print(A)\n",
    "\n",
    "# Element vector\n",
    "# f*psi[r]*h/2, f=1\n",
    "d = 2\n",
    "b = np.empty(d+1, dtype=object)\n",
    "for r in range(d+1):\n",
    "    integrand = -psi[r]*h/2\n",
    "    b[r] = sym.integrate(integrand, (X, -1, 1))\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The formatted element matrix and vector output becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        [[7/(3*h) -8/(3*h) 1/(3*h)]\n",
    "         [-8/(3*h) 16/(3*h) -8/(3*h)]\n",
    "         [1/(3*h) -8/(3*h) 7/(3*h)]]\n",
    "        [-h/6 -2*h/3 -h/6]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "or in mathematical notation:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde A^{(e)}=\\frac{1}{3h}\n",
    "\\left(\\begin{array}{ccc}\n",
    "7 & -8 & 1\\\\ \n",
    "-8 & 16 & -8\\\\ \n",
    "1 & -8 & 7\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(e)} = - \\frac{h}{6}\n",
    "\\left(\\begin{array}{c}\n",
    "1\\\\ \n",
    "4\\\\ \n",
    "1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Method 1: Excluding the unknown at $x=0$.**\n",
    "The expansion is $u=c_0{\\varphi}_1(x) + c_1{\\varphi}_2(x)$. The element matrix\n",
    "corresponding to the first element excludes contributions associated with\n",
    "the unknown at the left node, i.e., we exclude row and column 0. In the\n",
    "element vector, we exclude the first entry."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde A^{(0)}=\\frac{1}{3h}\n",
    "\\left(\\begin{array}{cc}\n",
    "16 & -8\\\\ \n",
    "-8 & 7\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(e)} = - \\frac{h}{6}\n",
    "\\left(\\begin{array}{c}\n",
    "4\\\\ \n",
    "1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, $h=1$.\n",
    "The solution of the linear system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{1}{3h}\n",
    "\\left(\\begin{array}{cc}\n",
    "16 & -8\\\\ \n",
    "-8 & 7\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "c_1\\\\ \n",
    "c_2\n",
    "\\end{array}\\right)\n",
    "=\n",
    "- \\frac{h}{6}\n",
    "\\left(\\begin{array}{c}\n",
    "4\\\\ \n",
    "1\n",
    "\\end{array}\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is $c_1=3/8$ and $c_2=1/2$. As for P1 elements in [Exercise 5: Compute the deflection of a cable with 2 P1 elements](#fem:deq:exer:cable:2P1), the values at the nodes are exact, but\n",
    "this time the variation between the nodes is quadratic, i.e., exact.\n",
    "One P2 element produces the complete, exact solution.\n",
    "\n",
    "**Method 2: Modifying the linear system.**\n",
    "This time the expansion reads $u=\\sum_{i=0}^2 c_i{\\varphi}_i(x)$ with\n",
    "three unknowns $c_0$, $c_1$, and $c_2$. The linear system consists\n",
    "of the complete $3\\times 3$ element matrix and the corresponding\n",
    "element vector:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{1}{3h}\n",
    "\\left(\\begin{array}{ccc}\n",
    "7 & -8 & 1\\\\ \n",
    "-8 & 16 & -8\\\\ \n",
    "1 & -8 & 7\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "c_0\\\\ \n",
    "c_1\\\\ \n",
    "c_2\n",
    "\\end{array}\\right)\n",
    "=\n",
    "- \\frac{h}{6}\n",
    "\\left(\\begin{array}{c}\n",
    "1\\\\ \n",
    "4\\\\ \n",
    "1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The boundary condition is incorporated by replacing the first equation\n",
    "by $c_0=0$, but prior to taking that action, we multiply by $3h$ and\n",
    "insert $h=1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\\begin{array}{ccc}\n",
    "1 & 0 & 0   \\\\ \n",
    "-8 & 16 & -8\\\\ \n",
    "1 & -8 & 7\n",
    "\\end{array}\\right)\n",
    "\\left(\\begin{array}{c}\n",
    "c_0\\\\ \n",
    "c_1\\\\ \n",
    "c_2\n",
    "\\end{array}\\right)\n",
    "=\n",
    "\\left(\\begin{array}{c}\n",
    "0\\\\ \n",
    "- 2\\\\ \n",
    "- \\frac{1}{2}\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Realizing that $c_0=0$, which means we can remove the first column of\n",
    "the system, shows that the equations are the same as above and hence\n",
    "that the solution is identical.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `cable_1P2`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 7: Compute the deflection of a cable with a step load\n",
    "<div id=\"fem:deq:exer:cable:stepload\"></div>\n",
    "\n",
    "We consider the deflection of a tension cable as described in\n",
    "[Exercise 2: Compute the deflection of a cable with sine functions](#fem:deq:exer:tension:cable): $w''=\\ell$, $w(0)=w(L)=0$.\n",
    "Now the load is discontinuous:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\ell (x) =\\left\\lbrace\\begin{array}{ll}\n",
    "\\ell_1, & x < L/2,\\\\ \n",
    "\\ell_2, & x \\geq L/2\n",
    "\\end{array}\\right.\\quad x\\in [0,L]\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This load is not symmetric\n",
    "with respect to the midpoint $x=L/2$ so the solution loses its symmetry.\n",
    "Scaling the problem by introducing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\bar x = \\frac{x}{L/2},\\quad u = \\frac{w}{w_c},\\quad\\bar\\ell = \\frac{\\ell - \\ell_1}{\\ell_2 - \\ell_1}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This leads to a scaled problem on $[0,2]$ (we rename $\\bar x$ as $x$\n",
    "for convenience):"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u'' = \\bar\\ell(x) = \\left\\lbrace\\begin{array}{ll}\n",
    "1, & x < 1,\\\\ \n",
    "0, & x \\geq 1\n",
    "\\end{array}\\right.\n",
    "\\quad x\\in (0,1),\\quad u(0)=0,\\ u(2)=0\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)**\n",
    "Find the analytical solution of the problem.\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint.**\n",
    "Integrate the equation separately for $x < 1$ and $x>1$. Use\n",
    "the conditions that $u$ and $u'$ must be continuous at $x=1$.\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "For $x < 1$ we get $u_1(x) = C_1x + C_2$, and the boundary condition\n",
    "$u_1(0)=0$ implies $C_2=0$. For $x>1$ we get $u_2(x)=\\frac{1}{2} x^2 + C_3x + C_4$.\n",
    "Continuity of $u'(1)$ leads to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C_1 = 1 + C_3,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and continuity of $u(1)$ means"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "C_1 = \\frac{1}{2} + C_3 + C_4,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "while the condition $u_2(2)=0$ gives the third equation we need:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "2 + 2C_3 + C_4 = 0{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We use `sympy` to solve them:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sympy import symbols, Rational, solve\n",
    "C1, C3, C4 = symbols('C1 C3 C4')\n",
    "solve([C1 - 1 - C3,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u(x) = \\left\\lbrace\\begin{array}{ll}\n",
    "-\\frac{1}{4}x, & x\\leq 1,\\\\ \n",
    "\\frac{1}{2} x^2 - \\frac{5}{4}x + \\frac{1}{2}, & x\\geq 1\n",
    "\\end{array}\\right.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:FIGURE: [fig/cable_discont_load_u_exact.png, width=500 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/cable_discont_load_u_exact.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**b)**\n",
    "Use ${\\psi}_i = \\sin((i+1)\\frac{\\pi x}{2})$,\n",
    "$i=0,\\ldots,N$ and the Galerkin method to find an approximate\n",
    "solution $u=\\sum_j c_j{\\psi}_j$.\n",
    "Plot how fast the coefficients $c_j$ tend to zero (on a log scale).\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The Galerkin formulation of the problem becomes"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "(u',v') = -(\\bar\\ell, v) = \\left\\lbrace\\begin{array}{ll}\n",
    "0, & x\\leq 1,\\\\ \n",
    "-(1,v), & x\\geq 1\n",
    "\\end{array}\\right.\\quad\\forall v\\in V{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A requirement is that $v(0)=v(2)=0$ because of the boundary conditions\n",
    "on $u$. The chosen basis functions\n",
    "fulfill this requirement for any integer $i$. Inserting\n",
    "$u=\\sum_{j=0}^N c_j{\\psi}_j$ and $v={\\psi}_i$, $i=0,\\ldots,N$,\n",
    "gives as usual the linear system $\\sum_j A_{i,j}c_j = b_i$, $i=0,\\ldots,N$,\n",
    "where"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i,j} = (i+1)(j+1)\\frac{\\pi^2}{4}\\int_0^2 \\cos((i+1)\\frac{\\pi x}{2})\n",
    "\\cos((j+1)\\frac{\\pi x}{2})dx{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The cosine functions are orthogonal on $[0,2]$ so $A_{i,j}=0$ for\n",
    "$i\\neq j$, while $A_{i,i}$ is computed (e.g., by `sympy`) as in\n",
    "[Exercise 2: Compute the deflection of a cable with sine functions](#fem:deq:exer:tension:cable), part e. The result is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "A_{i,i} = (i+1)^2\\frac{\\pi^2}{4}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The right-hand side is"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "b_i =\n",
    "-\\int_1^2 \\sin((i+1)\\frac{\\pi x}{2})dx\n",
    "= \\frac{2}{\\pi (i+1)}(\\cos((i+1)\\pi) - \\cos((i+1)\\pi/2)){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(Trying to do the integral in `sympy` gives a complicated expression that\n",
    "needs discussion - it\n",
    "is easier to do all calculations by hand.)\n",
    "We have that $\\cos((i+1)\\pi = -1$ for $i$ even and\n",
    "$\\cos((i+1)\\pi = 1$ for $i$ odd, while $\\cos((i+1)\\pi/2)$ is\n",
    "discussed in\n",
    "[Exercise 2: Compute the deflection of a cable with sine functions](#fem:deq:exer:tension:cable), part d. The values\n",
    "of $\\cos((i+1)\\pi) - \\cos((i+1)\\pi/2)$ can be summarized in the following\n",
    "table:\n",
    "\n",
    "<table border=\"1\">\n",
    "<thead>\n",
    "<tr><th align=\"center\">$i\\hbox{ mod } 4 = 0$</th> <th align=\"center\">$(i-1)\\hbox{ mod } 4 = 0$</th> <th align=\"center\">$(i-2)\\hbox{ mod } 4 = 0$</th> <th align=\"center\">$(i-3)\\hbox{ mod } 4 = 0$</th> </tr>\n",
    "</thead>\n",
    "<tbody>\n",
    "<tr><td align=\"center\">   $-1 -0$                  </td> <td align=\"center\">   $1 - (-1)$                   </td> <td align=\"center\">   $-1 - 0$                     </td> <td align=\"center\">   $1-1$                        </td> </tr>\n",
    "</tbody>\n",
    "</table>\n",
    "The following function computes the approximate solution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sine_solution(x, N):\n",
    "    from numpy import pi, sin\n",
    "    s = 0\n",
    "    u = [] # u[i] is the solution for N=i\n",
    "    for i in range(N+1):\n",
    "        if i % 4 == 0:\n",
    "            cos_min_cos = -1\n",
    "        elif (i-1) % 4 == 0:\n",
    "            cos_min_cos = 2\n",
    "        elif (i-2) % 4 == 0:\n",
    "            cos_min_cos = -1\n",
    "        elif (i-1) % 4 == 0:\n",
    "            cos_min_cos = 0\n",
    "\n",
    "        b_i = 2/(pi*(i+1))*cos_min_cos\n",
    "        A_ii = (i+1)**2*pi**2/4\n",
    "        c_i = b_i/A_ii\n",
    "        s += c_i*sin((i+1)*x*pi/2)\n",
    "        u.append(s.copy())\n",
    "    return u"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The exact solution is a function defined in a piecewise way. Below we make an implementation that works both for array and scalar arguments:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_solution(x):\n",
    "    if isinstance(x, np.ndarray):\n",
    "        return np.where(x < 1, -1./4*x, 0.5*x**2 - 5./4*x + 0.5)\n",
    "    else:\n",
    "        return -1./4*x if x < 1 else 0.5*x**2 - 5./4*x + 0.5"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can make a plot of the exact solution and approximate solutions for\n",
    "various $N$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_sine_solution():\n",
    "    x = np.linspace(0, 2, 101)\n",
    "    u = sine_solution(x, N=20)\n",
    "    plt.figure()\n",
    "    x = np.linspace(0, 2, 101)\n",
    "    plt.plot(x, exact_solution(x), '--')\n",
    "    N_values = 0, 1, 5\n",
    "    for N in 0, 1, 5, 10:\n",
    "        plt.plot(x, u[N])\n",
    "    plt.legend(['exact'] + ['N=%d' % N for N in N_values])\n",
    "    plt.savefig('tmp2.png');  plt.savefig('tmp2.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:FIGURE: [fig/cable_discont_load_sines.png, width=500 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/cable_discont_load_sines.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "**c)**\n",
    "Solve the problem with P1 finite elements.\n",
    "Plot the solution for $N_e=2,4,8$ elements.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The element matrices and vectors are as for the well-known model\n",
    "problem $u''=1$, except that the element vectors vanish for all\n",
    "elements in $[0,1]$. The following function defines a uniform mesh\n",
    "of P1 elements and runs a finite element algorithm where we use\n",
    "ready-made/known formulas for the element matrix and vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "def P1_solution():\n",
    "    plt.figure()\n",
    "    from fe1D import mesh_uniform, u_glob\n",
    "    N_e_values = [2, 4, 8]\n",
    "    d = 1\n",
    "    legends = []\n",
    "    for N_e in N_e_values:\n",
    "        vertices, cells, dof_map = mesh_uniform(\n",
    "            N_e=N_e, d=d, Omega=[0,2], symbolic=False)\n",
    "        h = vertices[1] - vertices[0]\n",
    "        Ae = 1./h*np.array(\n",
    "            [[1, -1],\n",
    "             [-1, 1]])\n",
    "        N = N_e + 1\n",
    "        A = np.zeros((N, N))\n",
    "        b = np.zeros(N)\n",
    "        for e in range(N_e):\n",
    "            if vertices[e] >= 1:\n",
    "                be = -h/2.*np.array(\n",
    "                    [1, 1])\n",
    "            else:\n",
    "                be = h/2.*np.array(\n",
    "                    [0, 0])\n",
    "            for r in range(d+1):\n",
    "                for s in range(d+1):\n",
    "                    A[dof_map[e][r], dof_map[e][s]] += Ae[r,s]\n",
    "                b[dof_map[e][r]] += be[r]\n",
    "        # Enforce boundary conditions\n",
    "        A[0,:] = 0; A[0,0] = 1; b[0] = 0\n",
    "        A[-1,:] = 0; A[-1,-1] = 1; b[-1] = 0\n",
    "        c = np.linalg.solve(A, b)\n",
    "\n",
    "        # Plot solution\n",
    "        print(('c:', c))\n",
    "        print(('vertices:', vertices))\n",
    "        print(('cells:', cells))\n",
    "        print(('len(cells):', len(cells)))\n",
    "        print(('dof_map:', dof_map))\n",
    "        xc, u, nodes = u_glob(c, vertices, cells, dof_map)\n",
    "        plt.plot(xc, u)\n",
    "        legends.append('$N_e=%d$' % N_e)\n",
    "    plt.plot(xc, exact_solution(xc), '--')\n",
    "    legends.append('exact')\n",
    "    plt.legend(legends, loc='lower left')\n",
    "    plt.savefig('tmp3.png'); plt.savefig('tmp3.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- dom:FIGURE: [fig/cable_discont_load_P1.png, width=500 frac=0.8] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/cable_discont_load_P1.png\" width=500>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "Filename: `cable_discont_load`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 8: Compute with a non-uniform mesh\n",
    "<div id=\"fem:deq:exer:1D:mesh:nonuniform\"></div>\n",
    "\n",
    "\n",
    "**a)**\n",
    "Derive the linear system for the problem $-u''=2$ on $[0,1]$, with\n",
    "$u(0)=0$ and $u(1)=1$, using P1 elements and a *non-uniform* mesh. The\n",
    "vertices have coordinates $x_{0}=0 < x_{1} <\\cdots <\n",
    "x_{N_n-1}=1$, and the length of cell number $e$ is $h_e = x_{e+1}\n",
    "-x_{e}$.\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The element matrix and vector for this problem is given by\n",
    "([71](#fem:deq:1D:ex1:Ab:elm)). The change in this exercise is that\n",
    "$h$ is not a constant element length, but varying with the element\n",
    "number $e$. We therefore write"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\tilde A^{(e)} =\\frac{1}{h_e}\\left(\\begin{array}{rr}\n",
    "1 & -1\\\\ \n",
    "-1 & 1\n",
    "\\end{array}\\right),\\quad\n",
    "\\tilde b^{(e)} = h_e\\left(\\begin{array}{c}\n",
    "1\\\\ \n",
    "1\n",
    "\\end{array}\\right){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assembling such element matrices yields"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\n",
    "\\begin{array}{ccccccccc}\n",
    "h_0^{-1} & -h_0^{-1} & 0 &\\cdots & \\cdots & \\cdots & \\cdots & \\cdots & 0 \\\\ \n",
    "-h_0^{-1} & h_0^{-1}+h_1^{-1} & -h_1^{-1} & \\ddots &   & &  & &  \\vdots \\\\ \n",
    "0 & -h_1^{-1} & h_1^{-1} + h_2^{-1} & -h_2^{-1} &\n",
    "\\ddots & &  &  & \\vdots \\\\ \n",
    "\\vdots & \\ddots &  & \\ddots & \\ddots & 0 &  & & \\vdots \\\\ \n",
    "\\vdots &  & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\ \n",
    "\\vdots & &  & 0 & -h_{i-1}^{-1} & h_{i-1}^{-1} + h_i^{-1} & -h_i^{-1} & \\ddots & \\vdots \\\\ \n",
    "\\vdots & & &  & \\ddots & \\ddots & \\ddots &\\ddots  & 0 \\\\ \n",
    "\\vdots & & & &  &\\ddots  & \\ddots &\\ddots  & -h_{N_e}^{-1} \\\\ \n",
    "0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots  & 0 & -h_{N_e}^{-1} & h_{N_e}^{-1}\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The element vectors assemble to"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\left(\n",
    "\\begin{array}{c}\n",
    "h_0 \\\\ \n",
    "h_0 + h_1\\\\ \n",
    "\\vdots\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots \\\\ \n",
    "h_{i-1} + h_i\\\\ \n",
    "\\vdots \\\\ \n",
    "\\vdots\\\\ \n",
    "h_{N_e}\n",
    "\\end{array}\n",
    "\\right)\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "\n",
    "**b)**\n",
    "It is of interest to compare the discrete equations for the finite\n",
    "element method in a non-uniform mesh with the corresponding discrete\n",
    "equations arising from a finite difference method. Go through the\n",
    "derivation of the finite difference formula $u''(x_i) \\approx [D_x D_x\n",
    "u]_i$ and modify it to find a natural discretization of $u''(x_i)$ on\n",
    "a non-uniform mesh. Compare the finite element and difference\n",
    "discretizations\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "Using the definition of the centered, 2nd-order finite difference approximation\n",
    "to $u''$ we can set up"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[D_xD_x u]_i = [D_x(D_x u)]_i = \\frac{\\frac{u_{i+1}-u_i}{x_{i+1}-x_i} -\n",
    "\\frac{u_{i}-u_{i-1}}{x_{i}-x_{i-1}}}{x_{i+1/2} - x_{i-1/2}}{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now,"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "x_{i+1/2} - x_{i-1/2} = \\frac{1}{2} (x_i - x_{i-1}) + \\frac{1}{2}(x_{i+1}-x_i) = \\frac{1}{2}(x_{i+1}-x_{i-1}){\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We then get the difference equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "u''(x_i)\\approx \\frac{2}{h_i + h_{i-1}}\\left(\n",
    "\\frac{u_{i+1}-u_i}{h_{i}} - \\frac{u_{i}-u_{i-1}}{h_{i-1}}\\right) = 2{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The factor 2 on either side cancels.\n",
    "\n",
    "Looking at the finite element equations in a), the equation for a general\n",
    "row $i$ reads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{1}{h_{i-1}}c_{i-1} - (\\frac{1}{h_{i-1}} + \\frac{1}{h_{i}})c_i\n",
    "+ \\frac{1}{h_{i}}c_{i+1} = h_{i-1} + h_i{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Replacing $c_i$ by $u_i$ (assuming we keep unknowns at all nodes) and\n",
    "rearranging gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\frac{1}{h_{i-1}}(u_i - u_{i-1}) + \\frac{1}{h_{i}}(u_{i+1}-u_i)\n",
    " = h_{i-1} + h_i{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dividing by the right-hand side gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\frac{1}{h_{i-1} + h_i}\\left(\\frac{1}{h_{i-1}}(u_i - u_{i-1}) - \\frac{1}{h_{i}}(u_{i+1}-u_i)\\right) = 1{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is the same difference equation as we have in the finite difference\n",
    "method.\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "\n",
    "Filename: `nonuniform_P1`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Problem 9: Solve a 1D finite element problem by hand\n",
    "<div id=\"fem:deq:exer:1D:gen:problem1\"></div>\n",
    "\n",
    "The following scaled 1D problem is a very simple, yet relevant, model\n",
    "for convective transport in fluids:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto47\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "u' = \\epsilon u'' ,\\quad u(0)=0,\\ u(1)=1,\\ x\\in [0,1]\n",
    "{\\thinspace .}\n",
    "\\label{_auto47} \\tag{119}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)**\n",
    "Find the analytical solution to this problem.\n",
    "(Introduce $w=u'$, solve the first-order differential equation for $w(x)$,\n",
    "and integrate once more.)\n",
    "\n",
    "**b)**\n",
    "Derive the variational form of this problem.\n",
    "\n",
    "**c)**\n",
    "Introduce a finite element mesh with uniform partitioning.\n",
    "Use P1 elements and compute the element matrix and vector for\n",
    "a general element.\n",
    "\n",
    "**d)**\n",
    "Incorporate the boundary conditions and\n",
    "assemble the element contributions.\n",
    "\n",
    "**e)**\n",
    "Identify the resulting linear system as a finite difference discretization\n",
    "of the differential equation using"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "[D_{2x}u = \\epsilon D_xD_x u]_i {\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**f)**\n",
    "Compute the numerical solution and plot it together with the exact solution\n",
    "for a mesh with 20 elements and\n",
    "$\\epsilon=10, 1, 0.1, 0.01$.\n",
    "\n",
    "Filename: `convdiff1D_P1`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 10: Investigate exact finite element solutions\n",
    "<div id=\"fem:deq:exer:1D:exact_numerics\"></div>\n",
    "\n",
    "Consider"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-u''(x)=x^m,\\quad x\\in (0,L),\\quad u'(0)=C,\\ u(L)=D,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where $m\\geq 0$ is an integer, and $L$, $C$, and $D$ are given numbers.\n",
    "Utilize a mesh with two (non-uniform) elements: $\\Omega^{(0)}=[0,3]$ and\n",
    "$\\Omega^{(0)}=[3,4]$.\n",
    "Plot the exact solution and the finite element solution for\n",
    "polynomial degree $d=1,2,3,4$ and $m=0, 1, 2, 3, 4$. Find values of $d$ and $m$\n",
    "that make the finite element solution exact at the nodes in the mesh.\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint.**\n",
    "Use the `mesh_uniform`, `finite_element1D`, and `u_glob2` functions\n",
    "from the `fe1D.py` module.\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "\n",
    "\n",
    "<!-- --- begin solution of exercise --- -->\n",
    "**Solution.**\n",
    "The `model2` function from the section [Simple model problems and their solutions](#fem:deq:1D:models:simple)\n",
    "can find the exact solution by `model2(x**m, L, C, D)`.\n",
    "We fix, for simplicity, the values of $L$, $C$, and $D$ as\n",
    "$L=4$, $C=5$, and $D=2$.\n",
    "After calculating a symbolic solution, we can convert the expression\n",
    "to a Python function with `sympy.lambdify`.\n",
    "For each $d$ value we then create a uniform mesh and displace the\n",
    "vertex with number 1 to the value 3.\n",
    "The various functions for specifying the element matrix and vector\n",
    "entries are as given in the section [Utilizing a sparse matrix](#fem:deq:1D:code:fe_sparse),\n",
    "since the model problem is the same. Our code then becomes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "from u_xx_f_sympy import model2, x\n",
    "import sympy as sym\n",
    "import numpy as np\n",
    "from fe1D import finite_element1D, mesh_uniform, u_glob\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "C = 5\n",
    "D = 2\n",
    "L = 4\n",
    "\n",
    "m_values = [0, 1, 2, 3, 4]\n",
    "d_values = [1, 2, 3, 4]\n",
    "for m in m_values:\n",
    "    u = model2(x**m, L, C, D)\n",
    "    print(('\\nm=%d, u: %s' % (m, u)))\n",
    "    u_exact = sym.lambdify([x], u)\n",
    "\n",
    "    for d in d_values:\n",
    "        vertices, cells, dof_map = mesh_uniform(\n",
    "            N_e=2, d=d, Omega=[0,L], symbolic=False)\n",
    "        vertices[1] = 3  # displace vertex\n",
    "        essbc = {}\n",
    "        essbc[dof_map[-1][-1]] = D\n",
    "\n",
    "        c, A, b, timing = finite_element1D(\n",
    "            vertices, cells, dof_map,\n",
    "            essbc,\n",
    "            ilhs=lambda e, phi, r, s, X, x, h:\n",
    "            phi[1][r](X, h)*phi[1][s](X, h),\n",
    "            irhs=lambda e, phi, r, X, x, h:\n",
    "            x**m*phi[0][r](X),\n",
    "            blhs=lambda e, phi, r, s, X, x, h: 0,\n",
    "            brhs=lambda e, phi, r, X, x, h:\n",
    "            -C*phi[0][r](-1) if e == 0 else 0,\n",
    "            intrule='GaussLegendre')\n",
    "\n",
    "        # Visualize\n",
    "        # (Recall that x is a symbol, use xc for coordinates)\n",
    "        xc, u, nodes = u_glob(c, vertices, cells, dof_map)\n",
    "        u_e = u_exact(xc)\n",
    "        print(('Max diff at nodes, d=%d:' % d, \\\n",
    "              np.abs(u_exact(nodes) - c).max()))\n",
    "        plt.figure()\n",
    "        plt.plot(xc, u, 'b-', xc, u_e, 'r--')\n",
    "        plt.legend(['finite elements, d=%d' %d, 'exact'],\n",
    "                   loc='lower left')\n",
    "        figname = 'tmp_%d_%d' % (m, d)\n",
    "        plt.savefig(figname + '.png'); plt.savefig(figname + '.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First we look at the numerical solution at the nodes:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "        m=0, u: -x**2/2 + 5*x - 10\n",
    "        Max diff at nodes, d=1: 2.22044604925e-16\n",
    "        Max diff at nodes, d=2: 3.5527136788e-15\n",
    "        Max diff at nodes, d=3: 1.7763568394e-15\n",
    "        Max diff at nodes, d=4: 2.46913600677e-13\n",
    "        \n",
    "        m=1, u: -x**3/6 + 5*x - 22/3\n",
    "        Max diff at nodes, d=1: 8.881784197e-16\n",
    "        Max diff at nodes, d=2: 1.7763568394e-15\n",
    "        Max diff at nodes, d=3: 7.9936057773e-15\n",
    "        Max diff at nodes, d=4: 3.01092484278e-13\n",
    "        \n",
    "        m=2, u: -x**4/12 + 5*x + 10/3\n",
    "        Max diff at nodes, d=1: 3.10862446895e-15\n",
    "        Max diff at nodes, d=2: 0.084375\n",
    "        Max diff at nodes, d=3: 0.0333333333333\n",
    "        Max diff at nodes, d=4: 5.20472553944e-13\n",
    "        \n",
    "        m=3, u: -x**5/20 + 5*x + 166/5\n",
    "        Max diff at nodes, d=1: 1.35555555556\n",
    "        Max diff at nodes, d=2: 0.3796875\n",
    "        Max diff at nodes, d=3: 0.185714285714\n",
    "        Max diff at nodes, d=4: 0.0254255022334\n",
    "        \n",
    "        m=4, u: -x**6/30 + 5*x + 1778/15\n",
    "        Max diff at nodes, d=1: 4.8\n",
    "        Max diff at nodes, d=2: 1.4428125\n",
    "        Max diff at nodes, d=3: 0.719047619047\n",
    "        Max diff at nodes, d=4: 0.16865583147\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that all elements are capable of computing the exact values\n",
    "at the nodes for $m=0$ and $m=1$. With $m=0$, the solution is quadratic\n",
    "in $x$, and P2, P3, and P4 will be exact. It is more of a surprise that\n",
    "also the P1 elements are exact in this case.\n",
    "A peculiar feature is that P1\n",
    "elements are also exact at the nodes $m=2$, but not P2 and P3 elements\n",
    "(the solution goes like $x^4$ so it is not surprising that P2 and P3\n",
    "elements give a numerical error also at the nodes).\n",
    "Clearly, P4 elements produce the exact solution for $m=4$ since $u$\n",
    "is a polynomial of degree 4. For larger $m$ values we have\n",
    "discrepancy between the numerical and exact values at the nodes.\n",
    "\n",
    "% for m in [0, 1, 2, 3, 4]:\n",
    "\n",
    "**Plots for m=${m}.**\n",
    "\n",
    "<!-- dom:FIGURE: [fig/u_xx_xm${m}_P1to4, width=800 frac=1] -->\n",
    "<!-- begin figure -->\n",
    "\n",
    "<p></p>\n",
    "<img src=\"fig/u_xx_xm${m}_P1to4\" width=800>\n",
    "\n",
    "<!-- end figure -->\n",
    "\n",
    "% endfor\n",
    "\n",
    "<!-- --- end solution of exercise --- -->\n",
    "Filename: `u_xx_xm_P1to4`.\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- Could have shooting method as a project -->\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 11: Compare finite elements and differences for a radially symmetric Poisson equation\n",
    "<div id=\"fem:deq:exer:1D:Poisson:polar\"></div>\n",
    "\n",
    "We consider the Poisson problem in a disk with radius $R$ with\n",
    "Dirichlet conditions at the boundary.\n",
    "Given that the solution is radially symmetric and hence dependent only on\n",
    "the radial coordinate ($r=\\sqrt{x^2+y^2}$), we can reduce the problem\n",
    "to a 1D Poisson equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:exer:1D:Poisson:polar:eq\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "-\\frac{1}{r}\\frac{d}{dr}\\left( r\\frac{du}{dr}\\right) = f(r),\\quad r\\in (0,R),\\ \n",
    "u'(0)=0,\\ u(R)=U_R\n",
    "{\\thinspace .}\n",
    "\\label{fem:deq:exer:1D:Poisson:polar:eq} \\tag{120}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**a)**\n",
    "Derive a variational form of ([120](#fem:deq:exer:1D:Poisson:polar:eq))\n",
    "by integrating over the whole disk, or posed equivalently: use\n",
    "a weighting function $2\\pi r v(r)$ and integrate $r$ from $0$ to $R$.\n",
    "\n",
    "**b)**\n",
    "Use a uniform mesh partition with P1 elements and show what the\n",
    "resulting set of equations becomes. Integrate the matrix entries\n",
    "exact by hand, but use a Trapezoidal rule to integrate the $f$ term.\n",
    "\n",
    "**c)**\n",
    "Explain that an intuitive\n",
    "finite difference method applied to ([120](#fem:deq:exer:1D:Poisson:polar:eq))\n",
    "gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\frac{1}{r_i}\\frac{1}{h^2}\\left( r_{i+\\frac{1}{2}}(u_{i+1}-u_i) -\n",
    "r_{i-\\frac{1}{2}}(u_{i}-u_{i-1})\\right) = f_i,\\quad i=rh\n",
    "{\\thinspace .}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For $i=0$ the factor $1/r_i$ seemingly becomes problematic. One must always\n",
    "have $u'(0)=0$, because of the radial symmetry, which implies\n",
    "$u_{-1}=u_1$, if we allow introduction of a fictitious value $u_{-1}$.\n",
    "Using this $u_{-1}$ in the difference equation for $i=0$ gives"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "\\begin{align*}\n",
    "&\\frac{1}{r_0}\\frac{1}{h^2}\\left( r_{\\frac{1}{2}}(u_{1}-u_0) -\n",
    "r_{-\\frac{1}{2}}(u_{0}-u_{1})\\right) = \\\\ \n",
    "& \\qquad\n",
    "\\frac{1}{r_0}\\frac{1}{2h^2}\\left( (r_0 + r_1)(u_{1}-u_0) -\n",
    "(r_{-1} + r_0)(u_{0}-u_{1})\\right) \\approx\n",
    "2(u_1-u_0),\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "if we use $r_{-1}+r_1\\approx 2r_0$.\n",
    "\n",
    "Set up the complete set of equations for the finite difference method\n",
    "and compare to the finite element method in case a Trapezoidal rule\n",
    "is used to integrate the $f$ term in the latter method.\n",
    "\n",
    "Filename: `radial_Poisson1D_P1`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 12: Compute with variable coefficients and P1 elements by hand\n",
    "<div id=\"fem:deq:exer:1D:gen:problem2\"></div>\n",
    "\n",
    "Consider the problem"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"fem:deq:1D:model4\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation}\n",
    "-\\frac{d}{dx}\\left( {\\alpha}(x)\\frac{du}{dx}\\right) + \\gamma u = f(x),\n",
    "\\quad x\\in\\Omega=[0,L],\\quad u(0)={\\alpha},\\ u'(L)=\\beta{\\thinspace .}\n",
    "\\label{fem:deq:1D:model4} \\tag{121}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We choose ${\\alpha}(x)=1+x^2$. Then"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<!-- Equation labels as ordinary links -->\n",
    "<div id=\"_auto48\"></div>\n",
    "\n",
    "$$\n",
    "\\begin{equation} u(x) = {\\alpha} + \\beta(1+L^2)\\tan^{-1}(x),\n",
    "\\label{_auto48} \\tag{122}\n",
    "\\end{equation}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "is an exact solution if $f(x) = \\gamma u$.\n",
    "\n",
    "Derive a variational formulation and compute general expressions for the\n",
    "element matrix and vector in an arbitrary element, using P1 elements\n",
    "and a uniform partitioning of $[0,L]$. The right-hand side\n",
    "integral is challenging and can be computed by a numerical integration\n",
    "rule. The Trapezoidal rule ([fem:approx:fe:numint1:trapez](#fem:approx:fe:numint1:trapez))\n",
    "gives particularly simple expressions.\n",
    "Filename: `atan1D_P1`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 13: Solve a 2D Poisson equation using polynomials and sines\n",
    "<div id=\"fem:deq:exer:2D:torsion:xy:sin\"></div>\n",
    "\n",
    "The classical problem of applying a torque to the ends of a rod\n",
    "can be modeled by a Poisson equation defined in the cross section $\\Omega$:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$\n",
    "-\\nabla^2 u = 2,\\quad (x,y)\\in\\Omega,\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $u=0$ on $\\partial\\Omega$. Exactly the same problem arises for\n",
    "the deflection of a membrane with shape $\\Omega$ under a constant load.\n",
    "\n",
    "For a circular cross section one can readily\n",
    "find an analytical solution. For a rectangular cross section the analytical\n",
    "approach ends up with a sine series. The idea in this exercise is to\n",
    "use a single basis function to obtain an approximate answer.\n",
    "\n",
    "We assume for simplicity that the cross section is the unit square:\n",
    "$\\Omega = [0,1]\\times [0,1]$.\n",
    "\n",
    "\n",
    "**a)**\n",
    "We consider the basis\n",
    "${\\psi}_{p,q}(x,y) = \\sin((p+1)\\pi x)\\sin (q\\pi y)$, $p,q=0,\\ldots,n$.\n",
    "These basis functions fulfill the Dirichlet condition.\n",
    "Use a Galerkin method and $n=0$.\n",
    "\n",
    "**b)**\n",
    "The basis function involving sine functions are orthogonal.\n",
    "Use this property in the Galerkin method\n",
    "to derive the coefficients $c_{p,q}$ in a\n",
    "formula $u=\\sum_p\\sum_q c_{p,q}{\\psi}_{p,q}(x,y)$.\n",
    "\n",
    "**c)**\n",
    "Another possible basis is\n",
    "${\\psi}_i(x,y) = (x(1-x)y(1-y))^{i+1}$, $i=0,\\ldots,N$.\n",
    "Use the Galerkin method to compute the solution for $N=0$.\n",
    "Which choice of a single basis function is best,\n",
    "$u\\sim x(1-x)y(1-y)$ or $u\\sim \\sin(\\pi x)\\sin(\\pi y)$?\n",
    "In order to answer the question,\n",
    "it is necessary to search the web or the literature for an accurate\n",
    "estimate of the maximum $u$ value at $x=y=1/2$.\n",
    "\n",
    "Filename: `torsion_sin_xy`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 14: Solve a 3D Laplace problem with FEniCS\n",
    "<div id=\"fem:fenics:borehole:exer:3D\"></div>\n",
    "\n",
    "Solve the problem in the section [Implementation in 2D and 3D via FEniCS](#fem:varform:fenics) in 3D.\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint.**\n",
    "Use `BoxMesh` as starting point. Consult the FEniCS tutorial [[ftut1]](#ftut1)\n",
    "if necessary.\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "Filename: `borehole_fenics3D`.\n",
    "\n",
    "<!-- --- end exercise --- -->\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "<!-- --- begin exercise --- -->\n",
    "\n",
    "## Exercise 15: Solve a 1D Laplace problem with FEniCS\n",
    "<div id=\"fem:fenics:borehole:exer:1D\"></div>\n",
    "\n",
    "Solve the problem in the section [Implementation in 2D and 3D via FEniCS](#fem:varform:fenics) in 1D, using\n",
    "the radial formulation $(ru')'=0$.\n",
    "\n",
    "<!-- --- begin hint in exercise --- -->\n",
    "\n",
    "**Hint.**\n",
    "Use `IntervalMesh` for generating the mesh and introduce a stretching if\n",
    "desired. Consult the FEniCS tutorial [[ftut1]](#ftut1)\n",
    "so you can extract the solution in an array\n",
    "and make your own curve plot of it.\n",
    "\n",
    "This problem can be solved without markers, see the section on\n",
    "multiple Dirichlet conditions in the tutorial [[ftut1]](#ftut1) (it is even\n",
    "easier to solve by saying that the solution at the boundary\n",
    "obeys a linear function from $u_a$ to $u_b$ and use this as the only\n",
    "Dirichlet condition).\n",
    "\n",
    "<!-- --- end hint in exercise --- -->\n",
    "Filename: `borehole_fenics1D`.\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
}