Finite Differences¶

The Nature of Mathematical Modeling (draft)¶

$ \def\%#1{\mathbf{#1}} \def\mat#1{\mathbf{#1}} \def\*#1{\vec{#1}} \def\ve#1{\vec{#1}} \def\ds#1#2{\cfrac{d #1}{d #2}} \def\dd#1#2{\cfrac{d^2 #1}{d #2^2}} \def\ps#1#2{\cfrac{\partial #1}{\partial #2}} \def\pp#1#2{\cfrac{\p^2 #1}{\p #2^2}} \def\p{\partial} \def\ba{\begin{eqnarray*}} \def\ea{\end{eqnarray*}} \def\eps{\epsilon} \def\del{\nabla} \def\disp{\displaystyle} $

ODEs¶

Taylor series expansion¶

Solve general first-order differential equation

$$\ds{y}{x} = f(x,y)$$

Generalizes to system of equations, and therefore to higher order:

$$ \ba \frac{dy_1}{dx} &=& f_1(x,y_1,\ldots,y_N) \\ \frac{dy_2}{dx} &=& f_2(x,y_1,\ldots,y_N) \\ &\vdots& \\ \frac{dy_N}{dx} &=& f_N(x,y_1,\ldots,y_N) \ea $$

Do a Taylor expansion of the solution, and substitute in the differential equation:

$$ \ba y(x+h) &=& y(x) + \left. h\frac{dy}{dx} \right|_x + \left. \frac{h^2}{2}\frac{d^2y}{dx^2} \right|_x + {\cal O}(h^3) \\ &=& y(x) + h f(x,y(x)) + \frac{h^2}{2} \frac{d}{dx} f(x,y(x)) + {\cal O}(h^3) \\ &=& y(x) + h f(x,y(x)) + \frac{h^2}{2} \left[\frac{\p f}{\p x} + \frac{\p f}{\p y} \frac{\p y}{\p x}\right] + {\cal O}(h^3) \\ &=& y(x) + h f(x,y(x)) + \frac{h^2}{2} \left[\frac{\p f}{\p x} + f \frac{\p f}{\p y}\right] + {\cal O}(h^3) \ea $$

Solvers must match terms up to the desired order

Euler's method¶

