{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Time-dependent variational forms\n",
"<div id=\"ch:femtime\"></div>\n",
"\n",
"There are at least three different strategies for performing\n",
"a discretization in time:\n",
"\n",
"1. Use *finite differences* for time derivatives to arrive at\n",
" a recursive set of spatial problems that can be discretized by\n",
" the finite element method.\n",
"\n",
"2. Discretize in space by finite elements first, and then solve\n",
" the resulting system of ordinary differential equations (ODEs) by\n",
" some *standard library* for ODEs.\n",
"\n",
"3. Discretize in space and time simultaneously by space-time finite elements.\n",
"\n",
"With the first strategy, we discretize in time prior to the space\n",
"discretization, while the second strategy consists of doing exactly\n",
"the opposite. It should come as no surprise that in many situations\n",
"these two strategies end up in exactly the same systems to be solved, but\n",
"this is not always the case. Also the third approach often reproduces standard\n",
"finite difference schemes such as the Backward Euler and the Crank-Nicolson\n",
"schemes for lower-order elements, but offers an interesting framework for deriving higher-order\n",
"methods. In this chapter we shall be concerned with\n",
"the first strategy,\n",
"which is the most common strategy as it turns the time-dependent\n",
"PDE problem to a sequence of stationary problems for which efficient\n",
"finite element solution strategies often are available.\n",
"The second strategy would\n",
"naturally employ well-known ODE software,\n",
"which are available as user-friendly routines\n",
"in Python. However, these routines are presently not efficient enough\n",
"for PDE problems in 2D and 3D. The first strategy gives complete hands-on\n",
"control of the implementation and the computational efficiency\n",
"in time and space.\n",
"\n",
"We shall use a simple diffusion problem to illustrate the basic\n",
"principles of how a time-dependent PDE is solved by finite differences\n",
"in time and finite elements in space. Of course, instead of finite elements,\n",
"we may employ other types of basis functions, such as global polynomials.\n",
"Our model problem reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:eq\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial u}{\\partial t} = {\\alpha}\\nabla^2 u + f(\\boldsymbol{x}, t),\\quad\n",
"\\boldsymbol{x}\\in\\Omega,\\ t\\in (0,T],\n",
"\\label{fem:deq:diffu:eq} \\tag{1}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:ic\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(\\boldsymbol{x}, 0) = I(\\boldsymbol{x}),\\quad \\boldsymbol{x}\\in\\Omega,\n",
"\\label{fem:deq:diffu:ic} \\tag{2}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:bcN\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\frac{\\partial u}{\\partial n} = 0,\\quad \\boldsymbol{x}\\in\\partial\\Omega,\\ t\\in (0,T]\n",
"\\label{fem:deq:diffu:bcN} \\tag{3}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, $u(\\boldsymbol{x},t)$ is the unknown function, ${\\alpha}$ is a constant, and\n",
"$f(\\boldsymbol{x},t)$ and $I(\\boldsymbol{x})$ are given functions. We have assigned the particular\n",
"boundary condition ([3](#fem:deq:diffu:bcN)) to minimize\n",
"the details on handling boundary conditions in the finite element method.\n",
"\n",
"**Remark.** For systems of PDEs the strategy for discretization in time may have great impact on\n",
"overall efficiency and accuracy. The Navier-Stokes equations for \n",
"an incompressible Newtonian fluid is a prime example where many methods have been proposed\n",
"and where there are notable differences between the different methods. Furthermore, \n",
"the differences often depend significantly on the application. \n",
"Discretization in time *before* discretization in space allows for manipulations\n",
"of the equations and schemes that are very efficient compared to \n",
"schemes based on discretizing in space first. \n",
"The schemes are so-called operator-splitting schemes or projection based schemes. These schemes do, however, \n",
"suffer from loss of accuracy particularly in terms of errors associated with the boundaries. \n",
"The numerical error is caused by the splitting of the equations which leads to non-trivial splitting\n",
"of the boundary conditions. \n",
"\n",
"\n",
"# Discretization in time by a Forward Euler scheme\n",
"<div id=\"fem:deq:diffu:FE\"></div>\n",
"\n",
"The discretization strategy is to first apply a simple finite difference\n",
"scheme in time and derive a recursive set of spatially continuous PDE\n",
"problems, one at each time level. For each spatial PDE problem we can\n",
"set up a variational formulation and employ the finite element method\n",
"for solution.\n",
"\n",
"## Time discretization\n",
"\n",
"We can apply a finite difference method in time to ([1](#fem:deq:diffu:eq)).\n",
"First we need 'a mesh' in time, here taken as uniform with\n",
"mesh points $t_n = n\\Delta t$, $n=0,1,\\ldots,N_t$.\n",
"A Forward Euler scheme consists of sampling ([1](#fem:deq:diffu:eq))\n",
"at $t_n$ and approximating the time derivative by a forward\n",
"difference $[D_t^+ u]^n\\approx\n",
"(u^{n+1}-u^n)/\\Delta t$.\n",
"This approximation turns ([1](#fem:deq:diffu:eq))\n",
"into a differential equation that is discrete in time, but still\n",
"continuous in space.\n",
"With a finite difference operator notation we can write the\n",
"time-discrete problem as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:FE:eq:FEop\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"[D_t^+ u = {\\alpha}\\nabla^2 u + f]^n,\n",
"\\label{fem:deq:diffu:FE:eq:FEop} \\tag{4}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $n=1,2,\\ldots,N_t-1$.\n",
"Writing this equation out in detail and\n",
"isolating the unknown $u^{n+1}$ on the left-hand side, demonstrates that\n",
"the time-discrete problem is a recursive set of problems that are\n",
"continuous in space:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:FE:eq:unp1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u^{n+1} = u^n + \\Delta t \\left( {\\alpha}\\nabla^2 u^n + f(\\boldsymbol{x}, t_n)\\right)\n",
"\\label{fem:deq:diffu:FE:eq:unp1} \\tag{5}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given $u^0=I$, we can use ([5](#fem:deq:diffu:FE:eq:unp1)) to compute\n",
"$u^1,u^2,\\dots,u^{N_t}$.\n",
"\n",
"**More precise notation.**\n",
"\n",
"For absolute clarity in the various stages of the discretizations, we\n",
"introduce ${u_{\\small\\mbox{e}}}(\\boldsymbol{x},t)$ as the exact solution of the space-and time-continuous\n",
"partial differential equation ([1](#fem:deq:diffu:eq)) and\n",
"${u_{\\small\\mbox{e}}}^n(\\boldsymbol{x})$ as the time-discrete approximation, arising from the finite\n",
"difference method in time ([4](#fem:deq:diffu:FE:eq:FEop)).\n",
"More precisely, ${u_{\\small\\mbox{e}}}$ fulfills"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:eq:uex\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial {u_{\\small\\mbox{e}}}}{\\partial t} = {\\alpha}\\nabla^2 {u_{\\small\\mbox{e}}} + f(\\boldsymbol{x}, t)\n",
"\\label{fem:deq:diffu:eq:uex} \\tag{6},\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"while ${u_{\\small\\mbox{e}}}^{n+1}$, with a superscript,\n",
"is the solution of the time-discrete equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:FE:eq:uex:n\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"{u_{\\small\\mbox{e}}}^{n+1} = {u_{\\small\\mbox{e}}}^n + \\Delta t \\left( {\\alpha}\\nabla^2 {u_{\\small\\mbox{e}}}^n + f(\\boldsymbol{x}, t_n)\\right)\n",
"\\label{fem:deq:diffu:FE:eq:uex:n} \\tag{7}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The ${u_{\\small\\mbox{e}}}^{n+1}$ quantity is then discretized in space and approximated\n",
"by $u^{n+1}$.\n",
"\n",
"\n",
"\n",
"\n",
"## Space discretization\n",
"\n",
"We now introduce a finite element approximation to ${u_{\\small\\mbox{e}}}^n$ and ${u_{\\small\\mbox{e}}}^{n+1}$\n",
"in ([7](#fem:deq:diffu:FE:eq:uex:n)), where the coefficients depend on the\n",
"time level:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:femapprox:n\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"{u_{\\small\\mbox{e}}}^n \\approx u^n = \\sum_{j=0}^{N} c_j^{n}{\\psi}_j(\\boldsymbol{x}),\n",
"\\label{fem:deq:diffu:femapprox:n} \\tag{8}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:femapprox:np1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"{u_{\\small\\mbox{e}}}^{n+1} \\approx u^{n+1} = \\sum_{j=0}^{N} c_j^{n+1}{\\psi}_j(\\boldsymbol{x})\n",
"\\label{fem:deq:diffu:femapprox:np1} \\tag{9}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Note that, as before, $N$ denotes the number of degrees of freedom\n",
"in the spatial domain. The number of time points is denoted by $N_t$.\n",
"We define a space $V$ spanned by the basis functions $\\left\\{ {{\\psi}}_i \\right\\}_{i\\in{\\mathcal{I}_s}}$.\n",
"<!-- Also note that we use $u^n$ as the numerical solution we want -->\n",
"<!-- to compute in a program, while ${u_{\\small\\mbox{e}}}$ and ${u_{\\small\\mbox{e}}}^n$ are used when -->\n",
"<!-- we occasionally -->\n",
"<!-- need to refer to the exact solution and the time-discrete solution, -->\n",
"<!-- respectively. -->\n",
"\n",
"\n",
"## Variational forms\n",
"\n",
"A Galerkin method or a\n",
"weighted residual method with weighting functions $w_i$ can\n",
"now be formulated. We insert ([8](#fem:deq:diffu:femapprox:n)) and\n",
"([9](#fem:deq:diffu:femapprox:np1)) in\n",
"([7](#fem:deq:diffu:FE:eq:uex:n)) to obtain the residual"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"R = u^{n+1} - u^n - \\Delta t \\left( {\\alpha}\\nabla^2 u^n + f(\\boldsymbol{x}, t_n)\\right)\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The weighted residual principle,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_\\Omega Rw{\\, \\mathrm{d}x} = 0,\\quad \\forall w\\in W,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"results in"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_\\Omega\n",
"\\left\\lbrack\n",
"u^{n+1} - u^n - \\Delta t \\left( {\\alpha}\\nabla^2 u^n + f(\\boldsymbol{x}, t_n)\\right)\n",
"\\right\\rbrack w {\\, \\mathrm{d}x} =0, \\quad\\forall w \\in W{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From now on we use the Galerkin method so $W=V$.\n",
"Isolating the unknown $u^{n+1}$ on the left-hand side gives"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{\\Omega} u^{n+1}v{\\, \\mathrm{d}x} = \\int_{\\Omega}\n",
"\\left\\lbrack u^n + \\Delta t \\left( {\\alpha}\\nabla^2 u^n + f(\\boldsymbol{x}, t_n)\\right)\n",
"\\right\\rbrack v{\\, \\mathrm{d}x},\\quad \\forall v\\in V\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As usual in spatial finite element problems involving second-order\n",
"derivatives, we apply integration by parts on the term\n",
"$\\int (\\nabla^2 u^n)v{\\, \\mathrm{d}x}$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\int_{\\Omega}{\\alpha}(\\nabla^2 u^n)v {\\, \\mathrm{d}x} =\n",
"-\\int_{\\Omega}{\\alpha}\\nabla u^n\\cdot\\nabla v{\\, \\mathrm{d}x} +\n",
"\\int_{\\partial\\Omega}{\\alpha}\\frac{\\partial u^n}{\\partial n}v {\\, \\mathrm{d}x}\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The last term vanishes because we have the Neumann condition\n",
"$\\partial u^n/\\partial n=0$ for all $n$. Our discrete problem in\n",
"space and time then reads"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:FE:vf:u:np1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{\\Omega} u^{n+1}v{\\, \\mathrm{d}x} =\n",
"\\int_{\\Omega} u^n v{\\, \\mathrm{d}x} -\n",
"\\Delta t \\int_{\\Omega}{\\alpha}\\nabla u^n\\cdot\\nabla v{\\, \\mathrm{d}x} +\n",
"\\Delta t\\int_{\\Omega}f^n v{\\, \\mathrm{d}x},\\quad \\forall v\\in V{\\thinspace .}\n",
"\\label{fem:deq:diffu:FE:vf:u:np1} \\tag{10}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the variational formulation of our recursive set of spatial\n",
"problems.\n",
"\n",
"\n",
"**Nonzero Dirichlet boundary conditions.**\n",
"\n",
"As in stationary problems,\n",
"we can introduce a boundary function $B(\\boldsymbol{x},t)$ to take care\n",
"of nonzero Dirichlet conditions:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:femapprox:n:B\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"{u_{\\small\\mbox{e}}}^n \\approx u^n = B(\\boldsymbol{x},t_n) + \\sum_{j=0}^{N} c_j^{n}{\\psi}_j(\\boldsymbol{x}),\n",
"\\label{fem:deq:diffu:femapprox:n:B} \\tag{11}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:femapprox:np1:B\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"{u_{\\small\\mbox{e}}}^{n+1} \\approx u^{n+1} = B(\\boldsymbol{x},t_{n+1}) +\n",
"\\sum_{j=0}^{N} c_j^{n+1}{\\psi}_j(\\boldsymbol{x})\n",
"\\label{fem:deq:diffu:femapprox:np1:B} \\tag{12}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Notation for the solution at recent time levels\n",
"\n",
"In a program it is only necessary to have the two variables $u^{n+1}$\n",
"and $u^n$ at the same time at a given time step. It is therefore\n",
"unnatural to use the index $n$ in computer code. Instead a natural\n",
"variable naming is `u` for $u^{n+1}$, the new unknown, and `u_n` for\n",
"$u^n$, the solution at the previous time level. When we have several\n",
"preceding (already computed) time levels, it is natural to number them\n",
"like `u_nm1`, `u_nm2`, `u_nm3`, etc., backwards in time, corresponding to\n",
"$u^{n-1}$, $u^{n-2}$, and $u^{n-3}$. Essentially, this means a one-to-one\n",
"mapping of notation in mathematics and software, except for $u^{n+1}$.\n",
"We shall therefore, to make the distance between mathematics and code\n",
"as small as possible, often introduce just $u$ for $u^{n+1}$ in the\n",
"mathematical notation. Equation\n",
"([10](#fem:deq:diffu:FE:vf:u:np1)) with this new naming convention is\n",
"consequently expressed as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:FE:vf:u\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{\\Omega} u v{\\, \\mathrm{d}x} =\n",
"\\int_{\\Omega} u^{n} v{\\, \\mathrm{d}x} -\n",
"\\Delta t \\int_{\\Omega}{\\alpha}\\nabla u^{n}\\cdot\\nabla v{\\, \\mathrm{d}x} +\n",
"\\Delta t\\int_{\\Omega}f^n v{\\, \\mathrm{d}x}\n",
"{\\thinspace .}\n",
"\\label{fem:deq:diffu:FE:vf:u} \\tag{13}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This variational form can alternatively be expressed by the inner\n",
"product notation:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:FE:vf:u:short\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(u,v) = (u^{n},v) -\n",
"\\Delta t ({\\alpha}\\nabla u^{n},\\nabla v) +\n",
"\\Delta t (f^n, v)\n",
"{\\thinspace .}\n",
"\\label{fem:deq:diffu:FE:vf:u:short} \\tag{14}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To simplify the notation for the solution at recent previous time steps \n",
"and avoid notation like `u_nm1`, `u_nm2`, `u_nm3`, etc., we will let $u_1$ denote the solution at previous time step, \n",
"$u_2$ is the solution two time steps ago, etc. \n",
"\n",
"## Deriving the linear systems\n",
"\n",
"In the following, we adopt the previously introduced convention that\n",
"the unknowns $c_j^{n+1}$ are written as $c_j$, while the known $c_j^n$\n",
"from the previous time level is simply written as $c_{j}^n$. To\n",
"derive the equations for the new unknown coefficients $c_j$, we insert"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u = \\sum_{j=0}^{N}c_j{\\psi}_j(\\boldsymbol{x}),\\quad\n",
"u^{n} = \\sum_{j=0}^{N} c_{j}^n{\\psi}_j(\\boldsymbol{x})\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"in ([13](#fem:deq:diffu:FE:vf:u)) or ([14](#fem:deq:diffu:FE:vf:u:short)),\n",
"let the equation hold for all $v={\\psi}_i$, $i=0,\\ldots,N$,\n",
"and order the terms as matrix-vector products:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto1\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\sum_{j=0}^{N} ({\\psi}_i,{\\psi}_j) c_j =\n",
"\\sum_{j=0}^{N} ({\\psi}_i,{\\psi}_j) c_{j}^n\n",
"-\\Delta t \\sum_{j=0}^{N} (\\nabla{\\psi}_i,{\\alpha}\\nabla{\\psi}_j) c_{j}^n\n",
"+ \\Delta t (f^n,{\\psi}_i),\\quad i=0,\\ldots,N\n",
"{\\thinspace .}\n",
"\\label{_auto1} \\tag{15}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is a linear system $\\sum_j A_{i,j}c_j = b_i$ with"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"A_{i,j} = ({\\psi}_i,{\\psi}_j)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"b_i = \\sum_{j=0}^{N} ({\\psi}_i,{\\psi}_j) c_{j}^n\n",
"-\\Delta t \\sum_{j=0}^{N} (\\nabla{\\psi}_i,{\\alpha}\\nabla{\\psi}_j) c_{j}^n\n",
"+ \\Delta t (f^n,{\\psi}_i){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It is instructive and convenient for implementations to write the linear\n",
"system on the form"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto2\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"Mc = Mc_1 - \\Delta t Kc_1 + \\Delta t f,\n",
"\\label{_auto2} \\tag{16}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"M &= \\{M_{i,j}\\},\\quad M_{i,j}=({\\psi}_i,{\\psi}_j),\\quad i,j\\in{\\mathcal{I}_s},\\\\\n",
"K &= \\{K_{i,j}\\},\\quad K_{i,j}=(\\nabla{\\psi}_i,{\\alpha}\\nabla{\\psi}_j),\n",
"\\quad i,j\\in{\\mathcal{I}_s},\\\\\n",
"f &= \\{f_i\\},\\quad f_i=(f(\\boldsymbol{x},t_n),{\\psi}_i),\\quad i\\in{\\mathcal{I}_s},\\\\\n",
"c &= \\{c_i\\},\\quad i\\in{\\mathcal{I}_s},\\\\\n",
"c_1 &= \\{c_{i}^n\\},\\quad i\\in{\\mathcal{I}_s}\n",
"{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We realize that $M$ is the matrix arising from a term with the\n",
"zero-th derivative of $u$, and called the mass matrix, while $K$ is\n",
"the matrix arising from a Laplace term $\\nabla^2 u$. The $K$ matrix\n",
"is often known as the *stiffness matrix*. (The terms mass and stiffness\n",
"stem from the early days of finite elements when applications to\n",
"vibrating structures dominated. The mass matrix arises from the\n",
"mass times acceleration term in Newton's second law, while the stiffness\n",
"matrix arises from the elastic forces (the \"stiffness\") in that law.\n",
"The mass and stiffness\n",
"matrix appearing in a diffusion have slightly different mathematical\n",
"formulas compared to the classic structure problem.)\n",
"\n",
"**Remark.** The mathematical symbol $f$ has two meanings, either the\n",
"function $f(\\boldsymbol{x},t)$ in the PDE or the $f$ vector in the linear system\n",
"to be solved at each time level.\n",
"\n",
"## Computational algorithm\n",
"\n",
"We observe that $M$ and $K$ can be precomputed so that we can avoid\n",
"computing the matrix entries at every time level. Instead, some\n",
"matrix-vector multiplications will produce the linear system to be solved.\n",
"The computational algorithm has the following steps:\n",
"\n",
"1. Compute $M$ and $K$.\n",
"\n",
"2. Initialize $u^0$ by interpolation or projection\n",
"\n",
"3. For $n=1,2,\\ldots,N_t$:\n",
"\n",
"a. compute $b = Mc_1 - \\Delta t Kc_1 + \\Delta t f$\n",
"\n",
"b. solve $Mc = b$\n",
"\n",
"c. set $c_1 = c$\n",
"\n",
"\n",
"In case of finite element basis functions, interpolation of the\n",
"initial condition at the nodes means $c_{j}^n = I(\\boldsymbol{x}_j)$. Otherwise\n",
"one has to solve the linear system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sum_j{\\psi}_j(\\boldsymbol{x}_i)c_{j}^n = I(\\boldsymbol{x}_i),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\boldsymbol{x}_i$ denotes an interpolation point. Projection\n",
"(or Galerkin's method) implies solving a linear system with $M$ as\n",
"coefficient matrix:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\sum_j M_{i,j}c_{j}^n = (I,{\\psi}_i),\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example using cosinusoidal basis functions\n",
"<div id=\"fem:deq:diffu:FE:cosex\"></div>\n",
"\n",
"Let us go through a computational example and demonstrate the\n",
"algorithm from the previous section. We consider a 1D problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:pde1D:eq\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\frac{\\partial u}{\\partial t} = {\\alpha}\\frac{\\partial^2 u}{\\partial x^2},\\quad\n",
"x\\in (0,L),\\ t\\in (0,T],\n",
"\\label{fem:deq:diffu:pde1D:eq} \\tag{17}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:pde1D:ic\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(x, 0) = A\\cos(\\pi x/L) + B\\cos(10\\pi x/L),\\quad x\\in[0,L],\n",
"\\label{fem:deq:diffu:pde1D:ic} \\tag{18}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:pde1D:bcN\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\frac{\\partial u}{\\partial x} = 0,\\quad x=0,L,\\ t\\in (0,T]\n",
"\\label{fem:deq:diffu:pde1D:bcN} \\tag{19}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We use a Galerkin method with basis functions"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"{\\psi}_i = \\cos(i\\pi x/L){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"These basis functions fulfill ([19](#fem:deq:diffu:pde1D:bcN)), which is\n",
"not a requirement (there are no Dirichlet conditions in this problem),\n",
"but helps to make the approximation good.\n",
"\n",
"Since the initial condition ([18](#fem:deq:diffu:pde1D:ic)) lies in the\n",
"space $V$ where we seek the approximation, we know that a Galerkin or\n",
"least squares approximation of the initial condition becomes exact.\n",
"Therefore, the initial condition can be expressed as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"c_{1}^n=A,\\quad c_{10}^n=B,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"while $c_{i}^n=0$ for $i\\neq 1,10$.\n",
"\n",
"The $M$ and $K$ matrices are easy to compute since the basis functions\n",
"are orthogonal on $[0,L]$. Hence, we\n",
"only need to compute the diagonal entries. We get"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"M_{i,i} = \\int_0^L cos^2(i x \\pi/L) {\\, \\mathrm{d}x},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which is computed as"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import sympy as sym\n",
"x, L = sym.symbols('x L')\n",
"i = sym.symbols('i', integer=True)\n",
"sym.integrate(sym.cos(i*x*sym.pi/L)**2, (x,0,L))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which means $L$ if $i=0$ and $L/2$ otherwise. Similarly,\n",
"the diagonal entries of the $K$ matrix are computed as"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"sym.integrate(sym.diff(cos(i*x*sym.pi/L),x)**2, (x,0,L))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"so"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"M_{0,0}=L,\\quad M_{i,i}=L/2,\\ i>0,\\quad K_{0,0}=0,\\quad K_{i,i}=\\frac{\\pi^2 i^2}{2L},\\ i>0{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The equation system becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"Lc_0 &= Lc_{0}^0 - \\Delta t \\cdot 0\\cdot c_{0}^0,\\\\\n",
"\\frac{L}{2}c_i &= \\frac{L}{2}c_{i}^n - \\Delta t\n",
"\\frac{\\pi^2 i^2}{2L} c_{i}^n,\\quad i>0{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first equation leads to $c_0=0$ for any $n$ since we start with $c_{0}^0=0$ and $K_{0,0}=0$. \n",
"The others imply"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"c_i = (1-\\Delta t (\\frac{\\pi i}{L})^2) c_{i}^n{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With the notation $c^n_i$ for $c_i$ at the $n$-th time level, we can apply\n",
"the relation above recursively and get"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"c^n_i = (1-\\Delta t (\\frac{\\pi i}{L})^2)^n c^0_i{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Since only two of the coefficients are nonzero at time $t=0$, we have\n",
"the closed-form discrete solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u^n_i = A(1-\\Delta t (\\frac{\\pi}{L})^2)^n \\cos(\\pi x/L)\n",
"+ B(1-\\Delta t (\\frac{10\\pi }{L})^2)^n \\cos(10\\pi x/L){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"\n",
"\n",
"\n",
"\n",
"# Discretization in time by a Backward Euler scheme\n",
"<div id=\"fem:deq:diffu:BE\"></div>\n",
"\n",
"## Time discretization\n",
"\n",
"The Backward Euler scheme in time applied to our diffusion problem\n",
"can be expressed as follows using the finite difference operator notation:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"[D_t^- u = {\\alpha}\\nabla^2 u + f(\\boldsymbol{x}, t)]^n\n",
"{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here $[D_t^- u]^n\\approx (u^{n}-u^{n-1})/\\Delta t$.\n",
"Written out, and collecting the unknown $u^n$ on the left-hand side\n",
"and all the known terms on the right-hand side,\n",
"the time-discrete differential equation becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:BE:eq:un\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u^{n} - \\Delta t {\\alpha}\\nabla^2 u^n =\n",
"u^{n-1} + \\Delta t f(\\boldsymbol{x}, t_{n})\n",
"\\label{fem:deq:diffu:BE:eq:un} \\tag{22}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"From equation ([22](#fem:deq:diffu:BE:eq:un)) we can compute\n",
"$u^1,u^2,\\dots,u^{N_t}$,\n",
"if we have a start $u^0=I$ from the initial condition.\n",
"However, ([22](#fem:deq:diffu:BE:eq:un)) is a partial differential\n",
"equation in space and needs a solution method based on discretization\n",
"in space. For this purpose we use an expansion as in\n",
"([8](#fem:deq:diffu:femapprox:n))-([9](#fem:deq:diffu:femapprox:np1)).\n",
"\n",
"## Variational forms\n",
"\n",
"Inserting ([8](#fem:deq:diffu:femapprox:n))-([9](#fem:deq:diffu:femapprox:np1))\n",
"in ([22](#fem:deq:diffu:BE:eq:un)), multiplying by any $v\\in V$\n",
"(or ${\\psi}_i\\in V$),\n",
"and integrating by parts, as we did in the Forward Euler case, results\n",
"in the variational form"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:BE:vf:u:n\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{\\Omega} \\left( u^{n}v\n",
"+ \\Delta t {\\alpha}\\nabla u^n\\cdot\\nabla v\\right){\\, \\mathrm{d}x}\n",
"= \\int_{\\Omega} u^{n-1} v{\\, \\mathrm{d}x} +\n",
"\\Delta t\\int_{\\Omega}f^n v{\\, \\mathrm{d}x},\\quad\\forall v\\in V\n",
"\\label{fem:deq:diffu:BE:vf:u:n} \\tag{23}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Expressed with $u$ for the unknown $u^n$ and $u^{n}$ for the previous\n",
"time level, as we have done before, the variational form becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:BE:vf:u\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int_{\\Omega} \\left( uv\n",
"+ \\Delta t {\\alpha}\\nabla u\\cdot\\nabla v\\right){\\, \\mathrm{d}x}\n",
"= \\int_{\\Omega} u^{n} v{\\, \\mathrm{d}x} +\n",
"\\Delta t\\int_{\\Omega}f^n v{\\, \\mathrm{d}x},\n",
"\\label{fem:deq:diffu:BE:vf:u} \\tag{24}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"or with the more compact inner product notation,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:BE:vf:u:short\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(u,v) + \\Delta t ({\\alpha}\\nabla u,\\nabla v)\n",
"= (u^{n},v) +\n",
"\\Delta t (f^n,v)\n",
"\\label{fem:deq:diffu:BE:vf:u:short} \\tag{25}\n",
"{\\thinspace .}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Linear systems\n",
"\n",
"Inserting $u=\\sum_j c_j{\\psi}_i$ and $u^{n}=\\sum_j c_{j}^n{\\psi}_i$,\n",
"and choosing $v$ to be the basis functions ${\\psi}_i\\in V$,\n",
"$i=0,\\ldots,N$, together with doing some algebra, lead\n",
"to the following linear system to be\n",
"solved at each time level:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:BE:vf:linsys\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"(M + \\Delta t K)c = Mc_1 + \\Delta t f,\n",
"\\label{fem:deq:diffu:BE:vf:linsys} \\tag{26}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $M$, $K$, and $f$ are as in the Forward Euler case\n",
"and we use the previously introduced notation $c = \\{c_i\\}$ and $c_1 = \\{c_{i}^n\\}$.\n",
"\n",
"\n",
"This time we really have to solve a linear system at each time level.\n",
"The computational algorithm goes as follows.\n",
"\n",
"1. Compute $M$, $K$, and $A=M + \\Delta t K$\n",
"\n",
"2. Initialize $u^0$ by interpolation or projection\n",
"\n",
"3. For $n=1,2,\\ldots,N_t$:\n",
"\n",
"a. compute $b = Mc_1 + \\Delta t f$\n",
"\n",
"b. solve $Ac = b$\n",
"\n",
"c. set $c_1 = c$\n",
"\n",
"\n",
"In case of finite element basis functions, interpolation of the\n",
"initial condition at the nodes means $c_{j}^n = I(\\boldsymbol{x}_j)$. Otherwise\n",
"one has to solve the linear system $\\sum_j{\\psi}_j(\\boldsymbol{x}_i)c_j =\n",
"I(\\boldsymbol{x}_i)$, where $\\boldsymbol{x}_i$ denotes an interpolation point. Projection\n",
"(or Galerkin's method) implies solving a linear system with $M$ as\n",
"coefficient matrix: $\\sum_j M_{i,j}c_{j}^n = (I,{\\psi}_i)$,\n",
"$i\\in{\\mathcal{I}_s}$.\n",
"\n"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Dirichlet boundary conditions\n",
"<div id=\"fem:deq:diffu:Dirichlet\"></div>\n",
"\n",
"\n",
"Suppose now that the boundary condition ([3](#fem:deq:diffu:bcN)) is\n",
"replaced by a mixed Neumann and Dirichlet condition,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto3\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u(\\boldsymbol{x},t) = u_0(\\boldsymbol{x},t),\\quad \\boldsymbol{x}\\in\\partial\\Omega_D,\n",
"\\label{_auto3} \\tag{29}\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",
"-{\\alpha}\\frac{\\partial}{\\partial n} u(\\boldsymbol{x},t) = g(\\boldsymbol{x},t),\\quad\n",
" \\boldsymbol{x}\\in\\partial{\\Omega}_N{\\thinspace .}\n",
"\\label{_auto4} \\tag{30}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Using a Forward Euler discretization in time, the variational\n",
"form at a time level becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto5\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int\\limits_\\Omega u^{n+1}v{\\, \\mathrm{d}x} =\n",
"\\int\\limits_\\Omega (u^n - \\Delta t{\\alpha}\\nabla u^n\\cdot\\nabla v){\\, \\mathrm{d}x} +\n",
"\\Delta t\\int\\limits_\\Omega fv {\\, \\mathrm{d}x} -\n",
"\\Delta t\\int\\limits_{\\partial\\Omega_N} gv{\\, \\mathrm{d}s},\\quad \\forall v\\in V{\\thinspace .}\n",
"\\label{_auto5} \\tag{31}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Boundary function\n",
"\n",
"\n",
"\n",
"The Dirichlet condition $u=u_0$ at $\\partial\\Omega_D$ can be incorporated\n",
"through a boundary function $B(\\boldsymbol{x})=u_0(\\boldsymbol{x})$ and demanding that the basis functions ${\\psi}_j=0$\n",
"at $\\partial\\Omega_D$. The expansion for $u^n$ is written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u^n(\\boldsymbol{x}) = u_0(\\boldsymbol{x},t_n) + \\sum_{j\\in{\\mathcal{I}_s}}c_j^n{\\psi}_j(\\boldsymbol{x}){\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Inserting this expansion in the variational formulation and letting it\n",
"hold for all test functions $v\\in V$, i.e., all basis functions ${\\psi}_i$ leads to the linear system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\sum_{j\\in{\\mathcal{I}_s}} \\left(\\int\\limits_\\Omega {\\psi}_i{\\psi}_j{\\, \\mathrm{d}x}\\right)\n",
"c^{n+1}_j &= \\sum_{j\\in{\\mathcal{I}_s}}\n",
"\\left(\\int\\limits_\\Omega\\left( {\\psi}_i{\\psi}_j -\n",
"\\Delta t{\\alpha}\\nabla {\\psi}_i\\cdot\\nabla{\\psi}_j\\right){\\, \\mathrm{d}x}\\right) c_j^n - \\\\\n",
"&\\quad \\int\\limits_\\Omega\\left( u_0(\\boldsymbol{x},t_{n+1}) - u_0(\\boldsymbol{x},t_n)\n",
"+ \\Delta t{\\alpha}\\nabla u_0(\\boldsymbol{x},t_n)\\cdot\\nabla\n",
"{\\psi}_i\\right){\\, \\mathrm{d}x} \\\\\n",
"& \\quad + \\Delta t\\int\\limits_\\Omega f{\\psi}_i{\\, \\mathrm{d}x} -\n",
"\\Delta t\\int\\limits_{\\partial\\Omega_N} g{\\psi}_i{\\, \\mathrm{d}s},\n",
"\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Finite element basis functions\n",
"\n",
"When using finite elements, each basis function ${\\varphi}_i$ is associated\n",
"with a node $\\boldsymbol{x}_{i}$. We have a collection of nodes\n",
"$\\{\\boldsymbol{x}_i\\}_{i\\in{I_b}}$ on the boundary $\\partial\\Omega_D$.\n",
"Suppose $U_k^n$ is the known\n",
"Dirichlet value at $\\boldsymbol{x}_{k}$ at time $t_n$ ($U_k^n=u_0(\\boldsymbol{x}_{k},t_n)$).\n",
"The appropriate boundary function is then"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"B(\\boldsymbol{x},t_n)=\\sum_{j\\in{I_b}} U_j^n{\\varphi}_j{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The unknown coefficients $c_j$ are associated with the rest of the nodes,\n",
"which have numbers $\\nu(i)$, $i\\in{\\mathcal{I}_s} = \\{0,\\ldots,N\\}$. The basis\n",
"functions of $V$ are chosen as ${\\psi}_i = {\\varphi}_{\\nu(i)}$, $i\\in{\\mathcal{I}_s}$,\n",
"and all of these vanish at the boundary nodes as they should.\n",
"The expansion for $u^{n+1}$ and $u^n$ become"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"u^n &= \\sum_{j\\in{I_b}} U_j^n{\\varphi}_j + \\sum_{j\\in{\\mathcal{I}_s}}c_{j}^n{\\varphi}_{\\nu(j)},\\\\\n",
"u^{n+1} &= \\sum_{j\\in{I_b}} U_j^{n+1}{\\varphi}_j +\n",
"\\sum_{j\\in{\\mathcal{I}_s}}c_{j}{\\varphi}_{\\nu(j)}{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The equations for the unknown coefficients $\\left\\{ {c}_j \\right\\}_{j\\in{\\mathcal{I}_s}}$ become"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\sum_{j\\in{\\mathcal{I}_s}} \\left(\\int\\limits_\\Omega {\\varphi}_i{\\varphi}_j{\\, \\mathrm{d}x}\\right)\n",
"c_j &= \\sum_{j\\in{\\mathcal{I}_s}}\n",
"\\left(\\int\\limits_\\Omega\\left( {\\varphi}_i{\\varphi}_j -\n",
"\\Delta t{\\alpha}\\nabla {\\varphi}_i\\cdot\\nabla{\\varphi}_j\\right){\\, \\mathrm{d}x}\\right) c_{j}^n\n",
"- \\\\\n",
"&\\quad \\sum_{j\\in{I_b}}\\int\\limits_\\Omega\\left( {\\varphi}_i{\\varphi}_j(U_j^{n+1} - U_j^n)\n",
"+ \\Delta t{\\alpha}\\nabla {\\varphi}_i\\cdot\\nabla\n",
"{\\varphi}_jU_j^n\\right){\\, \\mathrm{d}x} \\\\\n",
"&\\quad + \\Delta t\\int\\limits_\\Omega f{\\varphi}_i{\\, \\mathrm{d}x} -\n",
"\\Delta t\\int\\limits_{\\partial\\Omega_N} g{\\varphi}_i{\\, \\mathrm{d}s},\n",
"\\quad i\\in{\\mathcal{I}_s}{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Modification of the linear system\n",
"\n",
"Instead of introducing a boundary function $B$ we can work with\n",
"basis functions associated with all the nodes and incorporate the\n",
"Dirichlet conditions by modifying the linear system.\n",
"Let ${\\mathcal{I}_s}$ be the index set that counts all the nodes:\n",
"$\\{0,1,\\ldots,N=N_n-1\\}$. The\n",
"expansion for $u^n$ is then $\\sum_{j\\in{\\mathcal{I}_s}}c^n_j{\\varphi}_j$ and the\n",
"variational form becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"\\sum_{j\\in{\\mathcal{I}_s}} \\left(\\int\\limits_\\Omega {\\varphi}_i{\\varphi}_j{\\, \\mathrm{d}x}\\right)\n",
"c_j &= \\sum_{j\\in{\\mathcal{I}_s}}\n",
"\\left(\\int\\limits_\\Omega\\left( {\\varphi}_i{\\varphi}_j -\n",
"\\Delta t{\\alpha}\\nabla {\\varphi}_i\\cdot\\nabla{\\varphi}_j\\right){\\, \\mathrm{d}x}\\right) c_{1,j}\n",
" \\\\\n",
"&\\quad + \\Delta t\\int\\limits_\\Omega f{\\varphi}_i{\\, \\mathrm{d}x} -\n",
"\\Delta t\\int\\limits_{\\partial\\Omega_N} g{\\varphi}_i{\\, \\mathrm{d}s}{\\thinspace .}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We introduce the matrices $M$ and $K$ with entries\n",
"$M_{i,j}=\\int\\limits_\\Omega{\\varphi}_i{\\varphi}_j{\\, \\mathrm{d}x}$ and\n",
"$K_{i,j}=\\int\\limits_\\Omega{\\alpha}\\nabla{\\varphi}_i\\cdot\\nabla{\\varphi}_j{\\, \\mathrm{d}x}$,\n",
"respectively.\n",
"In addition, we define the vectors $c$, $c_1$, and $f$ with\n",
"entries $c_i$, $c_{1,i}$, and\n",
"$\\int\\limits_\\Omega f{\\varphi}_i{\\, \\mathrm{d}x} - \\int\\limits_{\\partial\\Omega_N}g{\\varphi}_i{\\, \\mathrm{d}s}$, respectively.\n",
"The equation system can then be written as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto6\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"Mc = Mc_1 - \\Delta t Kc_1 + \\Delta t f{\\thinspace .}\n",
"\\label{_auto6} \\tag{32}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"When $M$, $K$, and $f$ are assembled without paying attention to\n",
"Dirichlet boundary conditions, we need to replace equation $k$\n",
"by $c_k=U_k$ for $k$ corresponding to all boundary nodes ($k\\in{I_b}$).\n",
"The modification of $M$ consists in setting $M_{k,j}=0$, $j\\in{\\mathcal{I}_s}$, and\n",
"the $M_{k,k}=1$. Alternatively, a modification that preserves\n",
"the symmetry of $M$ can be applied. At each time level one forms\n",
"$b = Mc_1 - \\Delta t Kc_1 + \\Delta t f$ and sets $b_k=U^{n+1}_k$,\n",
"$k\\in{I_b}$, and solves the system $Mc=b$.\n",
"\n",
"In case of a Backward Euler method, the system becomes\n",
"([26](#fem:deq:diffu:BE:vf:linsys)). We can write the system\n",
"as $Ac=b$, with $A=M + \\Delta t K$ and $b = Mc_1 + f$.\n",
"Both $M$ and $K$ needs to be modified because of Dirichlet\n",
"boundary conditions, but the diagonal entries in $K$ should be\n",
"set to zero and those in $M$ to unity. In this way, for $k\\in{I_b}$ we\n",
"have $A_{k,k}=1$.\n",
"The right-hand side must read $b_k=U^n_k$ for $k\\in{I_b}$ (assuming\n",
"the unknown is sought at time level $t_n$).\n",
"\n",
"## Example: Oscillating Dirichlet boundary condition\n",
"<div id=\"fem:deq:diffu:Dirichlet:ex\"></div>\n",
"\n",
"We shall address the one-dimensional initial-boundary value problem"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:Dirichlet:ex:pde\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"u_t = ({\\alpha} u_x)_x + f,\\quad x\\in\\Omega =[0,L],\\ t\\in (0,T],\n",
"\\label{fem:deq:diffu:Dirichlet:ex:pde} \\tag{33} \n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:Dirichlet:ex:uic\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(x,0) = 0,\\quad x\\in\\Omega,\n",
"\\label{fem:deq:diffu:Dirichlet:ex:uic} \\tag{34}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:Dirichlet:ex:uL\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u(0,t) = a\\sin\\omega t,\\quad t\\in (0,T],\n",
"\\label{fem:deq:diffu:Dirichlet:ex:uL} \\tag{35}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:Dirichlet:ex:uR\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"u_x(L,t) = 0,\\quad t\\in (0,T]{\\thinspace .}\n",
"\\label{fem:deq:diffu:Dirichlet:ex:uR} \\tag{36}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"A physical interpretation may be that $u$ is the temperature\n",
"deviation from a constant mean temperature in a body $\\Omega$\n",
"that is subject to an oscillating temperature (e.g., day and\n",
"night, or seasonal, variations) at $x=0$.\n",
"\n",
"We use a Backward Euler scheme in time and P1 elements of\n",
"constant length $h$ in space.\n",
"Incorporation of the Dirichlet condition at $x=0$ through\n",
"modifying the linear system at each time level means that we\n",
"carry out the computations as explained in the section [Discretization in time by a Backward Euler scheme](#fem:deq:diffu:BE) and get a system ([26](#fem:deq:diffu:BE:vf:linsys)).\n",
"The $M$ and $K$ matrices computed without paying attention to\n",
"Dirichlet boundary conditions become"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto7\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"M = \\frac{h}{6}\n",
"\\left(\n",
"\\begin{array}{cccccccccc}\n",
"2 & 1 & 0\n",
"&\\cdots & \\cdots & \\cdots & \\cdots & \\cdots & 0 \\\\\n",
"1 & 4 & 1 & \\ddots & & & & & \\vdots \\\\\n",
"0 & 1 & 4 & 1 &\n",
"\\ddots & & & & \\vdots \\\\\n",
"\\vdots & \\ddots & & \\ddots & \\ddots & 0 & & & \\vdots \\\\\n",
"\\vdots & & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\\n",
"\\vdots & & & 0 & 1 & 4 & 1 & \\ddots & \\vdots \\\\\n",
"\\vdots & & & & \\ddots & \\ddots & \\ddots &\\ddots & 0 \\\\\n",
"\\vdots & & & & &\\ddots & 1 & 4 & 1 \\\\\n",
"0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots & 0 & 1 & 2\n",
"\\end{array}\n",
"\\right)\n",
"\\label{_auto7} \\tag{37}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto8\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"K = \\frac{{\\alpha}}{h}\n",
"\\left(\n",
"\\begin{array}{cccccccccc}\n",
"1 & -1 & 0 &\\cdots & \\cdots & \\cdots & \\cdots & \\cdots & 0 \\\\\n",
"-1 & 2 & -1 & \\ddots & & & & & \\vdots \\\\\n",
"0 & -1 & 2 & -1 & \\ddots & & & & \\vdots \\\\\n",
"\\vdots & \\ddots & & \\ddots & \\ddots & 0 & & & \\vdots \\\\\n",
"\\vdots & & \\ddots & \\ddots & \\ddots & \\ddots & \\ddots & & \\vdots \\\\\n",
"\\vdots & & & 0 & -1 & 2 & -1 & \\ddots & \\vdots \\\\\n",
"\\vdots & & & & \\ddots & \\ddots & \\ddots &\\ddots & 0 \\\\\n",
"\\vdots & & & & &\\ddots & -1 & 2 & -1 \\\\\n",
"0 &\\cdots & \\cdots &\\cdots & \\cdots & \\cdots & 0 & -1 & 1\n",
"\\end{array}\n",
"\\right)\n",
"\\label{_auto8} \\tag{38}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The right-hand side of the variational form contains no source term ($f$) and no boundary term from the\n",
"integration by parts ($u_x=0$ at $x=L$ and we compute as if $u_x=0$ at\n",
"$x=0$ too) and we are therefore left with $Mc_1$. However, we must incorporate the Dirichlet boundary\n",
"condition $c_0=a\\sin\\omega t_n$. Let us assume that our numbering of nodes is such that\n",
"${\\mathcal{I}_s} = \\{0,1,\\ldots,N=N_n-1\\}$.\n",
"The Dirichlet condition can then be incorporated\n",
" by ensuring that this is the\n",
"first equation in the linear system.\n",
"To this end,\n",
"the first row in $K$ and $M$ is set to zero, but the diagonal\n",
"entry $M_{0,0}$ is set to 1. The right-hand side is $b=Mc_1$,\n",
"and we set $b_0 = a\\sin\\omega t_n$.\n",
"We can write the complete linear system as"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto9\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"c_0 = a\\sin\\omega t_n,\n",
"\\label{_auto9} \\tag{39}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto10\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\frac{h}{6}(c_{i-1} + 4c_i + c_{i+1}) + \\Delta t\\frac{{\\alpha}}{h}(-c_{i-1}\n",
"+2c_i - c_{i+1}) = \\frac{h}{6}(c_{1,i-1} + 4c_{1,i} + c_{1,i+1}),\n",
"\\label{_auto10} \\tag{40}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\qquad i=1,\\ldots,N_n-2,\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"_auto11\"></div>\n",
"\n",
"$$\n",
"\\begin{equation} \n",
"\\frac{h}{6}(c_{i-1} + 2c_i) + \\Delta t\\frac{{\\alpha}}{h}(-c_{i-1}\n",
"+c_i) = \\frac{h}{6}(c_{1,i-1} + 2c_{1,i}),\n",
"\\label{_auto11} \\tag{41}\n",
"\\end{equation}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\qquad i=N_n-1{\\thinspace .}\\nonumber\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The Dirichlet boundary condition can alternatively be implemented\n",
"through a boundary function $B(x,t)=a\\sin\\omega t\\,{\\varphi}_0(x)$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u^n(x) = a\\sin\\omega t_n\\,{\\varphi}_0(x) +\n",
"\\sum_{j\\in{\\mathcal{I}_s}} c_j{\\varphi}_{\\nu(j)}(x),\\quad\n",
"\\nu(j) = j+1{\\thinspace .}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Now, $N=N_n-2$ and the $c$ vector contains values of $u$ at nodes\n",
"$1,2,\\ldots,N_n-1$. The right-hand side gets a contribution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"<!-- Equation labels as ordinary links -->\n",
"<div id=\"fem:deq:diffu:Dirichlet:ex:bterm\"></div>\n",
"\n",
"$$\n",
"\\begin{equation}\n",
"\\int\\limits_0^L \\left(\n",
"a(\\sin\\omega t_n - \\sin\\omega t_{n-1}){\\varphi}_0{\\varphi}_i\n",
"- \\Delta t{\\alpha} a\\sin\\omega t_n\\nabla{\\varphi}_0\\cdot\\nabla{\\varphi}_i\\right){\\, \\mathrm{d}x}\n",
"{\\thinspace .}\n",
"\\label{fem:deq:diffu:Dirichlet:ex:bterm} \\tag{42}\n",
"\\end{equation}\n",
"$$"
]
}
],
"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.8.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}