{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- 2DO -->\n",
"<!-- Must illustrate how to make weak form of continuous problem and -->\n",
"<!-- discretize. Do that in time-dependent problems too. -->\n",
"\n",
"<!-- Maybe <,> \\langle, \\rangle as inner product -->\n",
"\n",
"<!-- Must say something about error estimates! -->\n",
"\n",
"\n",
"<!-- no (au')', have (\\alpha u')' - it solves all the problems with a and a(.,.) -->\n",
"\n",
"# Variational formulations with global basis functions\n",
"<div id=\"ch:varform:global\"></div>\n",
"\n",
"The finite element method is a very flexible approach for solving partial\n",
"differential equations. Its two most attractive features are the ease\n",
"of handling domains of complex shape in two and three dimensions and\n",
"the variety of polynomials (with different properties and orders) \n",
"that are available.\n",
"The latter feature typically leads to errors proportional to\n",
"$h^{d+1}$, where $h$ is the element length and $d$ is the polynomial\n",
"degree. When the solution is sufficiently smooth, the ability to use\n",
"larger $d$ creates methods that are much more computationally efficient\n",
"than standard finite difference methods (and equally efficient finite\n",
"difference methods are technically much harder to construct).\n",
"\n",
"However, before we attack finite element methods, with localized basis\n",
"functions, it can be easier from a pedagogical point of view to study\n",
"approximations by global functions because the mathematics in this\n",
"case gets simpler.\n",
"\n",
"# Basic principles for approximating differential equations\n",
"<div id=\"fem:deq:1D:principles\"></div>\n",
"\n",
"The finite element method is usually applied for discretization in\n",
"space, and therefore spatial problems will be our focus in the coming\n",
"sections. Extensions to time-dependent problems usually employs\n",
"finite difference approximations in time.\n",
"\n",
"The coming sections address at how global basis functions and the least\n",
"squares, Galerkin, and collocation principles can be used to solve\n",
"differential equations.\n",
"\n",
"## Differential equation models\n",
"<div id=\"fem:deq:1D:models\"></div>\n",
"\n",
"Let us consider an abstract differential equation for a function $u(x)$ of\n",
"one variable, written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\mathcal{L}(u) = 0,\\quad x\\in\\Omega{\\thinspace .} \\label{_auto1} \\tag{1}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here are a few examples on possible choices of $\\mathcal{L}(u)$, of\n",
"increasing complexity:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:L1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\mathcal{L}(u) = \\frac{d^2u}{dx^2} - f(x),\n",
"\\label{fem:deq:1D:L1} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:L2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\mathcal{L}(u) = \\frac{d}{dx}\\left({\\alpha}(x)\\frac{du}{dx}\\right) + f(x),\n",
"\\label{fem:deq:1D:L2} \\tag{3}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:L3\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\mathcal{L}(u) = \\frac{d}{dx}\\left({\\alpha}(u)\\frac{du}{dx}\\right) - au + f(x),\n",
"\\label{fem:deq:1D:L3} \\tag{4}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:L4\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\mathcal{L}(u) = \\frac{d}{dx}\\left({\\alpha}(u)\\frac{du}{dx}\\right) + f(u,x)\n",
"\\label{fem:deq:1D:L4} \\tag{5}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Both ${\\alpha}(x)$ and $f(x)$ are considered as specified functions,\n",
"while $a$ is a prescribed parameter. Differential equations\n",
"corresponding to ([2](#fem:deq:1D:L1))-([3](#fem:deq:1D:L2)) arise in\n",
"diffusion phenomena, such as stationary (time-independent)\n",
"transport of heat in solids and\n",
"flow of viscous fluids between flat plates. The form\n",
"([4](#fem:deq:1D:L3)) arises when transient diffusion or wave\n",
"phenomena are discretized in time by finite differences. The equation\n",
"([5](#fem:deq:1D:L4)) appears in chemical models when diffusion of a\n",
"substance is combined with chemical reactions. Also in biology,\n",
"([5](#fem:deq:1D:L4)) plays an important role, both for spreading of\n",
"species and in models involving generation and\n",
"propagation of electrical signals.\n",
"\n",
"Let $\\Omega =[0,L]$ be the domain in one space dimension.\n",
"In addition to the differential equation, $u$ must fulfill\n",
"boundary conditions at the boundaries of the domain, $x=0$ and $x=L$.\n",
"When $\\mathcal{L}$ contains up to second-order derivatives, as in the\n",
"examples above, we need one boundary condition at each of\n",
"the (two) boundary points, here abstractly specified as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\mathcal{B}_0(u)=0,\\ x=0,\\quad \\mathcal{B}_1(u)=0,\\ x=L\n",
"\\label{_auto2} \\tag{6}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are three common choices of boundary conditions:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto3\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\mathcal{B}_i(u) = u - g,\\quad \\hbox{Dirichlet condition}\n",
"\\label{_auto3} \\tag{7}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto4\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\mathcal{B}_i(u) = -{\\alpha} \\frac{du}{dx} - g,\\quad \\hbox{Neumann condition}\n",
"\\label{_auto4} \\tag{8}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto5\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\mathcal{B}_i(u) = -{\\alpha} \\frac{du}{dx} - H(u-g),\\quad \\hbox{Robin condition}\n",
"\\label{_auto5} \\tag{9}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, $g$ and $H$ are specified quantities.\n",
"\n",
"From now on we shall use ${u_{\\small\\mbox{e}}}(x)$ as symbol for the *exact* solution,\n",
"fulfilling"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto6\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\mathcal{L}({u_{\\small\\mbox{e}}})=0,\\quad x\\in\\Omega,\n",
"\\label{_auto6} \\tag{10}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"while $u(x)$ is our notation for an *approximate* solution of the differential\n",
"equation.\n",
"\n",
"**Remark on notation.**\n",
"\n",
"In the literature about the finite element method,\n",
"it is common to use $u$ as the exact solution and $u_h$ as the\n",
"approximate solution, where $h$ is a discretization parameter. However,\n",
"the vast part of the present text is about the approximate solutions,\n",
"and having a subscript $h$ attached all the time\n",
"is cumbersome. Of equal importance is the close correspondence between\n",
"implementation and mathematics that we strive to achieve in this text:\n",
"when it is natural to use `u` and not `u_h` in\n",
"code, we let the mathematical notation be dictated by the code's\n",
"preferred notation. In the relatively few cases where we need to work\n",
"with the exact solution of the PDE problem we call it ${u_{\\small\\mbox{e}}}$ in\n",
"mathematics and `u_e` in the code (the function for computing\n",
"`u_e` is named `u_exact`).\n",
"<!-- After all, it is the powerful computer implementations -->\n",
"<!-- of the finite element method that justifies studying the mathematical -->\n",
"<!-- formulation and aspects of the method. -->\n",
"\n",
"\n",
"\n",
"## Simple model problems and their solutions\n",
"<div id=\"fem:deq:1D:models:simple\"></div>\n",
"\n",
"A common model problem used much in the forthcoming examples is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:model1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"-u''(x) = f(x),\\quad x\\in\\Omega=[0,L],\\quad u(0)=0,\\ u(L)=D\n",
"{\\thinspace .}\n",
"\\label{fem:deq:1D:model1} \\tag{11}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A closely related problem with a different boundary condition at\n",
"$x=0$ reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:model2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"-u''(x) = f(x),\\quad x\\in\\Omega=[0,L],\\quad u'(0)=C,\\ u(L)=D{\\thinspace .}\n",
"\\label{fem:deq:1D:model2} \\tag{12}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A third variant has a variable coefficient,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:model3\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"-({\\alpha}(x)u'(x))' = f(x),\\quad x\\in\\Omega=[0,L],\\quad u'(0)=C,\\ u(L)=D{\\thinspace .}\n",
"\\label{fem:deq:1D:model3} \\tag{13}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution $u$ to the model problem ([11](#fem:deq:1D:model1))\n",
"can be determined as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"u'(x) &= -\\int_0^x f(x) + c_0, \\\\ \n",
"u(x) &= \\int_0^x u'(x) + c_1,\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $c_0$ and $c_1$ are determined by the boundary conditions\n",
"such that $u'(0) = C$ and $u(L) = D$.\n",
"\n",
"Computing the solution is easily done\n",
"using `sympy`. Some common code is defined first:"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import sympy as sym\n",
"x, L, C, D, c_0, c_1, = sym.symbols('x L C D c_0 c_1')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The following function computes the solution\n",
"symbolically for the model problem ([11](#fem:deq:1D:model1)):"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def model1(f, L, D):\n",
" \"\"\"Solve -u'' = f(x), u(0)=0, u(L)=D.\"\"\"\n",
" # Integrate twice\n",
" u_x = - sym.integrate(f, (x, 0, x)) + c_0\n",
" u = sym.integrate(u_x, (x, 0, x)) + c_1\n",
" # Set up 2 equations from the 2 boundary conditions and solve\n",
" # with respect to the integration constants c_0, c_1\n",
" r = sym.solve([u.subs(x, 0)-0, # x=0 condition\n",
" u.subs(x,L)-D], # x=L condition\n",
" [c_0, c_1]) # unknowns\n",
" # Substitute the integration constants in the solution\n",
" u = u.subs(c_0, r[c_0]).subs(c_1, r[c_1])\n",
" u = sym.simplify(sym.expand(u))\n",
" return u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Calling `model1(2, L, D)` results in the solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:model1:sol\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(x) = \\frac{1}{L}x \\left(D + L^{2} - L x\\right)\n",
"\\label{fem:deq:1D:model1:sol} \\tag{14}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The model problem ([12](#fem:deq:1D:model2)) can be solved by"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def model2(f, L, C, D):\n",
" \"\"\"Solve -u'' = f(x), u'(0)=C, u(L)=D.\"\"\"\n",
" u_x = - sym.integrate(f, (x, 0, x)) + c_0\n",
" u = sym.integrate(u_x, (x, 0, x)) + c_1\n",
" r = sym.solve([sym.diff(u,x).subs(x, 0)-C, # x=0 cond.\n",
" u.subs(x,L)-D], # x=L cond.\n",
" [c_0, c_1])\n",
" u = u.subs(c_0, r[c_0]).subs(c_1, r[c_1])\n",
" u = sym.simplify(sym.expand(u))\n",
" return u"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"to yield"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:model2:sol\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(x) = - x^{2} + C x - C L + D + L^{2},\n",
"\\label{fem:deq:1D:model2:sol} \\tag{15}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"if $f(x)=2$. Model ([13](#fem:deq:1D:model3)) requires a bit more involved\n",
"code,"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"('model1:', x*(D + L*(L - x))/L, 0, D)\n",
"\\frac{x \\left(D + L \\left(L - x\\right)\\right)}{L}\n",
"('model2:', -C*L + C*x + D + L**2 - x**2, C, D)\n",
"- C L + C x + D + L^{2} - x^{2}\n",
"('model3:', (C*atan(L) - C*atan(x) + D*atan(x))/atan(L), C, D)\n",
"\\frac{C \\operatorname{atan}{\\left(L \\right)} - C \\operatorname{atan}{\\left(x \\right)} + D \\operatorname{atan}{\\left(x \\right)}}{\\operatorname{atan}{\\left(L \\right)}}\n"
]
}
],
"source": [
"def model3(f, a, L, C, D):\n",
" \"\"\"Solve -(a*u')' = f(x), u(0)=C, u(L)=D.\"\"\"\n",
" au_x = - sym.integrate(f, (x, 0, x)) + c_0\n",
" u = sym.integrate(au_x/a, (x, 0, x)) + c_1\n",
" r = sym.solve([u.subs(x, 0)-C,\n",
" u.subs(x,L)-D],\n",
" [c_0, c_1])\n",
" u = u.subs(c_0, r[c_0]).subs(c_1, r[c_1])\n",
" u = sym.simplify(sym.expand(u))\n",
" return u\n",
"\n",
"\n",
"def demo():\n",
" f = 2\n",
" u = model1(f, L, D)\n",
" print(('model1:', u, u.subs(x, 0), u.subs(x, L)))\n",
" print((sym.latex(u, mode='plain')))\n",
" u = model2(f, L, C, D)\n",
" #f = x\n",
" #u = model2(f, L, C, D)\n",
" print(('model2:', u, sym.diff(u, x).subs(x, 0), u.subs(x, L)))\n",
" print((sym.latex(u, mode='plain')))\n",
" u = model3(0, 1+x**2, L, C, D)\n",
" print(('model3:', u, u.subs(x, 0), u.subs(x, L)))\n",
" print((sym.latex(u, mode='plain')))\n",
"\n",
"if __name__ == '__main__':\n",
" demo()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With $f(x)=0$ and ${\\alpha}(x)=1+x^2$ we get"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x) =\n",
"\\frac{C \\tan^{-1}\\left (L \\right ) - C \\tan^{-1}\\left (x \\right ) + D \\tan^{-1}\\left (x \\right )}{\\tan^{-1}\\left (L \\right )}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Forming the residual\n",
"<div id=\"fem:deq:1D:residual:min\"></div>\n",
"\n",
"The fundamental idea is to seek an approximate solution\n",
"$u$ in some space $V$,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"V = \\hbox{span}\\{ {\\psi}_0(x),\\ldots,{\\psi}_N(x)\\},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which means that $u$ can always be expressed as a linear combination\n",
"of the basis functions $\\left\\{ {{\\psi}}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$, with ${\\mathcal{I}_s}$ as\n",
"the index set $\\{0,\\ldots,N\\}$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x) = \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(x){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The coefficients $\\left\\{ {c}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$ are unknowns to be computed.\n",
"\n",
"(Later, we will see that if we specify boundary values of $u$ different\n",
"from zero, we must look for an approximate solution\n",
"$u(x) = B(x) + \\sum_{j} c_j{\\psi}_j(x)$,\n",
"where $\\sum_{j}c_j{\\psi}_j\\in V$ and $B(x)$ is some function for\n",
"incorporating the right boundary values. Because of $B(x)$, $u$ will not\n",
"necessarily lie in $V$. This modification does not imply any difficulties.)\n",
"\n",
"We need principles for deriving $N+1$ equations to determine the\n",
"$N+1$ unknowns $\\left\\{ {c}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$.\n",
"When approximating a given function $f$ by $u=\\sum_jc_j{\\varphi}_j$,\n",
"a key idea is to minimize the square norm of the\n",
"approximation error $e=u-f$ or (equivalently) demand that $e$ is\n",
"orthogonal to $V$. Working with $e$ is not so useful here since\n",
"the approximation error in our case is $e={u_{\\small\\mbox{e}}} - u$ and ${u_{\\small\\mbox{e}}}$ is\n",
"unknown. The only general indicator we have on the quality of the approximate\n",
"solution is to what degree $u$ fulfills the differential equation.\n",
"Inserting $u=\\sum_j c_j {\\psi}_j$ into $\\mathcal{L}(u)$ reveals that the\n",
"result is not zero, because $u$ in general is an approximation and not identical to ${u_{\\small\\mbox{e}}}$.\n",
"The nonzero result,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto7\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"R = \\mathcal{L}(u) = \\mathcal{L}(\\sum_j c_j {\\psi}_j),\n",
"\\label{_auto7} \\tag{16}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"is called the *residual* and measures the\n",
"error in fulfilling the governing equation.\n",
"\n",
"Various principles for determining $\\left\\{ {c}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$ try to minimize\n",
"$R$ in some sense. Note that $R$ varies with $x$ and\n",
"the $\\left\\{ {c}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$ parameters. We may write this dependence\n",
"explicitly as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto8\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"R = R(x; c_0, \\ldots, c_N){\\thinspace .} \\label{_auto8} \\tag{17}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Below, we present three principles for making $R$ small:\n",
"a least squares method, a projection or Galerkin method, and\n",
"a collocation or interpolation method.\n",
"\n",
"## The least squares method\n",
"\n",
"The least-squares method aims to find $\\left\\{ {c}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$ such that\n",
"the square norm of the residual"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto9\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"||R|| = (R, R) = \\int_{\\Omega} R^2 {\\, \\mathrm{d}x}\n",
"\\label{_auto9} \\tag{18}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"is minimized. By introducing\n",
"an inner product of two functions $f$ and $g$\n",
"on $\\Omega$ as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto10\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(f,g) = \\int_{\\Omega} f(x)g(x) {\\, \\mathrm{d}x},\n",
"\\label{_auto10} \\tag{19}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"the least-squares method can be defined as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto11\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\min_{c_0,\\ldots,c_N} E = (R,R){\\thinspace .} \\label{_auto11} \\tag{20}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Differentiating with respect to the free parameters $\\left\\{ {c}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$\n",
"gives the $N+1$ equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:LS:eq1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{\\Omega} 2R\\frac{\\partial R}{\\partial c_i} {\\, \\mathrm{d}x} = 0\\quad\n",
"\\Leftrightarrow\\quad (R,\\frac{\\partial R}{\\partial c_i})=0,\\quad\n",
"i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"\\label{fem:deq:1D:LS:eq1} \\tag{21}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The Galerkin method\n",
"\n",
"The least-squares\n",
"principle is equivalent to demanding the error to be orthogonal to\n",
"the space $V$ when approximating a function $f$ by $u\\in V$.\n",
"With a differential equation\n",
"we do not know the true error so we must instead require the residual $R$\n",
"to be orthogonal to $V$. This idea implies\n",
"seeking $\\left\\{ {c}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$ such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:Galerkin0\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(R,v)=0,\\quad \\forall v\\in V{\\thinspace .}\n",
"\\label{fem:deq:1D:Galerkin0} \\tag{22}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the Galerkin method for differential equations.\n",
"\n",
"\n",
"<!-- As shown in ([fem:approx:vec:Np1dim:Galerkin](#fem:approx:vec:Np1dim:Galerkin)) and ([fem:approx:vec:Np1dim:Galerkin0](#fem:approx:vec:Np1dim:Galerkin0)), -->\n",
"The above abstract statement can be made concrete by choosing a concrete basis.\n",
"For example, the statement is equivalent to $R$ being orthogonal to the $N+1$\n",
"basis functions $\\{{\\psi}_i\\}$ spanning $V$ (and this is\n",
"the most convenient way to express ([22](#fem:deq:1D:Galerkin0)):"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:Galerkin\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(R,{\\psi}_i)=0,\\quad i\\in{\\mathcal{I}_s},\n",
"\\label{fem:deq:1D:Galerkin} \\tag{23}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"resulting in $N+1$ equations for determining $\\left\\{ {c}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$.\n",
"\n",
"## The method of weighted residuals\n",
"\n",
"\n",
"A generalization of the Galerkin method is to demand that $R$\n",
"is orthogonal to some space $W$, but not necessarily the same\n",
"space as $V$ where we seek the unknown function.\n",
"This generalization is called the *method of weighted residuals*:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:WRM0\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(R,v)=0,\\quad \\forall v\\in W{\\thinspace .}\n",
"\\label{fem:deq:1D:WRM0} \\tag{24}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If $\\{w_0,\\ldots,w_N\\}$ is a basis for $W$, we can equivalently\n",
"express the method of weighted residuals as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:WRM\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(R,w_i)=0,\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"\\label{fem:deq:1D:WRM} \\tag{25}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The result is $N+1$ equations for $\\left\\{ {c}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$.\n",
"\n",
"The least-squares method can also be viewed as a weighted residual\n",
"method with $w_i = \\partial R/\\partial c_i$.\n",
"\n",
"\n",
"**Variational formulation of the continuous problem.**\n",
"\n",
"Statements like ([22](#fem:deq:1D:Galerkin0)), ([23](#fem:deq:1D:Galerkin)),\n",
"([24](#fem:deq:1D:WRM0)), or\n",
"([25](#fem:deq:1D:WRM)))\n",
"are known as\n",
"[weak formulations](https://en.wikipedia.org/wiki/Weak_formulation)\n",
"or *variational formulations*.\n",
"These equations are in this text primarily used for a numerical approximation\n",
"$u\\in V$, where $V$ is a *finite-dimensional* space with dimension\n",
"$N+1$. However, we may also let the exact solution ${u_{\\small\\mbox{e}}}$ fulfill a\n",
"variational formulation $(\\mathcal{L}({u_{\\small\\mbox{e}}}),v)=0$ $\\forall v\\in V$,\n",
"but the exact solution lies in general in a space with infinite\n",
"dimensions (because an infinite number of parameters are needed to\n",
"specify the solution). The variational formulation for ${u_{\\small\\mbox{e}}}$\n",
"in an infinite-dimensional space $V$ is\n",
"a mathematical way of stating the problem and acts as an\n",
"alternative to the usual (strong) formulation of a differential equation with\n",
"initial and/or boundary conditions.\n",
"\n",
"Much of the literature on finite\n",
"element methods takes a differential equation problem and first\n",
"transforms it to a variational formulation in an infinite-dimensional space\n",
"$V$, before searching for an approximate solution in a finite-dimensional\n",
"subspace of $V$. However, we prefer the more intuitive approach with an\n",
"approximate solution $u$ in a finite-dimensional space $V$ inserted in\n",
"the differential equation, and then the resulting residual is demanded to be\n",
"orthogonal to $V$.\n",
"\n",
"\n",
"\n",
"\n",
"**Remark on terminology.**\n",
"\n",
"The terms weak or variational formulations often refer to a statement like\n",
"([22](#fem:deq:1D:Galerkin0)) or ([24](#fem:deq:1D:WRM0))\n",
"after *integration by parts* has been performed (the integration by\n",
"parts technique is\n",
"explained in the section [Integration by parts](#fem:deq:1D:varform)).\n",
"The result after\n",
"integration by parts is what is obtained after taking the *first\n",
"variation* of a minimization problem (see\n",
"the section [Variational problems and minimization of functionals](#fem:deq:1D:optimization)).\n",
"However, in this text we use variational formulation as a common term for\n",
"formulations which, in contrast to the differential equation $R=0$,\n",
"instead demand that an average of $R$ is zero: $(R,v)=0$ for all $v$ in some space.\n",
"\n",
"\n",
"\n",
"## Test and trial functions\n",
"\n",
"\n",
"In the context of the Galerkin method and the method of weighted residuals it is\n",
"common to use the name *trial function* for the approximate $u =\n",
"\\sum_j c_j {\\psi}_j$.\n",
"<!-- Sometimes the functions that spans the space where $u$ lies are also called -->\n",
"<!-- trial functions. -->\n",
"The space containing the trial function is known as the *trial space*.\n",
"The function $v$ entering the orthogonality requirement in\n",
"the Galerkin method and the method of weighted residuals is called\n",
"*test function*, and so are the ${\\psi}_i$ or $w_i$ functions that are\n",
"used as weights in the inner products with the residual. The space\n",
"where the test functions comes from is naturally called the\n",
"*test space*.\n",
"\n",
"We see that in the method of weighted residuals the test and trial spaces\n",
"are different and so are the test and trial functions.\n",
"In the Galerkin method the test and trial spaces are the same (so far).\n",
"<!-- Later in the section [Boundary conditions: specified nonzero value](#fem:deq:1D:essBC) we shall see that boundary -->\n",
"<!-- conditions may lead to a difference between the test and trial spaces -->\n",
"<!-- in the Galerkin method. -->\n",
"\n",
"\n",
"## The collocation method\n",
"\n",
"The idea of the collocation method is to demand that $R$ vanishes\n",
"at $N+1$ selected points $x_{0},\\ldots,x_{N}$ in $\\Omega$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:collocation\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"R(x_{i}; c_0,\\ldots,c_N)=0,\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"\\label{fem:deq:1D:collocation} \\tag{26}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The collocation method can also be viewed as a method of weighted residuals\n",
"with Dirac delta functions as weighting functions.\n",
"Let $\\delta (x-x_{i})$ be the Dirac delta function centered around\n",
"$x=x_{i}$ with the properties that $\\delta (x-x_{i})=0$ for $x\\neq x_{i}$\n",
"and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:Dirac\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{\\Omega} f(x)\\delta (x-x_{i}) {\\, \\mathrm{d}x} =\n",
"f(x_{i}),\\quad x_{i}\\in\\Omega{\\thinspace .}\n",
"\\label{fem:deq:1D:Dirac} \\tag{27}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Intuitively, we may think of $\\delta (x-x_{i})$ as a very peak-shaped\n",
"function around $x=x_{i}$ with an integral $\\int_{-\\infty}^\\infty \\delta(x-x_{i})dx$ that evaluates to unity. Mathematically, it can be shown that\n",
"$\\delta (x-x_{i})$ is the limit of a Gaussian function centered at\n",
"$x=x_{i}$ with a standard deviation that approaches zero.\n",
"Using this latter model, we can roughly visualize delta functions as\n",
"done in [Figure](#fem:deq:1D:fig:Dirac).\n",
"Because of ([27](#fem:deq:1D:Dirac)), we can let $w_i=\\delta(x-x_{i})$\n",
"be weighting functions in the method of weighted residuals,\n",
"and ([25](#fem:deq:1D:WRM)) becomes equivalent to\n",
"([26](#fem:deq:1D:collocation)).\n",
"\n",
"<!-- dom:FIGURE: [fig/delta_func_weight.png, width=400] Approximation of delta functions by narrow Gaussian functions. <div id=\"fem:deq:1D:fig:Dirac\"></div> -->\n",
"<!-- begin figure -->\n",
"<div id=\"fem:deq:1D:fig:Dirac\"></div>\n",
"\n",
"<p>Approximation of delta functions by narrow Gaussian functions.</p>\n",
"<img src=\"fig/delta_func_weight.png\" width=400>\n",
"\n",
"<!-- end figure -->\n",
"\n",
"\n",
"\n",
"### The subdomain collocation method\n",
"\n",
"The idea of this approach is to demand the integral of $R$ to vanish\n",
"over $N+1$ subdomains $\\Omega_i$ of $\\Omega$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto12\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{\\Omega_i} R\\, {\\, \\mathrm{d}x}=0,\\quad i\\in{\\mathcal{I}_s}{\\thinspace .} \\label{_auto12} \\tag{28}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This statement can also be expressed as a weighted residual method"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto13\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{\\Omega} Rw_i\\, {\\, \\mathrm{d}x}=0,\\quad i\\in{\\mathcal{I}_s}, \\label{_auto13} \\tag{29}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $w_i=1$ for $x\\in\\Omega_i$ and $w_i=0$ otherwise.\n",
"\n",
"\n",
"## Examples on using the principles\n",
"<div id=\"fem:deq:1D:ex:sines\"></div>\n",
"\n",
"Let us now apply global basis functions to illustrate the different\n",
"principles for making the residual $R$ small.\n",
"\n",
"### The model problem\n",
"\n",
"We consider the differential equation problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:model1b\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"-u''(x) = f(x),\\quad x\\in\\Omega=[0,L],\\quad u(0)=0,\\ u(L)=0\n",
"{\\thinspace .}\n",
"\\label{fem:deq:1D:model1b} \\tag{30}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Basis functions\n",
"\n",
"Our choice of basis functions ${\\psi}_i$\n",
"for $V$ is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:ex:sines:psi\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"{\\psi}_i(x) = {\\sin\\left((i+1)\\pi\\frac{x}{L}\\right)},\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"\\label{fem:deq:1D:ex:sines:psi} \\tag{31}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"An important property of these functions is that ${\\psi}_i(0)={\\psi}_i(L)=0$,\n",
"which means that the boundary conditions on $u$ are fulfilled:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(0) = \\sum_jc_j{\\psi}_j(0) = 0,\\quad u(L) = \\sum_jc_j{\\psi}_j(L) =0\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Another nice property is that the chosen sine functions\n",
"are orthogonal on $\\Omega$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto14\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int\\limits_0^L {\\sin\\left((i+1)\\pi\\frac{x}{L}\\right)}{\\sin\\left((j+1)\\pi\\frac{x}{L}\\right)}\\, {\\, \\mathrm{d}x} = \\left\\lbrace\n",
"\\begin{array}{ll} \\frac{1}{2} L & i=j \\\\ 0, & i\\neq j\n",
"\\end{array}\\right.\n",
"\\label{_auto14} \\tag{32}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"provided $i$ and $j$ are integers.\n",
"\n",
"<!-- Sympy can do this! -->\n",
"<!-- k, m, n = symbols('k m n', integer=True) -->\n",
"<!-- >>> integrate(sin(k*x)*sin(m*x), (x, 0, 2*pi)) -->\n",
"<!-- 0 -->\n",
"<!-- >>>integrate(sin(k*x)*sin(k*x), (x, 0, 2*pi)) -->\n",
"<!-- pi -->\n",
"\n",
"### The residual\n",
"\n",
"We can readily calculate the following explicit expression for the\n",
"residual:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"R(x;c_0, \\ldots, c_N) = u''(x) + f(x),\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"= \\frac{d^2}{dx^2}\\left(\\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(x)\\right)\n",
"+ f(x),\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:ex:sines:res\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"= \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j''(x) + f(x){\\thinspace .}\n",
"\\label{fem:deq:1D:ex:sines:res} \\tag{33}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The least squares method\n",
"\n",
"The equations ([21](#fem:deq:1D:LS:eq1))\n",
"in the least squares method require an expression for\n",
"$\\partial R/\\partial c_i$. We have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto15\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial R}{\\partial c_i} =\n",
"\\frac{\\partial}{\\partial c_i}\n",
"\\left(\\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j''(x) + f(x)\\right)\n",
"= \\sum_{j\\in{\\mathcal{I}_s}} \\frac{\\partial c_j}{\\partial c_i}{\\psi}_j''(x)\n",
"= {\\psi}_i''(x){\\thinspace .} \\label{_auto15} \\tag{34}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The governing equations for the unknown parameters $\\left\\{ {c}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$ are then"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto16\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(\\sum_j c_j {\\psi}_j'' + f,{\\psi}_i'')=0,\\quad i\\in{\\mathcal{I}_s},\n",
"\\label{_auto16} \\tag{35}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which can be rearranged as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto17\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\sum_{j\\in{\\mathcal{I}_s}}({\\psi}_i'',{\\psi}_j'')c_j = -(f,{\\psi}_i''),\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"\\label{_auto17} \\tag{36}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is nothing but a linear system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sum_{j\\in{\\mathcal{I}_s}}A_{i,j}c_j = b_i,\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The entries in the coefficient matrix are given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"A_{i,j} &= ({\\psi}_i'',{\\psi}_j'')\\nonumber\\\\ \n",
"& = \\pi^4(i+1)^2(j+1)^2L^{-4}\\int_0^L {\\sin\\left((i+1)\\pi\\frac{x}{L}\\right)}{\\sin\\left((j+1)\\pi\\frac{x}{L}\\right)}\\, {\\, \\mathrm{d}x}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The orthogonality of the sine functions simplify the coefficient matrix:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto18\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"A_{i,j} = \\left\\lbrace \\begin{array}{ll}\n",
"{1\\over2}L^{-3}\\pi^4(i+1)^4 & i=j \\\\ \n",
"0, & i\\neq j\n",
"\\end{array}\\right.\n",
"\\label{_auto18} \\tag{37}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The right-hand side reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto19\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"b_i = -(f,{\\psi}_i'') = (i+1)^2\\pi^2L^{-2}\\int_0^Lf(x){\\sin\\left((i+1)\\pi\\frac{x}{L}\\right)}\\, {\\, \\mathrm{d}x}\n",
"\\label{_auto19} \\tag{38}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the coefficient matrix is diagonal we can easily solve for"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:ex:sines:solution\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"c_i = \\frac{2L}{\\pi^2(i+1)^2}\\int_0^Lf(x){\\sin\\left((i+1)\\pi\\frac{x}{L}\\right)}\\, {\\, \\mathrm{d}x}{\\thinspace .}\n",
"\\label{fem:deq:1D:ex:sines:solution} \\tag{39}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the special choice of $f(x)=2$, the coefficients\n",
"can be calculated in `sympy` by"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Piecewise((4*L**2*((-1)**i + 1)/(pi**3*(i + 1)**3), Ne(pi*(i + 1)/L, 0)), (0, True))\n"
]
}
],
"source": [
"import sympy as sym\n",
"\n",
"i, j = sym.symbols('i j', integer=True)\n",
"x, L = sym.symbols('x L')\n",
"f = 2\n",
"a = 2*L/(sym.pi**2*(i+1)**2)\n",
"c_i = a*sym.integrate(f*sym.sin((i+1)*sym.pi*x/L), (x, 0, L))\n",
"c_i = sym.simplify(c_i)\n",
"print(c_i)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The answer becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"c_i = 4 \\frac{L^{2} \\left(\\left(-1\\right)^{i} + 1\\right)}{\\pi^{3}\n",
"\\left(i^{3} + 3 i^{2} + 3 i + 1\\right)}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, $1+(-1)^i=0$ for $i$ odd, so only the coefficients with even index\n",
"are nonzero. Introducing $i=2k$ for $k=0,\\ldots,N/2$ to count the\n",
"relevant indices (for $N$ odd, $k$ goes to $(N-1)/2$), we get the solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto20\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(x) = \\sum_{k=0}^{N/2} \\frac{8L^2}{\\pi^3(2k+1)^3}{\\sin\\left((2k+1)\\pi\\frac{x}{L}\\right)}{\\thinspace .} \\label{_auto20} \\tag{40}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The coefficients decay very fast: $c_2 = c_0/27$, $c_4=c_0/125$.\n",
"The solution will therefore be dominated by the first term,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x) \\approx \\frac{8L^2}{\\pi^3}\\sin\\left(\\pi\\frac{x}{L}\\right){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### The Galerkin method\n",
"\n",
"The Galerkin principle ([22](#fem:deq:1D:Galerkin0))\n",
"applied to ([30](#fem:deq:1D:model1b)) consists of inserting\n",
"our special residual ([33](#fem:deq:1D:ex:sines:res)) in\n",
"([22](#fem:deq:1D:Galerkin0))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(u''+f,v)=0,\\quad \\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto21\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(u'',v) = -(f,v),\\quad\\forall v\\in V{\\thinspace .} \\label{_auto21} \\tag{41}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the variational formulation, based on the Galerkin principle,\n",
"of our differential equation.\n",
"The $\\forall v\\in V$ requirement is equivalent to\n",
"demanding the equation $(u'',v) = -(f,v)$ to be fulfilled for all\n",
"basis functions $v={\\psi}_i$, $i\\in{\\mathcal{I}_s}$, see\n",
"([22](#fem:deq:1D:Galerkin0)) and ([23](#fem:deq:1D:Galerkin)).\n",
"We therefore have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto22\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(\\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j'', {\\psi}_i)=-(f,{\\psi}_i),\\quad i\\in{\\mathcal{I}_s}{\\thinspace .} \\label{_auto22} \\tag{42}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This equation can be rearranged to a form that explicitly shows\n",
"that we get a linear system for the unknowns $\\left\\{ {c}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto23\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\sum_{j\\in{\\mathcal{I}_s}} ({\\psi}_i,{\\psi}_j'')c_j = (f, {\\psi}_i),\\quad i\\in{\\mathcal{I}_s}{\\thinspace .} \\label{_auto23} \\tag{43}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For the particular choice of the basis functions ([31](#fem:deq:1D:ex:sines:psi))\n",
"we get in fact the same linear system\n",
"as in the least squares method\n",
"because ${\\psi}''= -(i+1)^2\\pi^2L^{-2}{\\psi}$.\n",
"Consequently, the solution $u(x)$ becomes identical to the one produced\n",
"by the least squares method.\n",
"\n",
"### The collocation method\n",
"\n",
"For the collocation method ([26](#fem:deq:1D:collocation)) we need to\n",
"decide upon a set of $N+1$ collocation points in $\\Omega$. A simple\n",
"choice is to use uniformly spaced points: $x_{i}=i\\Delta x$, where\n",
"$\\Delta x = L/N$ in our case ($N\\geq 1$). However, these points\n",
"lead to at least two rows in the matrix consisting of zeros\n",
"(since ${\\psi}_i(x_{0})=0$ and ${\\psi}_i(x_{N})=0$), thereby making the matrix\n",
"singular and non-invertible. This forces us to choose some other\n",
"collocation points, e.g., random points or points uniformly distributed\n",
"in the interior of $\\Omega$.\n",
"Demanding the residual to vanish\n",
"at these points leads, in our model problem ([30](#fem:deq:1D:model1b)), to\n",
"the equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto24\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"-\\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j''(x_{i}) = f(x_{i}),\\quad i\\in{\\mathcal{I}_s},\n",
"\\label{_auto24} \\tag{44}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is seen to be a linear system with entries"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"A_{i,j}=-{\\psi}_j''(x_{i})=\n",
"(j+1)^2\\pi^2L^{-2}\\sin\\left((j+1)\\pi \\frac{x_i}{L}\\right),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"in the coefficient matrix and entries\n",
"$b_i=2$ for the right-hand side (when $f(x)=2$).\n",
"\n",
"The special case of $N=0$\n",
"can sometimes be of interest. A natural choice is then the midpoint\n",
"$x_{0}=L/2$ of the domain, resulting in\n",
"$A_{0,0} = -{\\psi}_0''(x_{0}) = \\pi^2L^{-2}$, $f(x_0)=2$,\n",
"and hence $c_0=2L^2/\\pi^2$.\n",
"\n",
"\n",
"### Comparison\n",
"\n",
"In the present model problem, with $f(x)=2$, the exact solution is\n",
"$u(x)=x(L-x)$, while for $N=0$ the Galerkin and least squares method\n",
"result in $u(x)=8L^2\\pi^{-3}\\sin (\\pi x/L)$ and the\n",
"collocation method leads to $u(x)=2L^2\\pi^{-2}\\sin (\\pi x/L)$.\n",
"We can quickly use `sympy` to verify that the maximum error\n",
"occurs at the midpoint $x=L/2$ and find what the errors are.\n",
"First we set up the error expressions:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"If the derivative of the errors vanish at $x=L/2$, the errors reach\n",
"their maximum values here (the errors vanish at the boundary points)."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 0$"
],
"text/plain": [
"0"
]
},
"execution_count": 24,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"x, L = sym.symbols('x L')\n",
"e_Galerkin = x*(L-x) - 8*L**2*sym.pi**(-3)*sym.sin(sym.pi*x/L)\n",
"dedx_Galerkin = sym.diff(e_Galerkin, x)\n",
"dedx_Galerkin.subs(x, L/2)"
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 0$"
],
"text/plain": [
"0"
]
},
"execution_count": 25,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"e_colloc = x*(L-x) - 2*L**2*sym.pi**(-2)*sym.sin(sym.pi*x/L)\n",
"dedx_colloc = sym.diff(e_colloc, x)\n",
"dedx_colloc.subs(x, L/2)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we can compute the maximum error at $x=L/2$ and evaluate\n",
"the expressions numerically with three decimals:"
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle - 0.00812 L^{2}$"
],
"text/plain": [
"-0.00812*L**2"
]
},
"execution_count": 26,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sym.simplify(e_Galerkin.subs(x, L/2).evalf(n=3))"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {},
"outputs": [
{
"data": {
"text/latex": [
"$\\displaystyle 0.0473 L^{2}$"
],
"text/plain": [
"0.0473*L**2"
]
},
"execution_count": 27,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"sym.simplify(e_colloc.subs(x, L/2).evalf(n=3))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The error in the collocation method is about 6 times larger than\n",
"the error in the Galerkin or least squares method.\n",
"\n",
"\n",
"## Integration by parts\n",
"<div id=\"fem:deq:1D:varform\"></div>\n",
"\n",
"\n",
"A problem arises if we want to apply popular finite element functions\n",
"to solve our model problem ([30](#fem:deq:1D:model1b))\n",
"by the standard least squares, Galerkin, or collocation methods: the piecewise\n",
"polynomials ${\\psi}_i(x)$ have discontinuous derivatives at the\n",
"cell boundaries which makes it problematic to compute\n",
"the second-order derivative. This fact actually makes the least squares and\n",
"collocation methods less suitable for finite element approximation of\n",
"the unknown function. (By rewriting the equation $-u''=f$ as a\n",
"system of two first-order equations, $u'=v$ and $-v'=f$, the\n",
"least squares method can be applied. Also, differentiating discontinuous\n",
"functions can actually be handled by distribution theory in\n",
"mathematics.) The Galerkin method and the method of\n",
"weighted residuals can, however, be applied together with finite\n",
"element basis functions if we use *integration by parts*\n",
"as a means for transforming a second-order derivative to a first-order\n",
"one.\n",
"\n",
"Consider the model problem ([30](#fem:deq:1D:model1b)) and its\n",
"Galerkin formulation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-(u'',v) = (f,v)\\quad\\forall v\\in V{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using integration by parts in the Galerkin method,\n",
"we can \"move\" a derivative of $u$ onto $v$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_0^L u''(x)v(x) {\\, \\mathrm{d}x} = - \\int_0^Lu'(x)v'(x){\\, \\mathrm{d}x}\n",
"+ [vu']_0^L\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:intbyparts\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"= - \\int_0^Lu'(x)v'(x) {\\, \\mathrm{d}x}\n",
"+ u'(L)v(L) - u'(0)v(0){\\thinspace .}\n",
"\\label{fem:deq:1D:intbyparts} \\tag{45}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Usually, one integrates the problem at the stage where the $u$ and $v$\n",
"functions enter the formulation.\n",
"Alternatively, but less common, we can integrate by parts in the expressions for\n",
"the matrix entries:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_0^L{\\psi}_i(x){\\psi}_j''(x) {\\, \\mathrm{d}x} =\n",
"- \\int_0^L{\\psi}_i'(x){\\psi}_j'(x) dx\n",
"+ [{\\psi}_i{\\psi}_j']_0^L\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:intbyparts0\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"= - \\int_0^L{\\psi}_i'(x){\\psi}_j'(x) {\\, \\mathrm{d}x}\n",
"+ {\\psi}_i(L){\\psi}_j'(L) - {\\psi}_i(0){\\psi}_j'(0){\\thinspace .}\n",
"\\label{fem:deq:1D:intbyparts0} \\tag{46}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Integration by parts serves to reduce the order of the derivatives and\n",
"to make the coefficient matrix symmetric since\n",
"$({\\psi}_i',{\\psi}_j') = ({\\psi}_j',{\\psi}_i')$.\n",
"The symmetry property depends\n",
"on the type of terms that enter the differential equation.\n",
"As will be seen later,\n",
"integration by parts also provides a method for implementing\n",
"boundary conditions involving $u'$.\n",
"\n",
"With the choice ([31](#fem:deq:1D:ex:sines:psi)) of basis functions we see\n",
"that the \"boundary terms\"\n",
"${\\psi}_i(L){\\psi}_j'(L)$ and ${\\psi}_i(0){\\psi}_j'(0)$\n",
"vanish since ${\\psi}_i(0)={\\psi}_i(L)=0$.\n",
"<!-- A boundary term associated with -->\n",
"<!-- a location at the boundary where we have Dirichlet conditions will always -->\n",
"<!-- vanish because ${\\psi}_i=0$ at such locations. -->\n",
"We therefore end up with the following alternative Galerkin formulation:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-(u'',v) = (u', v') = (f,v)\\quad \\forall v\\in V{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Weak form\n",
"\n",
"Since the variational formulation after integration by parts make\n",
"weaker demands on the differentiability of $u$ and the basis\n",
"functions ${\\psi}_i$,\n",
"the resulting integral formulation is referred to as a *weak form* of\n",
"the differential equation problem. The original variational formulation\n",
"with second-order derivatives, or the differential equation problem\n",
"with second-order derivative, is then the *strong form*, with\n",
"stronger requirements on the differentiability of the functions.\n",
"\n",
"For differential equations with second-order derivatives, expressed as\n",
"variational formulations and solved by finite element methods, we will\n",
"always perform integration by parts to arrive at expressions involving\n",
"only first-order derivatives.\n",
"\n",
"\n",
"## Boundary function\n",
"<div id=\"fem:deq:1D:essBC:Bfunc\"></div>\n",
"\n",
"So far we have assumed zero Dirichlet boundary conditions, typically\n",
"$u(0)=u(L)=0$, and we have demanded that ${\\psi}_i(0)={\\psi}_i(L)=0$\n",
"for $i\\in{\\mathcal{I}_s}$. What about a boundary condition like $u(L)=D\\neq0$?\n",
"This condition immediately faces a problem:\n",
"$u = \\sum_j c_j{\\varphi}_j(L) = 0$ since all ${\\varphi}_i(L)=0$.\n",
"\n",
"We remark that we faced exactly the same problem where\n",
"we considered Fourier series approximations of functions that where non-zero at the boundaries.\n",
"We will use the same trick as we did earlier to get around this problem.\n",
"\n",
"A boundary condition of the form $u(L)=D$ can be implemented by\n",
"demanding that all ${\\psi}_i(L)=0$, but adding a\n",
"*boundary function* $B(x)$ with the right boundary value, $B(L)=D$, to\n",
"the expansion for $u$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x) = B(x) + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(x)\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This $u$ gets the right value at $x=L$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(L) = B(L) + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(L) = B(L) = D{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The idea is that for any boundary where $u$ is known we demand ${\\psi}_i$ to\n",
"vanish and construct a function $B(x)$ to attain the boundary value of $u$.\n",
"There are no restrictions on how $B(x)$ varies with $x$ in the interior of the\n",
"domain, so this variation needs to be constructed in some way. Exactly how\n",
"we decide the variation to be, is not important.\n",
"\n",
"For example, with $u(0)=0$ and\n",
"$u(L)=D$, we can choose $B(x)=x D/L$, since this form ensures that\n",
"$B(x)$ fulfills the boundary conditions: $B(0)=0$ and $B(L)=D$.\n",
"The unknown function is then sought on the form"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:essBC:Bfunc:u1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(x) = \\frac{x}{L}D + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(x),\n",
"\\label{fem:deq:1D:essBC:Bfunc:u1} \\tag{47}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with ${\\psi}_i(0)={\\psi}_i(L)=0$.\n",
"\n",
"The particular shape of the $B(x)$ function is not important\n",
"as long as its boundary\n",
"values are correct. For example, $B(x)=D(x/L)^p$ for any power $p$\n",
"will work fine in the above example. Another choice could be\n",
"$B(x)=D\\sin (\\pi x/(2L))$.\n",
"\n",
"As a more general example, consider a domain $\\Omega = [a,b]$\n",
"where the boundary conditions are $u(a)=U_a$ and $u(b)=U_b$. A class\n",
"of possible $B(x)$ functions is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:essBC:Bfunc:gen\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} B(x)=U_a + \\frac{U_b-U_a}{(b-a)^p}(x-a)^p,\\quad p>0\n",
"{\\thinspace .}\n",
"\\label{fem:deq:1D:essBC:Bfunc:gen} \\tag{48}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Real applications will most likely use the simplest version, $p=1$,\n",
"but here such a $p$ parameter was included to demonstrate that there\n",
"are many choices of $B(x)$ in a problem. Fortunately, there is a general, unique\n",
"technique for constructing $B(x)$ when we use finite element basis functions for\n",
"$V$.\n",
"\n",
"\n",
"**How to deal with nonzero Dirichlet conditions.**\n",
"\n",
"The general procedure of incorporating Dirichlet boundary\n",
"conditions goes as follows.\n",
"Let $\\partial\\Omega_E$ be the part(s) of the boundary\n",
"$\\partial\\Omega$ of the domain $\\Omega$ where $u$ is specified.\n",
"Set ${\\psi}_i=0$ at the points in $\\partial\\Omega_E$ and seek $u$\n",
"as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:1D:essBC:Bfunc:u2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(x) = B(x) + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(x),\n",
"\\label{fem:deq:1D:essBC:Bfunc:u2} \\tag{49}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $B(x)$ equals the boundary conditions on $u$ at $\\partial\\Omega_E$.\n",
"\n",
"\n",
"\n",
"**Remark.**\n",
"With the $B(x)$ term, $u$ does not in general lie in $V=\\hbox{span}\\,\n",
"\\{{\\psi}_0,\\ldots,{\\psi}_N\\}$ anymore. Moreover, when a prescribed value\n",
"of $u$ at the boundary, say $u(a)=U_a$ is different from zero, it does\n",
"not make sense to say that $u$ lies in a vector space, because\n",
"this space does not obey the requirements of addition and scalar multiplication.\n",
"For example,\n",
"$2u$ does not lie in the space since its boundary value is $2U_a$,\n",
"which is incorrect. It only makes sense to split $u$ in two parts,\n",
"as done above, and have the unknown part $\\sum_j c_j {\\psi}_j$ in a\n",
"proper function space.\n",
"<!-- Sometimes it is said that $u$ is in the *affine space* $B+V$. -->\n",
"\n",
"# Computing with global polynomials\n",
"\n",
"The next example uses global polynomials and shows\n",
"that if our solution, modulo boundary conditions, lies in the space spanned\n",
"by these polynomials, then the Galerkin method recovers the exact solution.\n",
"\n",
"\n",
"\n",
"## Computing with Dirichlet and Neumann conditions\n",
"<div id=\"fem:deq:1D:varform:ex:DN:case\"></div>\n",
"\n",
"<!-- ex_varform1D.py: case2 -->\n",
"\n",
"Let us perform the necessary calculations to solve"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-u''(x)=2,\\quad x\\in \\Omega=[0,1],\\quad u'(0)=C,\\ u(1)=D,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"using a global polynomial basis ${\\psi}_i\\sim x^i$.\n",
"The requirements on ${\\psi}_i$ is that ${\\psi}_i(1)=0$, because $u$ is\n",
"specified at $x=1$, so a proper set of polynomial basis functions can be"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"{\\psi}_i(x)=(1-x)^{i+1}, \\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A suitable $B(x)$ function\n",
"to handle the boundary condition $u(1)=D$ is $B(x)=Dx$.\n",
"The variational formulation becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(u',v') = (2,v) - Cv(0)\\quad\\forall v\\in V{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From inserting $u=B + \\sum_{j}c_j{\\psi}_j$ and choosing $v={\\psi}_i$ we get"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sum_{j\\in{\\mathcal{I}_s}} ({\\psi}_j',{\\psi}_i')c_j = (2,{\\psi}_i)\n",
"- (B',{\\psi}_i') - C{\\psi}_i(0),\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The entries in the linear system are then"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"A_{i,j} &= ({\\psi}_j',{\\psi}_i') = \\int_{0}^1 {\\psi}_i'(x){\\psi}_j'(x){\\, \\mathrm{d}x}\n",
"= \\int_0^1 (i+1)(j+1)(1-x)^{i+j}{\\, \\mathrm{d}x}\\\\ \n",
"&= \\frac{(i+1)(j+1)}{i + j + 1},\\\\ \n",
"b_i &= (2,{\\psi}_i) - (D,{\\psi}_i') -C{\\psi}_i(0)\\\\ \n",
"&= \\int_0^1\\left( 2{\\psi}_i(x) - D{\\psi}_i'(x)\\right){\\, \\mathrm{d}x} -C{\\psi}_i(0)\\\\ \n",
"&= \\int_0^1 \\left( 2(1-x)^{i+1} + D(i+1)(1-x)^i\\right){\\, \\mathrm{d}x} -C\\\\ \n",
"&= \\frac{(D-C)(i+2) + 2}{i+2} = D - C + \\frac{2}{i+2}\n",
"{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Relevant `sympy` commands to help calculate these expressions are"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"('A_ij:', (i + 1)*(j + 1)/(i + j + 1))\n",
"('b_i:', ((-C + D)*(i + 2) + 2)/(i + 2))\n"
]
}
],
"source": [
"from sympy import *\n",
"x, C, D = symbols('x C D')\n",
"i, j = symbols('i j', integer=True, positive=True)\n",
"psi_i = (1-x)**(i+1)\n",
"psi_j = psi_i.subs(i, j)\n",
"integrand = diff(psi_i, x)*diff(psi_j, x)\n",
"integrand = simplify(integrand)\n",
"A_ij = integrate(integrand, (x, 0, 1))\n",
"A_ij = simplify(A_ij)\n",
"print(('A_ij:', A_ij))\n",
"f = 2\n",
"b_i = integrate(f*psi_i, (x, 0, 1)) - \\\n",
" integrate(diff(D*x, x)*diff(psi_i, x), (x, 0, 1)) - \\\n",
" C*psi_i.subs(x, 0)\n",
"b_i = simplify(b_i)\n",
"print(('b_i:', b_i))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The output becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" A_ij: (i + 1)*(j + 1)/(i + j + 1)\n",
" b_i: ((-C + D)*(i + 2) + 2)/(i + 2)\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now choose some $N$ and form the linear system, say for $N=1$:"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"('fresh b:', Matrix([\n",
"[0, 0],\n",
"[0, 0]]))\n"
]
}
],
"source": [
"N = 1\n",
"A = zeros(N+1, N+1)\n",
"b = zeros(N+1)\n",
"print(('fresh b:', b))\n",
"for r in range(N+1):\n",
" for s in range(N+1):\n",
" A[r,s] = A_ij.subs(i, r).subs(j, s)\n",
" b[r,0] = b_i.subs(i, r)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The system becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left(\\begin{array}{cc}\n",
"1 & 1\\\\ \n",
"1 & 4/3\n",
"\\end{array}\\right)\n",
"\\left(\\begin{array}{c}\n",
"c_0\\\\ \n",
"c_1\n",
"\\end{array}\\right)\n",
"=\n",
"\\left(\\begin{array}{c}\n",
"1-C+D\\\\ \n",
"2/3 -C + D\n",
"\\end{array}\\right)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution (`c = A.LUsolve(b)`)\n",
"becomes $c_0=2 -C+D$ and $c_1=-1$, resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto25\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(x) = 1 -x^2 + D + C(x-1),\n",
"\\label{_auto25} \\tag{50}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can form this $u$ in `sympy` and check that the differential equation\n",
"and the boundary conditions are satisfied:"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"('u:', C*x - C + D - x**2 + 1)\n",
"(\"u'':\", -2)\n",
"('BC x=0:', C)\n",
"('BC x=1:', D)\n"
]
}
],
"source": [
"c = A.LUsolve(b)\n",
"u = sum(c[r,0]*psi_i.subs(i, r) for r in range(N+1)) + D*x\n",
"print(('u:', simplify(u)))\n",
"print((\"u'':\", simplify(diff(u, x, x))))\n",
"print(('BC x=0:', simplify(diff(u, x).subs(x, 0))))\n",
"print(('BC x=1:', simplify(u.subs(x, 1))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The output becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
" u: C*x - C + D - x**2 + 1\n",
" u'': -2\n",
" BC x=0: C\n",
" BC x=1: D\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The complete `sympy` code is found in [`u_xx_2_CD.py`](src/u_xx_2_CD.py).\n",
"\n",
"The exact solution is found by integrating twice and applying the\n",
"boundary conditions, either by hand or using `sympy` as shown in\n",
"the section [Simple model problems and their solutions](#fem:deq:1D:models:simple). It appears that the numerical\n",
"solution coincides with the exact one."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Abstract notation for variational formulations\n",
"<div id=\"fem:deq:1D:varform:abstract\"></div>\n",
"\n",
"We have seen that variational formulations end up with a formula involving\n",
"$u$ and $v$, such as $(u',v')$ and a formula involving $v$ and known\n",
"functions, such as $(f,v)$. A widely used notation is to introduce an abstract\n",
"variational statement written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a(u,v)=L(v)\\quad\\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $a(u,v)$ is a so-called *bilinear form* involving all the terms\n",
"that contain both the test and trial\n",
"function, while $L(v)$ is a *linear form* containing all the terms without\n",
"the trial function. For example, the statement"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{\\Omega} u' v' {\\, \\mathrm{d}x} =\n",
"\\int_{\\Omega} fv{\\, \\mathrm{d}x}\\quad\\hbox{or}\\quad (u',v') = (f,v)\n",
"\\quad\\forall v\\in V\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"can be written in abstract form: *find $u$ such that*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a(u,v) = L(v)\\quad \\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where we have the definitions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a(u,v) = (u',v'),\\quad L(v) = (f,v){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The term *linear* means that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"L(\\alpha_1 v_1 + \\alpha_2 v_2) =\\alpha_1 L(v_1) + \\alpha_2 L(v_2)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for two test functions $v_1$ and $v_2$, and\n",
"scalar parameters $\\alpha_1$ and $\\alpha_2$. Similarly, the term *bilinear*\n",
"means that $a(u,v)$ is linear in both its arguments:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"a(\\alpha_1 u_1 + \\alpha_2 u_2, v) &= \\alpha_1 a(u_1,v) + \\alpha_2 a(u_2, v),\n",
"\\\\ \n",
"a(u, \\alpha_1 v_1 + \\alpha_2 v_2) &= \\alpha_1 a(u,v_1) + \\alpha_2 a(u, v_2)\n",
"{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In nonlinear problems these linearity properties do not hold in general\n",
"and the abstract notation is then"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F(u;v)=0\\quad\\forall v\\in V{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The matrix system associated with $a(u,v)=L(v)$ can also be written in\n",
"an abstract form by inserting $v={\\psi}_i$ and $u=\\sum_j c_j{\\psi}_j$\n",
"in $a(u,v)=L(v)$. Using the linear properties, we get"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sum_{j\\in{\\mathcal{I}_s}} a({\\psi}_j,{\\psi}_i) c_j = L({\\psi}_i),\\quad i\\in{\\mathcal{I}_s},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is a linear system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sum_{j\\in{\\mathcal{I}_s}}A_{i,j}c_j = b_i,\\quad i\\in{\\mathcal{I}_s},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"A_{i,j} =a({\\psi}_j,{\\psi}_i), \\quad b_i = L({\\psi}_i){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In many problems, $a(u,v)$ is symmetric such that\n",
"$a({\\psi}_j,{\\psi}_i) = a({\\psi}_i,{\\psi}_j)$. In those cases the\n",
"coefficient matrix becomes symmetric, $A_{i,j}=A_{j,i}$, a property\n",
"that can simplify solution algorithms for linear systems\n",
"and make them more stable. The property also reduces memory\n",
"requirements and the computational work.\n",
"\n",
"\n",
"The abstract notation $a(u,v)=L(v)$ for linear differential equation problems\n",
"is much used in the literature and\n",
"in description of finite element software (in particular the\n",
"[FEniCS](http://fenicsproject.org) documentation). We shall\n",
"frequently summarize variational forms using this notation."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Examples on variational formulations\n",
"<div id=\"fem:deq:1D:varform:ex\"></div>\n",
"\n",
"The following sections derive variational formulations for some\n",
"prototype differential equations in 1D, and demonstrate how we with\n",
"ease can handle variable coefficients, mixed Dirichlet and Neumann\n",
"boundary conditions, first-order derivatives, and nonlinearities.\n",
"\n",
"## Variable coefficient\n",
"\n",
"Consider the problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto26\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"-\\frac{d}{dx}\\left( {\\alpha}(x)\\frac{du}{dx}\\right) = f(x),\\quad x\\in\\Omega =[0,L],\\ \n",
"u(0)=C,\\ u(L)=D{\\thinspace .}\n",
"\\label{_auto26} \\tag{51}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"There are two new features of this problem compared with\n",
"previous examples: a variable\n",
"coefficient ${\\alpha} (x)$ and nonzero Dirichlet conditions at both boundary points.\n",
"\n",
"Let us first deal with the boundary conditions. We seek"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x) = B(x) + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_i(x){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the Dirichlet conditions demand"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"{\\psi}_i(0)={\\psi}_i(L)=0,\\quad i\\in{\\mathcal{I}_s},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"the function $B(x)$\n",
"must fulfill $B(0)=C$ and $B(L)=D$. The we are guaranteed that $u(0)=C$\n",
"and $u(L)=D$. How $B$ varies in between\n",
"$x=0$ and $x=L$ is not of importance. One possible choice is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"B(x) = C + \\frac{1}{L}(D-C)x,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which follows from ([48](#fem:deq:1D:essBC:Bfunc:gen)) with $p=1$.\n",
"\n",
"We seek $(u-B)\\in V$. As usual,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"V = \\hbox{span}\\{{\\psi}_0,\\ldots,{\\psi}_N\\}{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that any $v\\in V$ has the property $v(0)=v(L)=0$.\n",
"\n",
"The residual arises by inserting our $u$ in the differential equation:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"R = -\\frac{d}{dx}\\left( {\\alpha}\\frac{du}{dx}\\right) -f{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Galerkin's method is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(R, v) = 0,\\quad \\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or written with explicit integrals,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{\\Omega} \\left(-\\frac{d}{dx}\\left( {\\alpha}\\frac{du}{dx}\\right) -f\\right)v {\\, \\mathrm{d}x} = 0,\\quad \\forall v\\in V {\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We proceed with integration by parts to lower the derivative from\n",
"second to first order:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-\\int_{\\Omega} \\frac{d}{dx}\\left( {\\alpha}(x)\\frac{du}{dx}\\right) v {\\, \\mathrm{d}x}\n",
"= \\int_{\\Omega} {\\alpha}(x)\\frac{du}{dx}\\frac{dv}{dx}{\\, \\mathrm{d}x} -\n",
"\\left[{\\alpha}\\frac{du}{dx}v\\right]_0^L\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The boundary term vanishes since $v(0)=v(L)=0$.\n",
"The variational formulation is then"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{\\Omega} {\\alpha}(x)\\frac{du}{dx}\\frac{dv}{dx}{\\, \\mathrm{d}x} = \\int_{\\Omega} f(x)v{\\, \\mathrm{d}x},\\quad\n",
"\\forall v\\in V{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The variational formulation can alternatively be written in a more\n",
"compact form:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"({\\alpha} u',v') = (f,v),\\quad \\forall v\\in V\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The corresponding abstract notation reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a(u,v)=L(v)\\quad\\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a(u,v)= ({\\alpha} u',v'),\\quad L(v)=(f,v) {\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We may insert $u=B + \\sum_jc_j{\\psi}_j$ and $v={\\psi}_i$ to\n",
"derive the linear system:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"({\\alpha} B' + {\\alpha} \\sum_{j\\in{\\mathcal{I}_s}} c_j {\\psi}_j', {\\psi}_i') =\n",
"(f,{\\psi}_i), \\quad i\\in{\\mathcal{I}_s} {\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Isolating everything with the $c_j$ coefficients on the left-hand side\n",
"and all known terms on the right-hand side\n",
"gives"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sum_{j\\in{\\mathcal{I}_s}} ({\\alpha}{\\psi}_j', {\\psi}_i')c_j =\n",
"(f,{\\psi}_i) + (\\alpha (D-C)L^{-1}, {\\psi}_i'), \\quad i\\in{\\mathcal{I}_s}\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is nothing but a linear system $\\sum_j A_{i,j}c_j=b_i$\n",
"with"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"A_{i,j} &= (\\alpha {\\psi}_j', {\\psi}_i') = \\int_{\\Omega} {\\alpha}(x){\\psi}_j'(x),\n",
"{\\psi}_i'(x){\\, \\mathrm{d}x},\\\\ \n",
"b_i &= (f,{\\psi}_i) + (\\alpha (D-C)L^{-1},{\\psi}_i')=\n",
"\\int_{\\Omega} \\left(f(x){\\psi}_i(x) + {\\alpha}(x)\\frac{D-C}{L}{\\psi}_i'(x)\\right) {\\, \\mathrm{d}x}\n",
"{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## First-order derivative in the equation and boundary condition\n",
"\n",
"The next problem to formulate in terms of a variational form reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto27\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"-u''(x) + bu'(x) = f(x),\\quad x\\in\\Omega =[0,L],\\ \n",
"u(0)=C,\\ u'(L)=E{\\thinspace .}\n",
"\\label{_auto27} \\tag{52}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The new features are a first-order derivative $u'$ in the equation\n",
"and the boundary\n",
"condition involving the derivative: $u'(L)=E$.\n",
"Since we have a Dirichlet condition at $x=0$,\n",
"we must force ${\\psi}_i(0)=0$ and use a boundary function\n",
"to take care of the condition $u(0)=C$.\n",
"Because there is no Dirichlet\n",
"condition on $x=L$ we do not make any requirements to ${\\psi}_i(L)$.\n",
"The simplest possible choice of $B(x)$ is $B(x)=C$.\n",
"\n",
"The expansion for $u$ becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u = C + \\sum_{j\\in{\\mathcal{I}_s}} c_j {\\psi}_i(x)\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The variational formulation arises from multiplying the equation by\n",
"a test function $v\\in V$ and integrating over $\\Omega$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(-u'' + bu' - f, v) = 0,\\quad\\forall v\\in V\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We apply integration by parts to the $u''v$ term only. Although we could\n",
"also integrate $u' v$ by parts, this is not common.\n",
"The result becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(u',v') + (bu',v) = (f,v) + [u' v]_0^L, \\quad\\forall v\\in V {\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, $v(0)=0$ so"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[u' v]_0^L = u'(L)v(L) = E v(L),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"because $u'(L)=E$.\n",
"Thus, integration by parts allows us to take care of the Neumann condition\n",
"in the boundary term.\n",
"\n",
"\n",
"\n",
"**Natural and essential boundary conditions.**\n",
"\n",
"A common mistake is to forget a boundary term like $[u'v]_0^L$ in\n",
"the integration by parts. Such a mistake implies that we actually\n",
"impose the condition $u'=0$ unless there is a Dirichlet condition\n",
"(i.e., $v=0$) at that point! This fact has great practical\n",
"consequences, because it is easy to forget the boundary term, and that\n",
"implicitly set a boundary condition!\n",
"\n",
"Since homogeneous Neumann conditions can be incorporated without\n",
"\"doing anything\" (i.e., omitting the boundary term), and\n",
"non-homogeneous Neumann conditions can just be inserted in the\n",
"boundary term, such conditions are known as *natural boundary\n",
"conditions*. Dirichlet conditions require more essential steps in the\n",
"mathematical formulation, such as forcing all ${\\varphi}_i=0$ on the\n",
"boundary and constructing a $B(x)$, and are therefore known as\n",
"*essential boundary conditions*.\n",
"\n",
"\n",
"\n",
"The final variational form reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(u',v') + (bu',v) = (f,v) + E v(L), \\quad\\forall v\\in V {\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In the abstract notation we have"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a(u,v)=L(v)\\quad\\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with the particular formulas"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"a(u,v)=(u',v') + (bu',v),\\quad L(v)= (f,v) + E v(L){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The associated linear system is derived by inserting $u=B+\\sum_jc_j{\\psi}_j$\n",
"and replacing $v$ by ${\\psi}_i$ for $i\\in{\\mathcal{I}_s}$. Some algebra results in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sum_{j\\in{\\mathcal{I}_s}} \\underbrace{(({\\psi}_j',{\\psi}_i') + (b{\\psi}_j',{\\psi}_i))}_{A_{i,j}} c_j = \\underbrace{(f,{\\psi}_i) + E {\\psi}_i(L)}_{b_i}\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Observe that in this problem, the coefficient matrix is not symmetric,\n",
"because of the term"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(b{\\psi}_j',{\\psi}_i)=\\int_{\\Omega} b{\\psi}_j'{\\psi}_i {\\, \\mathrm{d}x}\n",
" \\neq \\int_{\\Omega} b {\\psi}_i' {\\psi}_j {\\, \\mathrm{d}x} = ({\\psi}_i',b{\\psi}_j)\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Too early: -->\n",
"<!-- For finite element basis functions, it is worth noticing that the boundary term -->\n",
"<!-- $E{\\psi}_i(L)$ is nonzero only in the entry $b_N$ since all -->\n",
"<!-- ${\\psi}_i$, $i\\neq N$, are zero at $x=L$, provided the degrees of freedom -->\n",
"<!-- are numbered from left to right in 1D so that $x_{N}=L$. -->\n",
"\n",
"## Nonlinear coefficient\n",
"\n",
"Finally, we show that the techniques used above to derive variational\n",
"forms apply to nonlinear differential equation\n",
"problems as well. Here is a model problem with\n",
"a nonlinear coefficient $\\alpha(u)$ and a nonlinear right-hand side $f(u)$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto28\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"-({\\alpha}(u)u')' = f(u),\\quad x\\in [0,L],\\ u(0)=0,\\ u'(L)=E\n",
"{\\thinspace .}\n",
"\\label{_auto28} \\tag{53}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Our space $V$ has basis $\\left\\{ {{\\psi}}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$, and because of the\n",
"condition $u(0)=0$, we must require ${\\psi}_i(0)=0$, $i\\in{\\mathcal{I}_s}$.\n",
"\n",
"Galerkin's method is about inserting the approximate\n",
"$u$, multiplying the differential equation by $v\\in V$, and integrate,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-\\int_0^L \\frac{d}{dx}\\left({\\alpha}(u)\\frac{du}{dx}\\right)v {\\, \\mathrm{d}x} =\n",
"\\int_0^L f(u)v {\\, \\mathrm{d}x}\\quad\\forall v\\in V\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The integration by parts does not differ from the case where we have\n",
"${\\alpha}(x)$ instead of ${\\alpha}(u)$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_0^L {\\alpha}(u)\\frac{du}{dx}\\frac{dv}{dx}{\\, \\mathrm{d}x} =\n",
"\\int_0^L f(u)v{\\, \\mathrm{d}x} + [{\\alpha}(u)vu']_0^L\\quad\\forall v\\in V\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The term ${\\alpha}(u(0))v(0)u'(0)=0$ since $v(0)$.\n",
"The other term, ${\\alpha}(u(L))v(L)u'(L)$,\n",
"is used to impose the other boundary condition $u'(L)=E$, resulting in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_0^L {\\alpha}(u)\\frac{du}{dx}\\frac{dv}{dx}{\\, \\mathrm{d}x} =\n",
"\\int_0^L f(u)v{\\, \\mathrm{d}x} + {\\alpha}(u(L))v(L)E\\quad\\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or alternatively written more compactly as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"({\\alpha}(u)u', v') = (f(u),v) + {\\alpha}(u(L))v(L)E\\quad\\forall v\\in V\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since the problem is nonlinear, we cannot identify a bilinear\n",
"form $a(u,v)$ and a linear form $L(v)$.\n",
"An abstract formulation is typically *find $u$ such that*"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F(u;v) = 0\\quad\\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"F(u;v) = (a(u)u', v') - (f(u),v) - a(L)v(L)E\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By inserting $u=\\sum_j c_j{\\psi}_j$ and $v={\\psi}_i$ in $F(u;v)$,\n",
"we get a *nonlinear system of\n",
"algebraic equations* for the unknowns $c_i$, $i\\in{\\mathcal{I}_s}$. Such systems must\n",
"be solved by constructing a sequence of linear systems whose solutions\n",
"hopefully converge to the solution of the nonlinear system. Frequently applied\n",
"methods are Picard iteration and Newton's method.\n",
"\n",
"# Implementation of the algorithms\n",
"<div id=\"fem:global:deq:1D:code\"></div>\n",
"\n",
"Our hand calculations can benefit greatly by symbolic computing, as shown\n",
"earlier, so it is natural to extend our approximation programs based on\n",
"`sympy` to the problem domain of variational formulations.\n",
"\n",
"## Extensions of the code for approximation\n",
"<div id=\"fem:deq:1D:code:global\"></div>\n",
"\n",
"The user must prepare a function `integrand_lhs(psi, i, j)` for\n",
"returning the integrand of the integral that contributes to matrix\n",
"entry $(i,j)$ on the left-hand side. The `psi` variable is a Python dictionary holding the\n",
"basis functions and their derivatives in symbolic form. More\n",
"precisely, `psi[q]` is a list of"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\{\\frac{d^q{\\psi}_0}{dx^q},\\ldots,\\frac{d^q{\\psi}_{N_n-1}}{dx^q}\\}\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Similarly, `integrand_rhs(psi, i)` returns the integrand\n",
"for entry number $i$ in the right-hand side vector.\n",
"\n",
"Since we also have contributions to the right-hand side vector (and\n",
"potentially also the matrix) from boundary terms without any integral,\n",
"we introduce two additional functions, `boundary_lhs(psi, i, j)` and\n",
"`boundary_rhs(psi, i)` for returning terms in the variational\n",
"formulation that are not to be integrated over the domain $\\Omega$.\n",
"Examples, to be shown later, will explain in more detail how these\n",
"user-supplied functions may look like.\n",
"\n",
"The linear system can be computed and solved symbolically by\n",
"the following function:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"import sympy as sym\n",
"\n",
"def solver(integrand_lhs, integrand_rhs, psi, Omega,\n",
" boundary_lhs=None, boundary_rhs=None):\n",
" N = len(psi[0]) - 1\n",
" A = sym.zeros(N+1, N+1)\n",
" b = sym.zeros(N+1, 1)\n",
" x = sym.Symbol('x')\n",
" for i in range(N+1):\n",
" for j in range(i, N+1):\n",
" integrand = integrand_lhs(psi, i, j)\n",
" I = sym.integrate(integrand, (x, Omega[0], Omega[1]))\n",
" if boundary_lhs is not None:\n",
" I += boundary_lhs(psi, i, j)\n",
" A[i,j] = A[j,i] = I # assume symmetry\n",
" integrand = integrand_rhs(psi, i)\n",
" I = sym.integrate(integrand, (x, Omega[0], Omega[1]))\n",
" if boundary_rhs is not None:\n",
" I += boundary_rhs(psi, i)\n",
" b[i,0] = I\n",
" c = A.LUsolve(b)\n",
" u = sum(c[i,0]*psi[0][i] for i in range(len(psi[0])))\n",
" return u, c"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Fallback to numerical methods\n",
"\n",
"Not surprisingly, symbolic solution of differential\n",
"equations, discretized by a Galerkin or least squares method\n",
"with global basis functions,\n",
"is of limited interest beyond the simplest problems, because\n",
"symbolic integration might be very time consuming or impossible, not\n",
"only in `sympy` but also in\n",
"[WolframAlpha](http://wolframalpha.com)\n",
"(which applies the perhaps most powerful symbolic integration\n",
"software available today: Mathematica). Numerical integration\n",
"as an option is therefore desirable.\n",
"\n",
"The extended `solver` function below tries to combine symbolic and\n",
"numerical integration. The latter can be enforced by the user, or it\n",
"can be invoked after a non-successful symbolic integration (being\n",
"detected by an `Integral` object as the result of the integration\n",
"in `sympy`).\n",
"<!-- see the section [fem:approx:global:Lagrange](#fem:approx:global:Lagrange)). -->\n",
"Note that for a\n",
"numerical integration, symbolic expressions must be converted to\n",
"Python functions (using `lambdify`), and the expressions cannot contain\n",
"other symbols than `x`. The real `solver` routine in the\n",
"[`varform1D.py`](${fem_src}/varform1D.py)\n",
"file has error checking and meaningful error messages in such cases.\n",
"The `solver` code below is a condensed version of the real one, with\n",
"the purpose of showing how to automate the Galerkin or least squares\n",
"method for solving differential equations in 1D with global basis functions:"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [],
"source": [
"def solver(integrand_lhs, integrand_rhs, psi, Omega,\n",
" boundary_lhs=None, boundary_rhs=None, symbolic=True):\n",
" N = len(psi[0]) - 1\n",
" A = sym.zeros(N+1, N+1)\n",
" b = sym.zeros(N+1, 1)\n",
" x = sym.Symbol('x')\n",
" for i in range(N+1):\n",
" for j in range(i, N+1):\n",
" integrand = integrand_lhs(psi, i, j)\n",
" if symbolic:\n",
" I = sym.integrate(integrand, (x, Omega[0], Omega[1]))\n",
" if isinstance(I, sym.Integral):\n",
" symbolic = False # force num.int. hereafter\n",
" if not symbolic:\n",
" integrand_ = sym.lambdify([x], integrand, 'mpmath')\n",
" I = mpmath.quad(integrand_, [Omega[0], Omega[1]])\n",
" if boundary_lhs is not None:\n",
" I += boundary_lhs(psi, i, j)\n",
" A[i,j] = A[j,i] = I\n",
" integrand = integrand_rhs(psi, i)\n",
" if symbolic:\n",
" I = sym.integrate(integrand, (x, Omega[0], Omega[1]))\n",
" if isinstance(I, sym.Integral):\n",
" symbolic = False\n",
" if not symbolic:\n",
" integrand_ = sym.lambdify([x], integrand, 'mpmath')\n",
" I = mpmath.quad(integrand_, [Omega[0], Omega[1]])\n",
" if boundary_rhs is not None:\n",
" I += boundary_rhs(psi, i)\n",
" b[i,0] = I\n",
" c = A.LUsolve(b)\n",
" u = sum(c[i,0]*psi[0][i] for i in range(len(psi[0])))\n",
" return u, c"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example with constant right-hand side\n",
"\n",
"To demonstrate the code above, we address"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-u''(x)=b,\\quad x\\in\\Omega=[0,1],\\quad u(0)=1,\\ u(1)=0,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with $b$ as a (symbolic) constant. A possible basis for the space $V$\n",
"is ${\\psi}_i(x) = x^{i+1}(1-x)$, $i\\in{\\mathcal{I}_s}$. Note that\n",
"${\\psi}_i(0)={\\psi}_i(1)=0$ as required by the Dirichlet conditions.\n",
"We need a $B(x)$ function to take care of the known boundary\n",
"values of $u$. Any function $B(x)=1-x^p$, $p\\in\\mathbb{R}$, is a candidate,\n",
"and one arbitrary choice from this family\n",
"is $B(x)=1-x^3$. The unknown function is then written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x) = B(x) + \\sum_{j\\in{\\mathcal{I}_s}} c_j{\\psi}_j(x){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let us use the Galerkin method to derive the variational formulation.\n",
"Multiplying the differential\n",
"equation by $v$ and integrating by parts yield"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_0^1 u'v' {\\, \\mathrm{d}x} = \\int_0^1 fv {\\, \\mathrm{d}x}\\quad\\forall v\\in V,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and with $u=B + \\sum_jc_j{\\psi}_j$ we get the linear system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto29\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\sum_{j\\in{\\mathcal{I}_s}}\\left(\\int_0^1{\\psi}_i'{\\psi}_j' {\\, \\mathrm{d}x}\\right)c_j =\n",
"\\int_0^1(f{\\psi}_i-B'{\\psi}_i') {\\, \\mathrm{d}x},\n",
"\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"\\label{_auto29} \\tag{54}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The application can be coded as follows with `sympy`:"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"...evaluating matrix... (0,0): (1 - 2*x)**2\n",
"(0,1): (1 - 2*x)*(-x**2 + 2*x*(1 - x))\n",
"(0,2): (1 - 2*x)*(-x**3 + 3*x**2*(1 - x))\n",
"(0,3): (1 - 2*x)*(-x**4 + 4*x**3*(1 - x))\n",
"rhs: b*x*(1 - x) + 3*x**2*(1 - 2*x)\n",
"(1,1): (-x**2 + 2*x*(1 - x))**2\n",
"(1,2): (-x**2 + 2*x*(1 - x))*(-x**3 + 3*x**2*(1 - x))\n",
"(1,3): (-x**2 + 2*x*(1 - x))*(-x**4 + 4*x**3*(1 - x))\n",
"rhs: b*x**2*(1 - x) + 3*x**2*(-x**2 + 2*x*(1 - x))\n",
"(2,2): (-x**3 + 3*x**2*(1 - x))**2\n",
"(2,3): (-x**3 + 3*x**2*(1 - x))*(-x**4 + 4*x**3*(1 - x))\n",
"rhs: b*x**3*(1 - x) + 3*x**2*(-x**3 + 3*x**2*(1 - x))\n",
"(3,3): (-x**4 + 4*x**3*(1 - x))**2\n",
"rhs: b*x**4*(1 - x) + 3*x**2*(-x**4 + 4*x**3*(1 - x))\n",
"\n",
"A:\n",
" Matrix([[1/3, 1/6, 1/10, 1/15], [1/6, 2/15, 1/10, 8/105], [1/10, 1/10, 3/35, 1/14], [1/15, 8/105, 1/14, 4/63]]) \n",
"b:\n",
" Matrix([[b/6 - 1/2], [b/12 - 3/10], [b/20 - 1/5], [b/30 - 1/7]])\n",
"coeff: [b/2 - 1, -1, 0, 0]\n",
"approximation: -x**2*(1 - x) + x*(1 - x)*(b/2 - 1)\n",
"('solution u:', -b*x**2/2 + b*x/2 - x + 1)\n"
]
}
],
"source": [
"import sympy as sym\n",
"x, b = sym.symbols(\"x b\")\n",
"f = b\n",
"B = 1 - x**3\n",
"dBdx = sym.diff(B, x)\n",
"\n",
"# Compute basis functions and their derivatives\n",
"N = 3\n",
"psi = {0: [x**(i+1)*(1-x) for i in range(N+1)]}\n",
"psi[1] = [sym.diff(psi_i, x) for psi_i in psi[0]]\n",
"\n",
"def integrand_lhs(psi, i, j):\n",
" return psi[1][i]*psi[1][j]\n",
"\n",
"def integrand_rhs(psi, i):\n",
" return f*psi[0][i] - dBdx*psi[1][i]\n",
"\n",
"Omega = [0, 1]\n",
"\n",
"from src.varform1D import solver\n",
"u_bar, _ = solver(integrand_lhs, integrand_rhs, psi, Omega,\n",
" verbose=True, symbolic=True)\n",
"u = B + u_bar\n",
"print((\"solution u:\", sym.simplify(sym.expand(u))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The printout of `u` reads `-b*x**2/2 + b*x/2 - x + 1`. Note that\n",
"expanding `u`, before simplifying, is necessary in the present case to\n",
"get a compact, final expression with `sympy`. Doing `expand` before\n",
"`simplify` is a common strategy for simplifying expressions in\n",
"`sympy`. However, a non-expanded `u` might be preferable in other\n",
"cases - this depends on the problem in question.\n",
"\n",
"The exact solution ${u_{\\small\\mbox{e}}}(x)$ can be derived by some `sympy` code that\n",
"closely follows the examples in the section [Simple model problems and their solutions](#fem:deq:1D:models:simple). The idea is to integrate $-u''=b$ twice\n",
"and determine the integration constants from the boundary conditions:"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"# DO NOT RUN THIS CELL\n",
"C1, C2 = sym.symbols('C1 C2') # integration constants\n",
"f1 = sym.integrate(f, x) + C1\n",
"f2 = sym.integrate(f1, x) + C2\n",
"# Find C1 and C2 from the boundary conditions u(0)=0, u(1)=1\n",
"s = sym.solve([u_e.subs(x,0) - 1, u_e.subs(x,1) - 0], [C1, C2])\n",
"# Form the exact solution\n",
"u_e = -f2 + s[C1]*x + s[C2]\n",
"print(('analytical solution:', u_e))\n",
"print(('error:', sym.simplify(sym.expand(u - u_e))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The last line prints `0`, which is not surprising when\n",
"${u_{\\small\\mbox{e}}}(x)$ is a parabola and our approximate $u$ contains polynomials up to\n",
"degree 4. It suffices to have $N=1$, i.e., polynomials of degree\n",
"2, to recover the exact solution.\n",
"\n",
"We can play around with the code and test that with $f=Kx^p$, for\n",
"some constants $K$ and $p$,\n",
"the solution is a polynomial of degree $p+2$, and $N=p+1$ guarantees\n",
"that the approximate solution is exact.\n",
"\n",
"Although the symbolic code is capable of integrating many choices of $f(x)$,\n",
"the symbolic expressions for $u$ quickly become lengthy and non-informative,\n",
"so numerical integration in the code, and hence numerical answers,\n",
"have the greatest application potential."
]
}
],
"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",
"version": "3.7.4"
}
},
"nbformat": 4,
"nbformat_minor": 2
}