First two terms are Euler's method

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
scale = .6
width = 10*scale
height = 5*scale
marker = 2*height
linew = marker/2.5
font = 5.5*height
fig,ax = plt.subplots(figsize=(width,height))
ax.axis([0,width,0,height])
ax.set_axis_off()
def line(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k-',linewidth=0.8*linew)
def point(x,y):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=marker)
def text(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        fontsize=font,color='black')
def arrow(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    ax.annotate('',xy=(x1,y1),xytext=(x0,y0),
        arrowprops=dict(color='black',width=marker/7,headwidth=marker,headlength=marker))
arrow(1,1,1,4)
arrow(1,1,8,1)
line(2.5,.8,2.5,1.2)
line(6.5,.8,6.5,1.2)
point(2.5,2.5)
point(6.5,2.01)
point(6.5,3.5)
text(2.5,0.5,'$x$')
text(6.5,0.5,'$x+h$')
text(1,4.35,'$y$')
arrow(2.45,2.47,6.35,3.45)
def curve(x0,x1,a,b,c):
    x0 = scale*x0
    x1 = scale*x1
    x = np.linspace(x0,x1,100)
    y = a+b*x+c*x*x
    plt.plot(x,y,'k-',linewidth=0.8*linew)
curve(2,7,0.8,0.69,-0.15),
text(2.5,3,'$y(x)$')
text(7.4,2.4,'$y(x+h)$')
text(6.5,3.9,'$y(x)+hf(x,y(x))$')
plt.show()

Example: $dy/dx = Ay$, solved by $y(x) = \exp(Ax)$
Euler approximation: $y(x+h) = (1+hA)y(x)$
Euler approximation can itself be solved: $y(x) = (1+hA)^{x/h}y(0)$
- If $A > 0$ diverges, as expected
- If $0 > hA > -1$ solution decays to zero, as expected
- If $-1 > hA > -2$ solution incorrectly oscillates
- If $-2 > hA$ solution incorrectly diverges
Issue is jumping between solution curves

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
scale = .6
width = 12*scale
height = 8*scale
marker = 1.5*height
linew = marker/4
font = 3.5*height
fig,ax = plt.subplots(figsize=(width,height))
ax.axis([0,width,0,height])
ax.set_axis_off()
def line(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k-',linewidth=0.8*linew)
def point(x,y):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=marker)
def text(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        fontsize=font,color='black')
def arrow(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    ax.annotate('',xy=(x1,y1),xytext=(x0,y0),
        arrowprops=dict(color='black',width=marker/7,headwidth=marker,headlength=marker))
x0 = 1*scale
x1 = 12*scale
h = 0.1
A = -0.5
def curve(y0):
    x = np.linspace(x0,x1,100)
    y = scale*(4+(1+h*A)**(x/h)*y0)
    plt.plot(x,y,'-',color='0.75',linewidth=linew)
for y0 in range(-20,21,1):
    curve(y0)
def curvepoint(x,y0):
    x = x*scale
    y = scale*(4+(1+h*A)**(x/h)*y0)
    plt.plot(x,y,'ko',markersize=marker)
#for x in range(1,12,1):
#    curvepoint(x,4)
def step(h,N):
    x = x0
    h = h*scale
    y = 3
    for i in range(N):
        plt.plot(x,scale*(4+y),'ko',markersize=marker/1.2)
        xnew = x+h
        ynew = (1+h*A)*y
        plt.plot([x,xnew],[scale*(4+y),scale*(4+ynew)],'k--')
        x = xnew
        y = ynew
        plt.plot(x,scale*(4+y),'ko',markersize=marker/1.2)
step(1,10)
step(3,3)
step(5.4,2)
arrow(1,4,1,8)
arrow(1,4,1,0)
arrow(1,4,12,4)
text(12.25,4,'$x$')
text(.65,7,'$y$')
plt.show()

Midpoint method¶

Rather than evaluating the slope at the start of the interval, an Euler step can be used to evaluate the slope in the middle: $$ y(x+h) = y(x) + h f\left[x+\frac{h}{2},~y(x)+\frac{h}{2}f(x,y(x))\right] $$
This is called the midpoint method, and improves the approximation to second order (Problem 1)

Runge-Kutta method¶

Family of expansions, most common is fourth-order:

$$ \ba k_1 &=& hf(x,y(x)) \\ k_2 &=& hf\left(x+\frac{h}{2},~y(x)+\frac{k_1}{2}\right) \\ k_3 &=& hf\left(x+\frac{h}{2},~y(x)+\frac{k_2}{2}\right) \\ k_4 &=& hf(x+h,~y(x)+k_3) \\ y(x+h) &=& y(x) + \frac{k_1}{6} + \frac{k_2}{3} + \frac{k_3}{3} + \frac{k_4}{6} + {\cal O}(h^5) \ea $$

For a system of equations:

$$ \ba k_{1,i} &=& hf_i(x,y_1,\ldots,y_N) \\ k_{2,i} &=& hf_i\left(x+\frac{h}{2},~y_1+\frac{k_{1,1}}{2}, \ldots, y_N+\frac{k_{1,N}}{2}\right) \\ k_{3,i} &=& hf_i\left(x+\frac{h}{2},~y_1+\frac{k_{2,1}}{2}, \ldots, y_N+\frac{k_{2,N}}{2}\right) \\ k_{4,i} &=& hf_i(x+h,~y_1+k_{3,1}, \ldots, y_N+k_{3,N}) \\ y_i(x+h) &=& y_i(x) + \frac{k_{1,i}}{6} + \frac{k_{2,i}}{3} + \frac{k_{3,i}}{3} + \frac{k_{4,i}}{6} + {\cal O}(h^5) \ea $$

Problem 2 shows the impact of the order on error scaling.

Error estimation¶

Step size error in a fourth-order method:

$$ y(x+h) - y_{\rm true}(x+h) = h^5\phi(x) + {\cal O}(h^6) $$

where $\phi$ which is approximately constant over the interval (to order $h^5$)

Half-step error:

$$ y(x+h/2) - y_{\rm true}(x+h/2) = \left(\frac{h}{2}\right)^5\phi(x) + {\cal O}(h^6) $$

Difference between the full step and two half steps gives an estimate of the local error:

$$ \ba y(x+h) - y(x+h/2+h/2) &=& h^5\phi(x) - 2\left(\frac{h}{2}\right)^5\phi(x) + {\cal O}(h^6) \\ &\approx& h^5\phi \ea $$

The local error estimate can be used to adjust the step size up and down to match a target error
It's also possible to eliminate $\phi$ and combine the estimates and improve the order:

$$ y_{\rm true}(x+h) = y(x+h/2+h/2) + \frac{y(x+h/2+h/2)-y(x+h)}{15} + {\cal O}(h^6) $$

PDEs¶

Because PDEs have more degrees of freedom, there are more alternatives for discretization
Each type of PDE introduces new issues and solutions

Hyperbolic equations¶

example: flux¶

$$ \ps{u}{t} = -v \ps{u}{x} $$

Solved by any function of the form $u=f(x-vt)$
Index solution as $u(j\Delta x,n\Delta t) = u_j^n$
Take finite differences symmetrically:

$$ \frac{u_j^{n+1}-u_j^n}{\Delta t} = -v\left(\frac{u_{j+1}^n-u_{j-1}^n}{2\Delta x}\right) $$$$ u_j^{n+1} = u_j^n - \frac{v\Delta t}{2\Delta x}(u_{j+1}^n-u_{j-1}^n) $$

Visualize algorithm as a computational cluster

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
scale = .6
width = 10*scale
height = 3*scale
marker = 15*scale
linew = marker/2.5
font = 25*scale
fig,ax = plt.subplots(figsize=(width,height))
ax.axis([0,width,0,height])
ax.set_axis_off()
def line(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k-',linewidth=0.8*linew)
def point(x,y):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=marker)
def text(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        fontsize=font,color='black')
text(0.2,0.2,'$u_{j-1}^n$')
text(1.6,0.2,'$u_j^n$')
text(3.0,0.2,'$u_{j+1}^n$')
point(0.2,0.8)
point(1.6,0.8)
point(3.0,0.8)
line(0.2,0.8,1.6,0.8)
line(1.6,0.8,3.0,0.8)
line(1.6,0.8,1.6,2.2)
point(0.2,2.2)
point(1.6,2.2)
point(3.0,2.2)
text(0.2,2.7,'$u_{j-1}^{n+1}$')
text(1.6,2.7,'$u_j^{n+1}$')
text(3.0,2.8,'$u_{j+1}^{n+1}$')
plt.show()

von Neumann stability analysis¶

Locally linearize, separate variables, look at growth of ansatz:

$$ u_j^n = A(k)^n e^{ijk\Delta x} $$

For the flux equation:

$$ \ba A^{n+1}e^{ijk\Delta x} &=& A^ne^{ijk\Delta x} - \frac{v\Delta t}{2\Delta x} \left(A^ne^{i(j+1)k\Delta x}-A^ne^{i(j-1)k\Delta x}\right) \\ A &=& 1 - \frac{v\Delta t}{2\Delta x} \left(e^{ik\Delta x}-e^{-ik\Delta x}\right) \\ &=& 1 - i\frac{v\Delta t}{\Delta x}\sin k\Delta x \ea $$

The absolute magnitude of this is always greater than 1, and so this scheme is always unstable. Any initial condition will diverge!

Courant-Friedrichs-Levy stability criterion¶

For the flux equation, average neighbors in the time derivation:

$$ u_j^{n+1} = \frac{1}{2}(u_{j+1}^n + u_{j-1}^n) - \frac{v\Delta t}{2\Delta x}(u_{j+1}^n-u_{j-1}^n) $$

Repeat the stability analysis:

$$ A = \cos k\Delta x - i \frac{v\Delta t}{\Delta x}\sin k\Delta x $$

Require that the magnitude be less than 1:

$$ |A|^2 = \cos^2k\Delta x + \left(\frac{v\Delta t}{\Delta x}\right)^2 \sin^2k\Delta x \le 1 $$$$ \Rightarrow \frac{|v|\Delta t}{\Delta x} \le 1 $$

This has a natural interpretation: numerical propagation must be faster than physical propagation

numerical dissipation¶

Rewrite by subtracting $u_j^n$ from both sides:

$$ \frac{u_j^{n+1}-u_j^n}{\Delta t} = -v \left( \frac{u_{j+1}^n-u_{j-1}^n}{2\Delta x} \right) +\frac{1}{2\Delta t}(u_{j+1}^n-2u_j^n+u_{j-1}^n) $$

This corresponds to our original flux equation with an extra fictitious diffusion term added that depends on the discretization:

$$ \ps{u}{t} = -v\ps{u}{x} + \frac{(\Delta x)^2}{2\Delta t} \pp{u}{x} $$

These kinds of terms can be intentionally added for stability; Problem 3 shows the impact.

Parabolic equations¶

example: diffusion¶

$$ \ps{u}{t} = D\pp{u}{x} $$

Straightforward discretization:

$$ u_j^{n+1} = u_j^n + \frac{D\Delta t}{(\Delta x)^2} \left[u_{j+1}^n-2u_j^n+u_{j-1}^n\right] $$

Stability analysis:

$$ \ba A &=& 1 + \frac{D\Delta t}{(\Delta x)^2} \underbrace{\underbrace{ \left[e^{ik\Delta x} - 2 + e^{-ik\Delta x}\right]}_{\displaystyle 2\cos k\Delta x - 2}}_{\displaystyle 2\left(2\cos^2\frac{k\Delta x}{2} - 1\right) - 2} \\ &=& 1 - \frac{4D\Delta t}{(\Delta x)^2} \sin^2 \frac{k\Delta x}{2} \\ |A| \le 1 &\Rightarrow& \frac{4D\Delta t}{(\Delta x)^2} \le 2 ~\Rightarrow~ \frac{2D\Delta t}{(\Delta x)^2} \le 1 \ea $$

The method is stable for small step sizes, but since for a diffusive process the time $t$ to expand a distance $L$ is roughly $t\sim L^2/D$, the number of time steps required to model this will be $\sim L^2/(\Delta x)^2$ (i.e., a very large number).

implicit methods¶

Improve by evaluating the space derivative forwards in time:

$$ \ba &\frac{\disp u_j^{n+1} - u_j^n}{\disp \Delta t} = D \left[ \frac{\disp u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}} {\disp (\Delta x)^2} \right] & \\ &u_j^{n+1} - \frac{\disp D\Delta t}{\disp (\Delta x)^2} \left[u_{j+1}^{n+1}-2u_j^{n+1}+u_{j-1}^{n+1}\right] = u_j^n \ea $$

This is now unconditionally stable:

$$ \ba A - \frac{D\Delta t}{(\Delta x)^2} \left[Ae^{ik\Delta x}-2A+Ae^{-ik\Delta x}\right] &=& 1 \\ A\left[1+\frac{4D\Delta t}{(\Delta x)^2}\sin^2\frac{k\Delta x}{2}\right] &=& 1 \\ A &=& \frac{1}{\disp 1+\frac{4D\Delta t}{(\Delta x)^2} \sin^2\frac{k\Delta x}{2}} \le 1 \ea $$

This scheme is stable for all step sizes, but might appear to be useless: how can we implement it since we don't know the forward values used in the space derivative? These future values are implicitly determined by the past values, and the trick is to recognize that the full set of equations can be inverted; the stability then follows because peeking into the future in this way helps move information through the solution more quickly. This can be written as a matrix problem:

$$ \left(\begin{array}{c c c c c c c} 1 & 0 & 0 & \cdots & & & 0 \\ -\alpha & 1+2\alpha & -\alpha & 0 & \cdots & & 0 \\ & \ddots & \ddots & \ddots & \\ \vdots & 0 & -\alpha & 1+2\alpha & -\alpha & 0 & \vdots \\ & & & \ddots & \ddots & \ddots & \\ 0 & & \cdots & 0 & -\alpha & 1+2\alpha & -\alpha \\ 0 & & & \cdots & 0 & 0 & 1 \\ \end{array}\right) \left(\begin{array}{c} u_1^{n+1}\\ u_2^{n+1}\\ \vdots\\ u_i^{n+1}\\ \vdots\\ u_{N-1}^{n+1}\\ u_N^{n+1} \end{array}\right) = \left(\begin{array}{c} u_1^n \\ u_2^n \\ \vdots \\ u_i^n \\ \vdots \\ u_{N-1}^n \\ u_N^n \end{array}\right) $$

This is a tridiagonal matrix that is easily inverted with a linear-time pass up down and up the matrix.

In higher dimensions the matrix becomes banded matrix. While these can still be inverted more efficiently, it can be possible to improve stability by alternating implicit and explicit directions with the Alternating-Direction Implicit (ADI) method.

Elliptic equations¶

example: Poisson's¶

$$\del^2 u = \rho$$

Appears in electric fields, fluid flow, heat transfer, ...
Can write as the steady-state solution of a diffusion equation:

$$ \ps{u}{t} = \del^2 u - \rho $$

Discretize:

$$ u_{j,k}^{n+1} = u_{j,k}^n + \frac{\Delta t}{(\Delta x)^2} (u_{j+1,k}^n+u_{j-1,k}^n+u_{j,k+1}^n+u_{j,k-1}^n - 4u_{j,k}^n) - \Delta t \rho_{j,k} $$

relaxation methods¶

In 2D the Courant condition is $\Delta t / (\Delta x)^2 \le 1/4$, leading to:

$$ u_{j,k}^{n+1} = \frac{1}{4} (u_{j+1,k}^n + u_{j-1,k}^n + u_{j,k+1}^n + u_{j,k-1}^n) - \frac{(\Delta x)^2}{4} \rho_{j,k} $$

This has a very natural interpretation: starting from a random guess, at each time step each lattice site is set to the average of its neighbors and then a source term is added. This process is repeated until the solution stops changing, a relaxation technique. Successive Over-Relaxation (SOR) accelerates convergence by over-shooting the relaxation steps.

Problems¶

Show that the midpoint method is correct to second order.

Implement the Euler and fourth-order Runge-Kutta differential equation methods, use them to solve $\ddot x+x=0$ for an initial-value problem, and compare the dependence on the step size of the average error from the analytical solution.

Numerically solve for the motion of a pendulum that is parametrically driven with periodic forcing $l\ddot\theta + (g + \ddot z)\sin \theta = 0$ and explore its dynamics.

import matplotlib.pyplot as plt
import matplotlib.patches as patches
import numpy as np
scale = .6
width = 3*scale
height = 3*scale
line = width/1
marker = 5*width
font = 8*width
fig,ax = plt.subplots(figsize=(width,height))
ax.axis([0,width,0,height])
ax.set_axis_off()
def resistor(x,y,width):
    x = scale*x
    y = scale*y
    width = scale*width
    dx = width/16
    dy = 2.5*dx
    xplot = np.linspace(x-width/2,x+width/2,17)
    yplot = y+np.array([0,0,0,dy,0,-dy,0,dy,0,-dy,0,dy,0,-dy,0,0,0])
    plt.plot(xplot,yplot,'k-',linewidth=line)
def sline(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k-',linewidth=line)
def dline(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    plt.plot([x0,x1],[y0,y1],'k--',linewidth=line)
def point(x,y):
    x = scale*x
    y = scale*y
    plt.plot(x,y,'ko',markersize=marker)
def text(x,y,text):
    x = scale*x
    y = scale*y
    ax.text(x,y,text,
        ha='center',va='center',
        math_fontfamily='cm',
        fontsize=font,color='black')
def arrow(x0,y0,x1,y1):
    x0 = scale*x0
    x1 = scale*x1
    y0 = scale*y0
    y1 = scale*y1
    ax.annotate('',xy=(x1,y1),xytext=(x0,y0),
        arrowprops=dict(color='black',width=marker/7,headwidth=marker,headlength=marker))
def rectangle(x,y,w,h):
    x = scale*x
    y = scale*y
    w = scale*w
    h = scale*h
    #rect = patches.Rectangle((x,y),w,h,facecolor='gray')
    rect = patches.Rectangle((x,y),w,h,facecolor=(0.7,0.7,0.7))
    ax.add_patch(rect)
rectangle(0.75,2.25,1.5,.5)
dline(1.5,2.25,1.5,0.1)
sline(1.5,2.25,2.5,0.25)
point(2.5,0.25)
text(2.9,0.25,'$m$')
text(2.3,1.25,'$l$')
text(1.7,1.25,'$\\theta$')
arrow(2.75,2.25,2.75,3)
text(2.75,2.05,'$z$')
plt.show()

Simulate a plucked string by numerically solving the wave equation $\pp{u}{t} = v^2\pp{u}{x}+\gamma\ps{}{t}\pp{u}{x}$ with fixed boundary conditions. Calculate the stability of the undpamped equation ($\gamma=0$), and explore its dynamics with and without dissipation. Sonify the string motion.