{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$\\newcommand{mb}[1]{\\mathbf{#1}}$\n",
"\n",
"\n",
"# Partial differential equations and finite difference methods.\n",
"\n",
" \n",
"**Anne Kværnø**\n",
"\n",
"Date: **Dec 3, 2018**\n",
"\n",
"# Introduction\n",
"In this note the finite difference method for solving partial differential equations\n",
"(PDEs) will be presented. \n",
"\n",
"Roughly speaking, a finite difference method is composed by the following steps: \n",
"1. Discretize the domain on which the equation is defined. \n",
"\n",
"2. On each grid point, replace the derivatives with an approximation, using the values in neighbouring grid points. \n",
"\n",
"3. Replace the exact solutions by its approximations.\n",
"\n",
"4. Solve the resulting system of equations. \n",
"\n",
"We will first see how to find approximations to the\n",
"derivative of a function, and then how this can be used to solve boundary value\n",
"problems like"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u'' + p(x) u' + q(x) u = r(x), \\qquad a \\leq x \\leq b, \\qquad u(a)=u_a, \\quad\n",
"u(b)=u_b,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"and time dependent partial differential equations with the heat equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u_t = u_{xx}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"as the example of choice. The technique described here is however applicable several other PDEs, so please concentrate on understanding the underlying idea. \n",
"\n",
"# Numerical differentiation.\n",
"This is the main tool for finite difference methods.\n",
"\n",
"Given a sufficiently smooth function $f$. How can we find an approximation to\n",
"$f'(x)$ or $f''(x)$ in some given point $x$, just by using evaluation of the\n",
"function itself?\n",
"\n",
"[The derivative of of\n",
"$f$](https://wiki.math.ntnu.no/tma4100/tema/differentiation?&#definisjonen_av_den_deriverte_gitt_som_en_grenseverd)\n",
"is defined by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f'(x) = \\lim_{h\\rightarrow 0} \\frac{f(x+h)-f(x)}{h}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Given a sufficiently small value of $h$, the right hand side can be used an\n",
"approximation to $f'(x)$. A small collection of the most used approximations to\n",
"$f'(x)$ is:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f'(x) \\approx \\left\\{\n",
" \\begin{array}{ll}\n",
" \\displaystyle \\frac{f(x+h)-f(x)}{h}, \\qquad & \\text{Forward difference,} \\\\ \n",
" \\displaystyle \\frac{f(x)-f(x-h)}{h}, & \\text{Bakward difference,} \\\\ \n",
" \\displaystyle \\frac{f(x+h)-f(x-h)}{2h}, & \\text{Central difference.}\n",
" \\end{array} \\right.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first one is taken directly from the definition, so is the second, and the\n",
"third is just the mean of the first two. A common approximation to the second\n",
"derivative is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"f''(x) \\approx \\frac{f(x+h)-2f(x)+f(x-h)}{h^2}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Numerical example 1:**\n",
"Test the method on the function $f(x)=\\sin(x)$ at the point $x=\\frac{\\pi}{4}$.\n",
"Compare with the exact derivative. Try different step sizes, e.g. $h=0.1, h=0.01, h=0.001$.\n",
"Notice how the error in each case change with $h$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"\n",
"from numpy import * \n",
"from scipy.sparse import diags\t # Greate diagonal matrices\n",
"from scipy.linalg import solve\t # Solve linear systems\n",
"from matplotlib.pyplot import * \t\n",
"from mpl_toolkits.mplot3d import Axes3D # For 3-d plot\n",
"from matplotlib import cm \n",
"newparams = {'figure.figsize': (8.0, 4.0), 'axes.grid': True,\n",
" 'lines.markersize': 8, 'lines.linewidth': 2,\n",
" 'font.size': 14}\n",
"rcParams.update(newparams)\n",
"from numpy import *"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Numerical differentiation\n",
"\n",
"# Forward difference\n",
"def diff_forward(f, x, h=0.1):\n",
" return (f(x+h)-f(x))/h\n",
"\n",
"# Backward difference\n",
"def diff_backward(f, x, h=0.1):\n",
" return (f(x)-f(x-h))/h\n",
" \n",
"# Central difference for f'(x):\n",
"def diff_central(f, x, h=0.1):\n",
" return (f(x+h)-f(x-h))/(2*h)\n",
"# end of diff_central\n",
"\n",
"# Central difference for f''(x):\n",
"def diff2_central(f, x, h=0.1):\n",
" return (f(x+h)-2*f(x)+f(x-h))/h**2\n",
"# end of diff2_central"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Numerical example 1\n",
"x = pi/4;\n",
"df_exact = cos(x)\n",
"ddf_exact = -sin(x)\n",
"h = 0.1\n",
"f = sin\n",
"df = diff_forward(f, x, h)\n",
"print('Approximations to the first derivative')\n",
"print('Forward difference: df = {:12.8f}, Error = {:10.3e} '.format(df, df_exact-df))\n",
"df = diff_backward(f, x, h)\n",
"print('Backward difference: df = {:12.8f}, Error = {:10.3e} '.format(df, df_exact-df))\n",
"df = diff_central(f, x, h)\n",
"print('Central difference: df = {:12.8f}, Error = {:10.3e} '.format(df, df_exact-df))\n",
"print('Approximation to the second derivative') \n",
"ddf = diff2_central(f, x, h)\n",
"print('Central difference: ddf= {:12.8f}, Error = {:10.3e} '.format(ddf, ddf_exact-ddf))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Error analysis\n",
"In this case the error analysis is quite simple: Do a Taylor expansion of the\n",
"error around $x$. The Taylor expansion becomes a power series in $h$.\n",
"\n",
"\n",
"The expansion for the error of the forward difference is:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e(x;h) = f'(x) - \\frac{f(x+h)-f(x)}{h} = f'(x) - \\frac{f(x)+f'(x)h + \\frac{1}{2}f''(\\xi)h^2 - f(x)}{h} = -\\frac{1}{2}f''(\\xi)h\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\xi\\in (x,x+h)$. \n",
"\n",
"The expansion for the error of the central difference is slightly more complicated:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"e(x; h) &= f'(x) - \\frac{f(x+h)-f(x-h)}{2h} \\\\ \n",
" &= f'(x) \\\\ &- \\frac{\\big(f(x)+f'(x)h + \\frac{1}{2} f''(x)h^2 + \\frac{1}{6} f'''(\\xi_1)h^2 \\big) - \\big(f(x)-f'(x)h + \\frac{1}{2} f''(x)h^2 - \\frac{1}{6} f'''(\\xi_2)h^2\\big)}{2h} \\\\ \n",
" &= -\\frac{1}{12}\\big(f'''(\\xi_1) + f'''(\\xi_2)\\big)h^2 \\\\ \n",
" &= -\\frac{1}{6}f'''(\\eta)h^2, \\qquad \\qquad \\eta \\in (x-h, x+h),\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where the two remainder terms have been combined by the intermediate value theorem\n",
"(Result 2 at the end of *Preliminaries*). The error for the approximation of the\n",
"second order derivative can be found similarly. \n",
"\n",
"The order of an approximation is $p$ if there exist a constant $C$ independent on $h$ such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"|e(h;x) \\leq C h^p,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"see *Preliminaries*. \n",
"\n",
"In practice, it is sufficient to show that the power expansion of the error satisfies"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"e(x,h)=C_ph^{p}+ C_{p+1}h^{p+1} + \\dotsm, \\qquad C_p \\not=0\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The forward and backward approximations\n",
"are of order 1, the central differences of order 2. \n",
"\n",
"We are going to use these formulas a lot in the sequel, so let us just summarize\n",
"the results, including the error terms: \n",
"\n",
"\n",
"\n",
"**Difference formulas for derivatives:**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
" f'(x) &= \\left\\{\n",
" \\begin{array}{ll} \\displaystyle\n",
" \\frac{f(x+h)-f(x)}{h} - \\frac{h}{2}f''(\\xi), \\ & \\text{Forward difference} \\\\ \\mbox{} \\\\ \n",
" \\displaystyle \n",
" \\frac{f(x)-f(x-h)}{h} + \\frac{h}{2}f''(\\xi), & \\text{Backward difference} \\\\ \\mbox{} \\\\ \n",
" \\displaystyle\n",
" \\frac{f(x+h)-f(x-h)}{2h} - \\frac{h^2}{6}f'''(\\xi).\\qquad & \\text{Central difference} \n",
" \\end{array} \\right. \\\\ \\mbox{} \\\\ \n",
" f''(x) & = \n",
" \\frac{f(x+h)-2f(x)+f(x-h)}{h^2} - \\frac{h^2}{12}f^{(4)}(\\xi), \\qquad \n",
"\\text{Central difference}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Two point boundary problems (BVP)\n",
"\n",
"Given a two point boundary value problem:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u'' + p(x) u' + q(x) u = r(x), \\qquad a \\leq x \\leq b, \\qquad u(a)=u_a, \\quad u(b)=u_b,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $p$, $q$ are given functions of $x$ and the boundary values $u_a$ and $u_a$ are given constants.\n",
"\n",
"A finite difference method for this problem is constructed by the following steps: \n",
"\n",
"**Step 1:**\n",
"Given the interval $[a,b]$. Choose N, let $h=(b-a)/N$ and let $x_i=a+ih$, $i=0,1,\\dotsc,N$. \n",
"\n",
"**Step 2:**\n",
"For each inner grid point $x_i$, $i=1,\\dotsc,N-1$, replace the derivatives by their approximations in the BVP. The result is:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{u(x_i+h)-2u(x_i)+u(x_i-h)}{h^2} + p(x_i) \\frac{u(x_i+h)-u(x_i-h)}{2h} + q(x_i) u(x_i) + \\mathcal{O}(h^2) = r(x_i)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for each $i=1,2\\dotsc,N-1$, and the term $\\mathcal{O}(h^2)$ represents the errors in the difference formulas. \n",
"\n",
"**Step 3:** Ignore the error term, and replace the exact solution $u(x_i)$ by its numerical (and still unknown) approximation $U_i$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{U_{i+1}-2U_i+U_{i-1}}{h^2} + p(x_i)\\frac{U_{i+1}-U_{i-1}}{2h} + q(x_i) U_i = r(x_i), \\qquad i=1,\\dotsc,N-1.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This is the *discretization* of the BVP. If we know include the two boundary values as equations, the discretization is a linear system of equations"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"A \\mb{U} = \\mb{b},\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $A$ er en $N+1\\times N+1$ matrix and $\\mb{U} = [U_0,\\dotsc,U_{N}]^T$. Or more specific, by multiplying the equations by $h^2$ we end up with:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"A = \\left[ \\begin{array}{ccccccc}\n",
" 1 & 0 & & \\\\ \n",
" v_1 & d_1 & w_1 & & & \\\\ \n",
" & v_2 & d_2 & w_2 & \\\\ \n",
" & & v_3 & \\ddots & \\ddots & \\\\ \n",
" & & & \\ddots & \\ddots & w_{N-2} \\\\ \n",
" & & & & v_{N-1} & d_{N-1} & w_{N-1} \\\\ \n",
" & & & & & 0 & 1\n",
" \\end{array} \\right]\n",
" \\qquad \\text{with} \\qquad\n",
" \\begin{array}{l}\n",
" \\displaystyle v_i =1-\\frac{h}{2}p(x_i) \\\\ \n",
" \\displaystyle d_i = -2 + h^2q(x_i) \\\\ \n",
" \\displaystyle w_i = 1+\\frac{h}{2}p(x_i)\n",
" \\end{array}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The right hand side $\\mb{b}$ is given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\mb{b} = [u_a, h^2r(x_1), \\dotsc, h^2r(x_{N-1}), u_b]^T.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Obviously, the first and last equations are trivial to solve, and is therefore often included in the right hand side.\n",
"\n",
"**Step 4:** Solve $A \\mb{U} = \\mb{b}$ with respect to $\\mb{U}$.\n",
"\n",
"\n",
"**Example 1:**\n",
"Given the equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u'' + 2u' - 3u = 9x, \\qquad u(0)=u_a = 1, \\quad u(1)=u_b = e^{-3}+2e-5=0.486351,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with exact solution $u(x)= e^{-3x}+2e^{x}-3x-2$.\n",
"\n",
"Choose $N$, let $h=1/N$. Use the central differences for $u'$ and $u''$, such that"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{u(x_i+h)-2u(x_i)+u(x_i-h)}{h^2} + 2 \\frac{u(x_i+h)-u(x_i-h)}{2h} -3\n",
"u(x_i) + \\mathcal{O}(h^2) = 9 x_i, \\qquad i=1,\\dotsc, N\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let $U_i \\approx u(x_i)$. Multiply by $h^2$ on both sides, include $U_0=u_a$ og $U_N=u_b$ and clean the mess:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
" U_0 &= 1 \\\\ \n",
" (1-h)U_{i-1} + (-2-3h^2)U_i + (1+h)U_{i+1} &= 9x_ih^2, && i=1, \\cdots N-1, \\\\ \n",
" U_N &= 0.486351\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To be even more concrete, for $N=4$, we get $h=0.25$. The linear system of equations becomes"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left( \\begin{array}{ccccc}\n",
" 1 & 0 & 0 & 0 & 0\\\\ \n",
" 0.75 & -2.1875 & 1.25 & 0 & 0\\\\ \n",
" 0 & 0.75 & -2.1875 & 1.25 & 0\\\\ \n",
" 0 & 0 & 0.75 & -2.1875 & 1.25 \\\\ \n",
" 0 & 0 & 0 & 0 & 1\n",
"\\end{array} \\right) \n",
"\\left(\\begin{array}{c} U_0 \\\\ U_1 \\\\ U_2 \\\\ U_3 \\\\ U_4\n",
"\\end{array} \\right)\n",
"=\n",
"\\left( \\begin{array}{c} 1. \\\\ 0.140625 \\\\ 0.28125 \\\\ 0.421875 \\\\ 0.48635073\n",
"\\end{array} \\right).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The first and the last equation is trivial to solve, so in practice you have a system of 3 equations in 3 unknowns,"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left(\\begin{array}{ccc}\n",
" -2.1875 & 1.25 & 0 \\\\ \n",
" 0.75 & -2.1875 & 1.25 \\\\ \n",
" 0 & 0.75 & -2.1875\n",
"\\end{array} \\right)\n",
"\\left( \\begin{array}{c} U_1 \\\\ U_2 \\\\ U_3 \\end{array} \\right) =\n",
"\\left( \\begin{array}{c} 0.140625-0.75\\cdot 1 \\\\ 0.28125 \\\\ 0.421875-1.25 \\cdot\n",
"0.48635073 \\end{array} \\right),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with the solution"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"U_1 = 0.293176, \\qquad U_2= 0.025557, \\qquad 0.093820.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"For comparison, the exact solution in these points are:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(0.25) = 0.290417, \\qquad u(0.5) = 0.020573, \\qquad u(0.75) = 0.089400.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Implementation\n",
"For simplicity, the implementation below is only done for BVPs with constant\n",
"coefficients, that is $p(x)=p$ and $q(x)=q$. This makes the diagonal, sub- and\n",
"super-diagonals constant, except at the first and the last row. \n",
"An extra function is included to construct matrices of the form $A =\n",
"\\text{tridiag}\\{v,d,w\\}$. \n",
"\n",
"The implementation consist of\n",
"1. Choose $N$, let $h=(b-a)/N$ and $x_i=a+ih$, $i=0,\\dotsc,N$.\n",
"\n",
"2. Construct the matrix $A\\in \\mathbb{R}^{N+1\\times N+1}$ and the vector $b\\in\\mathbb{R}^{N+1}$. The matrix $A$ is tridiagonal, and except from the first and last row, has the elements $v=1-\\frac{h}{2}p$ below the diagonal, $d = -2 + h^2 q$ as diagonal elements and $w = 1+\\frac{h}{2}p$ above the diagonal. \n",
"\n",
"3. Construct the vector $\\mb{b} = [b_0,\\dotsc,b_N]^T$ with elements $b_i=h^2r(x_i)$ for $i=1,\\dotsc,N-1$.\n",
"\n",
"4. Modify the first and the last row of the matrix $A$ and the first and last element of the vector $\\mb{b}$, depending on the boundary conditions. \n",
"\n",
"5. Solve the system $A\\mb{U} = \\mb{b}$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def tridiag(v, d, w, N):\n",
" # Help function \n",
" # Returns a tridiagonal matrix A=tridiag(v, d, w) of dimension N x N.\n",
" e = ones(N) # array [1,1,...,1] of length N\n",
" A = v*diag(e[1:],-1)+d*diag(e)+w*diag(e[1:],1)\n",
" return A"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Example 1, BVP\n",
"\n",
"# Define the equation \n",
"# u'' + p*u' + q*u = r(x) on the interval [a,b]\n",
"# Boundary condition: u(a)=ua and u(b)=ub\n",
"\n",
"p = 2\n",
"q = -3\n",
"def r(x):\n",
" return 9*x\n",
"a, b = 0, 1\n",
"ua, ub = 1, exp(-3)+2*exp(1)-5\n",
"\n",
"# The exact solution (if known)\n",
"def u_eksakt(x):\n",
" return exp(-3*x)+2*exp(x)-3*x-2\n",
"\n",
"# Set up the discrete system\n",
"N = 4 # Number of intervals \n",
"\n",
"# Start the discretization \n",
"h = (b-a)/N # Stepsize\n",
"x = linspace(a, b, N+1) # The gridpoints x_0=a, x_1=a+h, .... , x_N=b \n",
"\n",
"# Make a draft of the A-matrix (first and last row have to be adjusted)\n",
"v = 1-0.5*h*p # Subdiagonal\n",
"d = -2+h**2*q # Diagonal\n",
"w = 1+0.5*h*p # Superdiagonal\n",
"A = tridiag(v, d, w, N+1) \n",
"\n",
"# Make a draft of the b-vector\n",
"b = h**2*r(x) \n",
"\n",
"# Modify the first equation (left boundary) \n",
"A[0,0] = 1\n",
"A[0,1] = 0\n",
"b[0] = ua\n",
" \n",
"# Modify the last equation (right boundary) \n",
"A[N,N] = 1 \n",
"A[N,N-1] = 0\n",
"b[N] = ub\n",
"\n",
"\n",
"U = solve(A, b) # Solve the equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"To verify the calculations done above, print the matrix $A$, the vector $\\mb{b}$ and the numerical solution $\\mb{U}$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Print the matrix A, the right hand side b the numerical and exact solution\n",
"print('A =\\n', A) \n",
"print('\\nb =\\n ', b)\n",
"print('\\nU =\\n ', U)\n",
"print('\\nu(x)=\\n', u_eksakt(x))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot the solution of the BVP\n",
"xe = linspace(0,1,101)\n",
"plot(x,U,'.-')\n",
"plot(xe, u_eksakt(xe),':') \n",
"xlabel('x')\n",
"ylabel('u')\n",
"legend(['Numerisk','Eksakt'])\n",
"title('Løsningen av et to-punkt randverdiproblem.');"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot the error |u(x)-U| in the gridpoints\n",
"error = abs(u_eksakt(x)-U)\n",
"plot(x, error,'.-')\n",
"xlabel('x')\n",
"ylabel('error')\n",
"title('Error: u(x)-U');\n",
"print('Max error = {:.3e}'.format(max(abs(error))))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Boundary conditions\n",
"\n",
"To get a unique solution of a BVP (or a PDE), some information about the\n",
"solution, usually given on the boundaries has to be known. The most common boundary conditions are:\n",
"1. Dirichlet condition: The solution is known at the boundary.\n",
"\n",
"2. Neumann condition: The derivative is known at the boundary.\n",
"\n",
"3. Robin (or mixed) condition: A combination of those.\n",
"\n",
"In the example above, Dirichlet conditions were used. We will now see\n",
"how to handle Neumann conditions. Robin conditions can be treated similarly. \n",
"\n",
"Given the BVP with a Neumann condition at the left boundary:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u'' + p(x) u' + q(x) u = r(x), \\qquad a \\leq x \\leq b, \\qquad u'(a)=u'_a,\n",
"\\quad u(b)=u_b.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Here, $u'_a$ is some given value. In this case, the solution $u(a)$ and its\n",
"corresponding approximation $U_0$ are unknown, and we need some difference\n",
"formula also for the point $a=x_0$. The simplest option is to use a forward\n",
"difference"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u'_a = \\frac{u(x_1)-u(x_0)}{h} + \\mathcal{O}(h) \\qquad \\Rightarrow \\qquad\n",
" \\frac{U_1-U_0}{h} = u'_a\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"but this is only a first order approximation, and thus lower accuracy is to be\n",
"expected. We could also use a second order approximation using the values in the\n",
"grid points $x_0$, $x_1$ and $x_2$, but this will ruin the nice tridiagonal\n",
"structure of the coefficient matrix. Instead, use the idea of a *false\n",
"boundary*: \n",
"\n",
"Assume that the solution can be stretched outside the boundary $x=a$, all the\n",
"way to a fictitious grid point $x_{-1}=a-h$, where we also assume there is an\n",
"approximate and equally fictitious approximation $U_{-1}$ to $u(x_{-1})$. Then\n",
"we have two difference formulas in the point $a$, one for the BVP itself\n",
"and a central difference for the boundary conditions:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
" \\frac{U_{1}-2U_0+U_{-1}}{h^2} + p(x_0)\\frac{U_{1}-U_{-1}}{2h} + q(x_0) U_0 & = r(x_0) \\\\ \n",
" \\frac{U_1 - U_{-1}}{2h} &= u'_a\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Solve the second equation with respect to $U_{-1}$, insert this into the first\n",
"equation which then becomes:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{2U_1-2U_0-2hu'_a}{h^2} + p(x_0)u'_a + q(x_0)U_0 = r(x_0).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So the only thing that has changed is the first equation. And since central differences have been used both for the BVP and the boundary condition, the overall order of this approximation can be proved to be 2. \n",
"\n",
"** Example 2:**\n",
"\n",
"Given the same example as before, but now with a Neumann condition at the left\n",
"boundary:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u'' + 2u' - 3u = 9x, \\qquad u'(0)=u'_a =-4, \\quad u(1)=u_b = -2e^{-3}+e-5 =\n",
"0.48635073,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"with exact solution $u(x)= e^{-3x}-2e^{x}-3x-2$.\n",
"\n",
"The modified difference equation at the boundary $x_0=0$ is:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{2U_1-2U_0-2u'_a h}{h^2} +2u'_a - 3U_0 = 0.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Multiply this equation by $h^2$, and include the equation as the \n",
"discretization for the grid point $x_0$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
" (-2-3h^2)U_0 - 2U_1 &= (2h-2h^2)u'_a \\\\ \n",
" (1-h)U_{i-1} + (-2-3h^2)U_i + (1+h)U_{i+1} &= 9 h^2 x_i, && i=1, \\cdots N-1. \\\\ \n",
" U_N &= u_b,\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"which, for $N=4$ og $h=0.25$ becomes:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\left( \\begin{array}{ccccc}\n",
" -2.1875 & 2 & 0 & 0 & 0 \\\\ \n",
" 0.75 & -2.1875 & 1.25 & 0 & 0 \\\\ \n",
" 0 & 0.75 & -2.1875 & 1.25 & 0 \\\\ \n",
" 0 & 0 & 0.75 & -2.1875 & 1.25 \\\\ \n",
" 0 & 0 & 0 & 0 & 1 \\\\ \n",
"\\end{array}\\right)\\left(\\begin{array}{c} U_{0} \\\\ U_1 \\\\ U_2 \\\\ U_3 \\\\ U_4\n",
"\\end{array}\\right) = \\left(\\begin{array}{c} -1.5 \\\\ 0.140625 \\\\ 0.28125 \\\\ \n",
"0.421875 \\\\ 0.48635073 \\end{array} \\right).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution of this is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"U_0 = 0.92103219, \\quad U_1 = 0.25737896, \\quad U_2 = 0.01029386, \\qquad U_3 = 0.08858688.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Numerical exercises:**\n",
"1. Modify the code above to solve this problem. Use $N=4$ to check your solution, but try also $N=10$ and $N=20$. \n",
"\n",
"2. Modify the code above to solve the same BVP, but now with the left boundary condition \\\\ $u'(a)+u(a)/4=0$. \n",
"\n",
"# The heat equation\n",
"\n",
"In this section we will see how to solve the heat equation by finite difference\n",
"methods. It should however be emphasized that the strategies applies to a lot of \n",
"time-dependent PDEs, the heat equation is just an example. \n",
"\n",
"Given the equation, well known from the first part of this course:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
" u_t & = u_{xx}, && 0 \\leq x \\leq 1 \\\\ \n",
" u(0,t) &=g_0(t), \\quad u(1,t)=g_1(t), && \\text{Boundary conditions} \\\\ \n",
" u(x,0) &= f(x) && \\text{Initial conditions}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The equation is solved from $t=0$ to $t=t_{end}$. \n",
"\n",
"## Semi-discretization\n",
"This is a technique which combines the discretization of boundary problems\n",
"explained above with the techniques for solving ordinary differential equations. \n",
"\n",
"The idea is as follows:\n",
"\n",
"**Step 1:**\n",
"Make a discretization of the interval in the $x$-direction. Choose some $M$, let\n",
"$\\Delta x=1/M$ (since the interval is $[0,1]$) and the grid points are $x_i=i\\Delta\n",
"x$, $i=0,1,\\dotsc,M$. \n",
"\n",
"Note that for each grid point $x_i$ the solution $u(x_i,t)$ is a function of $t$\n",
"alone. \n",
"\n",
"**Step 2:**\n",
"Fix some arbitrary point time $t$, and make a discretization of the right hand side of the\n",
"PDE. Using central differences to approximate $u_{xx}$, this will give"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\frac{\\partial u}{\\partial t}(x_i,t) =\n",
"\\frac{u(x_{i+1},t)-2u(x_i,t)+u(x_{i-1},t)}{\\Delta x^2} + \\mathcal{O}(\\Delta\n",
"x^2).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Step 3:** \n",
"Ignore the error term $\\mathcal{O}(\\Delta x^2)$ and replace $u(x_i,t)$ with the\n",
"approximation $U_i(t)$ in the formula above. The result is"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"U'_i(t) = \\frac{U_{i+1}(t) - 2 U_i(t) + U_{i-1}(t)}{\\Delta x^2}, \\qquad i=1,2,\\dotsc,M-1,\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $U'_i(t) = dU_i(t)/dt$. \n",
"And this, together with the boundary conditions $U_0(t)=g_0(t)$, $U_M(t)=g_1(t)$\n",
"and the initial condition $U_i(0) = f(x_i)$, $i=0,1,\\dotsc,M$ forms a well defined\n",
"system of ordinary differential equations. \n",
"\n",
"This system is usually called a *semi-discretization* of the PDE. \n",
"\n",
"**Step 4:**\n",
"Solve the system of ODEs by the method of your preference. \n",
"\n",
"The explicit Euler method with step size $\\Delta t$ applied to these ODEs is:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"U_{i}^{n+1} = U_{i}^n + r \\big(\\; U_{i+1}^n - 2U_{i}^n\n",
"+ U_{i-1}^n\\; \\big), \\quad i=1,2,\\dotsc,M-1, \\qquad \\text{where } r = \\frac{\\Delta t}{\\Delta x^2}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So $U_i^n \\approx u(x_i,t_n)$ where $t_n=n\\Delta t$, and to distinguish between\n",
"the indices representing space and time, we have used the time indices $n$ as\n",
"superscripts. \n",
"\n",
"Let us test this algorithm on two examples.\n",
"\n",
"**Numerical examples 1:** \n",
"Solve the heat equation $u_t=u_{xx}$ on the interval $[0,1]$ with"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\begin{align*}\n",
"u(x,0) & = \\sin(\\pi x), && \\text{Initial value} \\\\ \n",
"g_0(t) & =g_1(t)=0. && \\text{Boundary values}\n",
"\\end{align*}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The exact solution of this problem is given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x,t) = e^{-\\pi^2t}\\sin(\\pi x).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Numerical example 2:**\n",
"Repeat example 1, but now with the initial values"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u(x,0) = \\left\\{ \\begin{array}{ll} 2x & 0 \\leq x \\leq 0.5, \\\\ 2(1-x) & 0.5 < x\\leq 2-2x, \n",
"\\end{array} \\right.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"In this case, the exact solution is not given. \n",
"\n",
"Run the codes below with \n",
"1. $M=4$, $N=20$.\n",
"\n",
"2. $M=8$, $N=40$.\n",
"\n",
"3. $M=16$, $N=80$. \n",
"\n",
"Both initial values are already implemented. \n",
"\n",
"## Implementation\n",
"\n",
"We first include a function for plotting the solution."
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"def plot_heat_solution(x, t, U, txt='Solution'):\n",
" # Help function\n",
" # Plot the solution of the heat equation\n",
" fig = figure()\n",
" ax = fig.gca(projection='3d')\n",
" T, X = meshgrid(t,x)\n",
" # ax.plot_wireframe(T, X, U)\n",
" ax.plot_surface(T, X, U, cmap=cm.coolwarm)\n",
" ax.view_init(azim=30) # Rotate the figure\n",
" xlabel('t')\n",
" ylabel('x')\n",
" title(txt);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Define the problem, this time in terms of initial values and boundary conditions."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
" # Define the problem\n",
"\n",
" # Initial condition\n",
" def f1(x): # Example 1\n",
" return sin(pi*x)\n",
"\n",
" def f2(x): # Example 2\n",
" y = 2*x\n",
" y[x>0.5] = 2-2*x[x>0.5]\n",
" return y\n",
"\n",
" f = f1\n",
" \n",
" # Boundary conditions\n",
" def g0(t):\n",
" return 0\n",
" def g1(t):\n",
" return 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The main part of the code is:"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Solve the heat equation by a forward difference in time (forward Euler)\n",
"#\n",
"M = 4 # Number of intervals in the x-direction\n",
"Dx = 1/M\n",
"x = linspace(0,1,M+1) # Gridpoints in the x-direction\n",
"\n",
"tend = 0.5\n",
"N = 20 # Number of intervals in the t-direction\n",
"Dt = tend/N\n",
"t = linspace(0,tend,N+1) # Gridpoints in the t-direction\n",
"\n",
"# Array to store the solution\n",
"U = zeros((M+1,N+1))\n",
"U[:,0] = f(x) # Initial condition U_{i,0} = f(x_i)\n",
"\n",
"r = Dt/Dx**2 \n",
"print('r =',r)\n",
"\n",
"# Main loop \n",
"for n in range(N):\n",
" U[1:-1, n+1] = U[1:-1,n] + r*(U[2:,n]-2*U[1:-1,n]+U[0:-2,n]) \n",
" U[0, n+1] = g0(t[n+1])\n",
" U[M, n+1] = g1(t[n+1])"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot the numerical solution\n",
"plot_heat_solution(x, t, U)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Plot the error from example 1\n",
"def u_exact(x,t):\n",
" return exp(-pi**2*t)*sin(pi*x)\n",
"T, X = meshgrid(t, x)\n",
"error = u_exact(X, T) - U\n",
"plot_heat_solution(x, t, error, txt='Error')\n",
"print('Maximum error: {:.3e}'.format(max(abs(error.flatten())))) # Maksimal feil over hele arrayet"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The solution is stable for $M=4$, $N=20$, and apparently unstable for $M=16$, $N=80$. \n",
"Why? \n",
"\n",
"## Stability analysis\n",
"\n",
"The semi-discretized system"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\dot{U}_i(t) = \\frac{U_{i+1}(t) - 2 U_i(t) + U_{i-1}(t)}{\\Delta x^2}, \\qquad i=1,2,\\dotsc,M-1, \\qquad U_0(t)=g_0(t), \\qquad U_M(t)=g_1(t),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"is a linear ordinary differential equation:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\dot{\\mb{U}} = \\frac{1}{\\Delta x^2}\\big( A \\mb{U} + \\mb{g}(t)\\big),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\mb{U} = \\left(\\begin{array}{c} U_1 \\\\ U_2 \\\\ \\vdots \\\\ U_{M-1} \\end{array} \\right), \\qquad\n",
"A = \\left(\\begin{array}{cccc}\n",
" -2 & 1 & \\\\ \n",
" 1 & \\ddots & \\ddots \\\\ \n",
" & \\ddots & \\ddots & 1 \\\\ \n",
" & & 1 & -2\n",
" \\end{array}\\right)\n",
"\\qquad \\text{and} \\qquad\n",
"\\mb{g}(t) = \\left(\\begin{array}{c} g_0(t) \\\\ 0 \\\\ \\vdots \\\\ 0 \\\\ g_1(t) \\end{array} \\right)\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Stability requirements for such problems were discussed in the note on stiff ordinary differential equation. We proved there that the stability depends on the eigenvalues $\\lambda_k$ of the matrix $\\frac{1}{\\Delta x^2}A$. For the forward Euler method, it was proved that the step size has to chosen such that $-2 \\leq \\Delta t \\lambda_k \\leq 0$ for all $\\lambda_k$. Otherwise, the numerical solution will be unstable.\n",
"\n",
"It is possible to prove that the eigenvalues of the matrix $A$ is given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\lambda_k = -4\\sin^2 \\big( \\frac{k\\pi}{M}\\big), \\qquad k=1,\\cdots,M-1.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So all the eigenvalues $\\lambda_k$ of $\\frac{1}{\\Delta x^2}A$ satisfy"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"-\\frac{4}{\\Delta x^2} < \\lambda _k < 0.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The numerical solution is stable if $\\Delta t < -2/\\lambda_k$ for all $k$, that means whenever"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"r = \\frac{\\Delta t}{\\Delta x^2} \\leq \\frac{1}{2}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Exercise:**\n",
"Repeat the two experiments above (for the two different initial values) to justify the bound above. \n",
"Use $M=16$, and in each case find the corresponding $r$ and observe from the experiments whether the solution is stable or not. \n",
"\n",
"1. Let $N=256$.\n",
"\n",
"2. Let $N=128$. \n",
"\n",
"3. Let $N=250$. \n",
"\n",
"In the last case, it seems like the metod is unstable for the first initial value, and unstable for the second. Do you have any idea why? (Both solutions will be unstable if integrated over a longer time periode). \n",
"\n",
"**Hint:** Relate to the Fourier expansion solution of the heat equation from the first part of the course. \n",
"\n",
"## Implicit methods\n",
"\n",
"The semi-discretized system is an example of a stiff ODE, which can only be handled reasonable efficiently by $A(0)$-stable methods, like implicit Euler or the trapezoidal rule, see the note on stiff ODEs. \n",
"\n",
"\n",
"### Implicit Euler\n",
"\n",
"The implicit Euler method for the discritized system \n",
"$\\dot{\\mb{U}}=\\frac{1}{\\Delta x^2}\\;\\big(A\\mb{U} + \\mb{g}(t)\\big) $ \n",
"is given by"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\mb{U}^{n+1} = \\mb{U}^{n} + r \\,A \\mb{U}^{n+1} + r \\,\\mb{g}(t_{n+1}), \\qquad \\text{med} \\qquad r=\\frac{\\Delta t}{\\Delta x^2}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $\\mb{U}^n=[U_1^n, U_2^n, \\dotsc, U_{M-1}^n$ and $U_i^n \\approx u(x_i, t_n)$.\n",
"\n",
"For each time step, the following system of linear equations has to be solved:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(I_{M-1}- r \\,A)\\mb{U}^{n+1} = \\mb{U}^n + r\\, \\mb{g}(t_{n+1}),\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where $I_{M-1}$ is the identity matrix of dimension $(M-1)\\times (M-1)$. \n",
"\n",
"The error in the gridpoints can be shown to be $\\mathcal{O}(\\Delta t+ \\Delta x^2)$.\n",
"\n",
"### Crank-Nicolson (trapezoidal rule)\n",
"\n",
"The trapezoidal rule applied to the semi-discretized system is often referred to as *Crank-Nicolson's method*. The method is $A(0)$-stable and of order 2 in time, so we can expect better accuracy.\n",
"The method is written as:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"\\mb{U}^{n+1} = \\mb{U}^n + \\frac{\\Delta t}{2\\Delta x^2}A\\big(\\mb{U}^{n+1}+\\mb{U}^n\\big) + \\frac{\\Delta t}{2\\Delta x^2}\\big(\\mb{g}(t_n)+\\mb{g}(t_{n+1})\\big).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"So for each timestep the following system of equations has to be solved with respect to $\\mb{U}^n$:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"(I_{M-1} - \\frac{r}{2}A)\\mb{U}^{n+1} = (I_{M-1} + \\frac{r}{2}A)\\mb{U}^{n} + \\frac{r}{2}\\big(\\mb{g}(t_n)+\\mb{g}(t_{n+1}\\big), \\qquad r=\\frac{\\Delta t}{\\Delta x^2}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The error in the gridpoints can be shown to be $\\mathcal{O}(\\Delta t^2 + \\Delta x^2)$.\n",
"\n",
"\n",
"### Implementation\n",
"\n",
"It is possible to solve the system of ODEs directly by the methods developed in the note on stiff ODEs, or by using some other existing ODE solver. For nonlinear problems, this is often advisable (but not always). Mostly for the purpose of demonstration, the implicit Euler method as well as the Crank-Nicolson scheme is implemented directly in the following code. \n",
"\n",
"For each time step, a system of linear equation has to be solved:"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"K \\mb{U}_{n+1} = \\mb{b}\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"where\n",
"\n",
"**Implicit Euler:**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"K = I_{M-1} - rA, \\qquad \\mb{b} = \\mb{U}_n + r [g_0(t_{n+1}), 0, \\dotsc, 0, g_1(t_{n+1})]^T\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Crank-Nicolson:**"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"K = I_{M-1} - \\frac{r}{2}A, \\qquad \\mb{b} = (I_{M-1}+\\frac{r}{2}A)\\mb{U}_n + r \\bigg[\\frac{1}{2}(g_0(t_n)+g_0(t_{n+1})), 0, \\dotsc, 0,\\frac{1}{2}(g_1(t_n)+g_1(t_{n+1}))\\bigg]^T.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The methods can of course be applied to the problems from Numerical examples 1 and 2. But for the fun of it, we now include a problem with nontrivial boundary conditions. \n",
"\n",
"**Numerical example 3:**\n",
"Solve the equation"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"$$\n",
"u_t = u_{xx}, \\qquad u(0,t) = e^{-\\pi^2 t}, \\quad u(1,t) = -e^{-\\pi^2 t}, \\quad u(x,0) = \\cos(\\pi x).\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"up to $t_{end}=0.2$ by implicit Euler and Crank-Nicolson. Plot the solution and the error. \n",
"The exact solution is $u(x,t) = e^{-\\pi^2t}\\cos(\\pi x)$. \n",
"\n",
"Use $N=M$, and $M=10$ and $M=100$ (for example). \n",
"Notice that there are no stability issues, even for $r$ large. \n",
"Also notice the difference in accuracy for the two methods."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Apply implicit Euler and Crank-Nicolson on\n",
"# the heat equation u+t=u_{xx}\n",
"\n",
"# Define the problem of example 3\n",
"def f3(x):\n",
" return cos(pi*x)\n",
"\n",
"# Boundary values\n",
"def g0(t):\n",
" return exp(-pi**2*t)\n",
"def g1(t):\n",
" return -exp(-pi**2*t)\n",
"\n",
"# Exact solution \n",
"def u_exact(x,t):\n",
" return exp(-pi**2*t)*cos(pi*x)\n",
"\n",
"f = f3\n",
"\n",
"# Choose method\n",
"method = 'iEuler'\n",
"# method = 'CrankNicolson'\n",
"\n",
"\n",
"M = 100 # Number of intervals in the x-direction\n",
"Dx = 1/M\n",
"x = linspace(0,1,M+1) # Gridpoints in the x-direction\n",
"\n",
"tend = 0.5\n",
"N = M # Number of intervals in the t-direction\n",
"Dt = tend/N\n",
"t = linspace(0,tend,N+1) # Gridpoints in the t-direction\n",
"\n",
"# Array to store the solution\n",
"U = zeros((M+1,N+1))\n",
"U[:,0] = f(x) # Initial condition U_{i,0} = f(x_i)\n",
"\n",
"# Set up the matrix K: \n",
"A = tridiag(1, -2, 1, M-1)\n",
"r = Dt/Dx**2\n",
"print('r = ', r)\n",
"if method is 'iEuler':\n",
" K = eye(M-1) - r*A\n",
"elif method is 'CrankNicolson':\n",
" K = eye(M-1) - 0.5*r*A\n",
"\n",
"Utmp = U[1:-1,0] # Temporary solution for the inner gridpoints. \n",
"\n",
"# Main loop over the time steps. \n",
"for n in range(N):\n",
" # Set up the right hand side of the equation KU=b:\n",
" if method is 'iEuler':\n",
" b = copy(Utmp) # NB! Copy the array\n",
" b[0] = b[0] + r*g0(t[n+1])\n",
" b[-1] = b[-1] + r*g1(t[n+1]) \n",
" elif method is 'CrankNicolson':\n",
" b = dot(eye(M-1)+0.5*r*A, Utmp)\n",
" b[0] = b[0] + 0.5*r*(g0(t[n])+g0(t[n+1]))\n",
" b[-1] = b[-1] + 0.5*r*(g1(t[n])+g1(t[n+1])) \n",
" \n",
" Utmp = solve(K,b) # Solve the equation K*Utmp = b\n",
" \n",
" U[1:-1,n+1] = Utmp # Store the solution\n",
" U[0, n+1] = g0(t[n+1]) # Include the boundaries. \n",
" U[M, n+1] = g1(t[n+1])"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"plot_heat_solution(x, t, U)"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"T, X = meshgrid(t, x)\n",
"error = u_exact(X, T) - U\n",
"plot_heat_solution(x, t, error, txt='Error')\n",
"print('Maximum error: {:.3e}'.format(max(abs(error.flatten()))))"
]
}
],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}