CV asked me to cover some basic linear PDE solvers and two-point ODE boundary value problems, specifically:
I don't know exactly what he's covered thus far, but am assuming that you've seen (and mostly remember) some basic ODE initial value problem solvers and something about finite difference approximations of derivatives. I'm in the office next door to CV, and often we both have our doors open during his OH, so I think that's a fair assumption.
I'm also going to use this as an excuse to play with a teaching technique that I've contemplated for a while, namely real-time demonstrations with a notebook environment.
Let's start with a model two-point BVP: $$ -\frac{d^2 u}{dx^2} = f \mbox{ on } (0,1), \quad u(0) = u(1) = 0 $$ It's always good to start with problems that we know how to solve by hand, so let's use a simple right-hand side, $f \equiv 1$. Solutions to the differential equation (without boundary conditions) take the form $$ u(x) = -\frac{1}{2}x^2 + bx + c $$ and the boundary conditions take the form $$ \begin{bmatrix} u(0) \\ u(1) \end{bmatrix} = \begin{bmatrix} 0 \\ -0.5 \end{bmatrix} + \begin{bmatrix} 0 & 1 \\ 1 & 1 \end{bmatrix} \begin{bmatrix} b \\ c \end{bmatrix} = \begin{bmatrix} 0 \\ 0 \end{bmatrix}. $$ We can solve this by hand: $c = 0$ and $b = 0.5$.
And let's draw a pretty picture, just to show that we can.
using PyPlot
x = linspace(0,1)
uref(x) = -x.^2/2 + x/2
plot(x, uref(x))
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x11e399690>
Let's pretend that we never learned how to solve this differential equation in closed form. An alternative is to use an ODE solver. Most of the time, we need to put the equation in first-order form, so let's start there: $$ \frac{d}{dx} \begin{bmatrix} u \\ v \end{bmatrix} = \begin{bmatrix} v \\ -f \end{bmatrix}, \quad u(0) = u(1) = 0 $$ The problem is that the ODE solver requires $u(0)$ and $v(0)$ rather than $u(0)$ and $u(1)$. Let's initially pretend that we know $v(0)$ and work from there.
using Sundials
v0 = 0.0
function ode_rhs(t, uv, uvdot)
uvdot[1] = uv[2]
uvdot[2] = -1.0
end
res = Sundials.cvode(ode_rhs, [0.0; v0], x)
plot(x, uref(x))
plot(x, res[:,1], "--")
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x11961b150>
OK, this is clearly not right. But the only thing not right about is it $v(0)$!
So let's set up an iteration to find the right $v(0)$, using cvode
to
evaluate the mapping $v(0) \mapsto u(1)$. I'm going to do a simple
bisection solve, and augment the mapping function so that we can watch
the convergence of the iteration.
plot(x, uref(x))
function ufinal(v0)
res = Sundials.cvode(ode_rhs, [0.0; v0], x)
plot(x, res[:,1], "--")
return res[end,1]
end
function bisection_solver(f, xlo, xhi, tol=1e-3)
flo = f(xlo)
fhi = f(xhi)
print("Initial bracket: [$xlo, $xhi]\n")
while abs(xlo-xhi) > tol
xmid = (xlo+xhi)/2
fmid = f(xmid)
if flo*fmid > 0.0
xlo, flo = xmid, fmid
elseif flo*fmid < 0.0
xhi, fhi = xmid, fmid
end
print(" refine: [$xlo, $xhi]\n")
end
(xlo, xhi)
end
bisection_solver(ufinal, 0.0, 2.0)
Initial bracket: [0.0, 2.0] refine: [0.0, 1.0] refine: [0.5, 1.0] refine: [0.5, 0.75] refine: [0.5, 0.625] refine: [0.5, 0.5625] refine: [0.5, 0.53125] refine: [0.5, 0.515625] refine: [0.5, 0.5078125] refine: [0.5, 0.50390625] refine: [0.5, 0.501953125] refine: [0.5, 0.5009765625]
(0.5,0.5009765625)
It's possible to write down what's called a variational equation to compute derivatives of the mapping from the conditions at one end to conditions at the other; this is particularly helpful if the shooting problem leads to a system of equations rather than a scalar equation. Since the topic of this class is mostly not linear or nonlinear equation solver technology, I'll just make a plug here for CS 4220 next semester, where we'll talk about such things in detail.
The shooting approach involves solving a set of equations involving values of the solution at the boundary. An alternate approach is to form a set of equations involving values of the solution at the boundary and on the interior.
Let's set up a mesh of $N+2$ points (including boundary points) with spacing $h = 1/(N+1)$ on $[0,1]$: $$ 0 = x_0 < x_1 < x_2 < x_3 < \ldots < x_N < x_{N+1} = 1 $$ When we need to illustrate, we'll use $N = 4$.
You should remember that you can approximate $u_i'' \approx u''(x_i)$ from $u_i \approx u(x_i)$ by the finite difference formula $$ u_i'' = h^{-2} \left( u_{i-1}-2u_i+u_{i+1} \right). $$ Putting thing together, we can think about the (approximate) mapping from function values to (negative) second derivative at all the interior mesh points simultaneously: $$ -\begin{bmatrix} u_1'' \\ u_2'' \\ u_3'' \\ u_4'' \end{bmatrix} = h^{-2} \begin{bmatrix} -1 & 2 & -1 & 0 & 0 & 0 \\ 0 & -1 & 2 & -1 & 0 & 0 \\ 0 & 0 & -1 & 2 & -1 & 0 \\ 0 & 0 & 0 & -1 & 2 & -1 \end{bmatrix} \begin{bmatrix} u_0 \\ u_1 \\ u_2 \\ u_3 \\ u_4 \\ u_5 \end{bmatrix} $$
To solve the BVP with general boundary values and RHS, then, we would look at something like $$ h^{-2} \left( \begin{bmatrix} 2 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 2 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix} + \begin{bmatrix} -1 \\ 0 \\ 0 \\ 0 \end{bmatrix} u_0 + \begin{bmatrix} 0 \\ 0 \\ 0 \\ -1 \end{bmatrix} u_5 \right) = \begin{bmatrix} f_1 \\ f_2 \\ f_3 \\ f_4 \end{bmatrix} $$ Let's rearrange a bit: $$ \begin{bmatrix} 2 & -1 & 0 & 0 \\ -1 & 2 & -1 & 0 \\ 0 & -1 & 2 & -1 \\ 0 & 0 & -1 & 2 \end{bmatrix} \begin{bmatrix} u_1 \\ u_2 \\ u_3 \\ u_4 \end{bmatrix} = h^2 \begin{bmatrix} f_1 \\ f_2 \\ f_3 \\ f_4 \end{bmatrix} + \begin{bmatrix} u_0 \\ 0 \\ 0 \\ u_5 \end{bmatrix} $$ For the case of homogeneous Dirichlet conditions ($u_1 = u_5 = 0$), we have $Tu = h^2 f$, where $T$ is the standard tridiagonal matrix.
OK. Now let's try this out in practice. First, let's make sure we can construct the standard tridiagonal. We won't worry about being dense.
NN = length(x)
N = NN-2
Tdense(N) = 2*eye(N) - diagm(ones(N-1),1) - diagm(ones(N-1),-1)
T = Tdense(N)
T[1:5,1:5]
5x5 Array{Float64,2}: 2.0 -1.0 0.0 0.0 0.0 -1.0 2.0 -1.0 0.0 0.0 0.0 -1.0 2.0 -1.0 0.0 0.0 0.0 -1.0 2.0 -1.0 0.0 0.0 0.0 -1.0 2.0
OK! Now let's do a solve and plot.
h = 1.0/(N+1)
ufd = zeros(NN)
f = ones(NN)
ufd[2:N+1] = T\(h^2*f[2:N+1])
plot(x,uref(x))
plot(x,ufd, "--")
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x11e764850>
Just for kicks, let's do something a little more efficient and construct a sparse representation of this matrix.
Tsparse(N) = 2*speye(N) - spdiagm(ones(N-1),1,N,N) - spdiagm(ones(N-1),-1,N,N)
T = Tsparse(N)
T[1:5,1:5]
5x5 sparse matrix with 13 Float64 entries: [1, 1] = 2.0 [2, 1] = -1.0 [1, 2] = -1.0 [2, 2] = 2.0 [3, 2] = -1.0 [2, 3] = -1.0 [3, 3] = 2.0 [4, 3] = -1.0 [3, 4] = -1.0 [4, 4] = 2.0 [5, 4] = -1.0 [4, 5] = -1.0 [5, 5] = 2.0
And we're going to do the same thing with the sparse version.
ufd[2:N+1] = T\(h^2*f[2:N+1])
plot(x, uref(x))
plot(x, ufd, "--")
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x11a825a50>
This is great -- we have a good match in the "eyeball norm." Let's actually do a plot of the error.
plot(x, uref(x)-ufd)
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x119383190>
This is the type of error plot that should make you suspicious. It's too good! This is because the second-order finite difference formula is exact when applied to a quadratic. Let's take a moment to set up a more challenging (sort of) problem.
A nice function that is zero at $0$ and $1$ is $$ u(x) = \sin(\pi x) $$ and we have $$ -u'' = f = \pi^2 \sin(\pi x). $$ So let's set this up as a BVP and look at the error for different size meshes.
uref(x) = sin(pi*x)
fref(x) = pi^2 * sin(pi*x)
mmax = 8
hs = zeros(mmax)
errnorm = zeros(mmax)
for m = 3:mmax
h = 2.0^-m
N = 2^m-1
x = linspace(0,1,N+2)
ufd = Tsparse(N)\(h^2*fref(x[2:end-1]))
hs[m] = h
errnorm[m] = norm(uref(x[2:end-1])-ufd, Inf)
end
loglog(hs[3:mmax], errnorm[3:mmax])
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x119bd1090>
This looks pretty linear. Let's look at the step-by-step behavior. We should be cutting the error by roughly a factor of four each time we cut the mesh step by a factor of two ($O(h^2)$).
errnorm[4:mmax]./errnorm[3:mmax-1]
5-element Array{Float64,1}: 0.248554 0.249639 0.24991 0.249977 0.249994
So far, we've looked at a steady-state model problem (an elliptic differential equation). Let's move on now to one of the simplest time-dependent problems: the heat equation, a model parabolic differential equation: $$ \frac{\partial u}{\partial t} = \kappa \frac{\partial^2 u}{\partial x^2} + f , \quad u(0,t) = u(1,t) = 0 $$ with $u(x,0)$ given. This equation is common not only in physical models of diffusion, but also in financial models (e.g. the Black-Scholes equation is a parabolic differential equation).
What's the simplest thing that we might think about doing? Perhaps the simplest approach is to discretize in space in order to get a system of ODEs, then solve the system using a standard time stepper. This is sometimes called the method of lines. For example, if $u^h(t)$ represents the vector of (approximate) $u(x_i,t)$ values for some mesh with spacing $h$, then we have the equations $$ \frac{du^h}{dt} = \kappa h^{-2} T u^h(t) + f^h(t) $$ where $f^h$ represents the sampled vector of forcing function values.
Let's see what happens if we try to solve this differential equation using Euler's method in time.
N = 20
h = 1.0/(N+1)
x = linspace(0,1,N+2)
T = Tsparse(N)
u = zeros(N)
f = ones(N)
kappa = 1
dt = 1e-2
usteady(x) = (-x.^2+x)/2
plot(x[2:end-1], u)
for j = 1:10
u += (-(kappa/h^2)*(T*u) + f) * dt
plot(x[2:end-1], u)
end
(D,V) = eig(Tdense(100))
lambdas = D
print(lambdas[1], " ", lambdas[end], "\n")
0.0009674354160247115 3.999032564583976
This massively depends on the time step! If the time step is too large, the whole thing goes unstable. This should not surprise us, as the spectrum of $T$ has eigenvalues from nearly one down to $O(h^2)$. That is, the discretized system is very stiff, and the stiffness gets much worse as the mesh density increases. In order to maintain stability with Euler's method, we would need a time step of $\Delta t = O(h^2/\kappa)$. An alternative approach, of course, is to run backward Euler. This time, as a diagnostic, let's plot the solution at each time step and also the steady-state solution (which we computed earlier as the solution of an elliptic problem).
u = zeros(N)
plot(x[2:end-1], u)
A = speye(N) + dt*(kappa/h^2)*T
for j = 1:50
u = A\(u+dt*f)
plot(x[2:end-1], u)
end
plot(x, usteady(x), "k--")
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x11b87cc90>
u = zeros(N)
A0 = speye(N) + dt*(kappa/h^2/2)*T
A1 = speye(N) - dt*(kappa/h^2/2)*T
for j = 1:50
u = A0\(A1*u+dt*f)
plot(x[2:end-1], u)
end
plot(x, usteady(x), "k--")
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x11e2ad090>
Finally, let's consider the wave equation, a simple hyperbolic PDE that models wave propagation. The PDE is $$ \frac{\partial^2 u}{\partial t^2} = c^2 \frac{\partial^2 u}{\partial x^2}, \quad u(0) = u(1) = 0 $$ with given initial conditions. We could add a forcing term, but we're not going to for this lecture.
As in the parabolic case, we first discretize in space, then in time. Discretizing in space gives $$ \frac{d^2 u^h(t)}{dt} = -(c/h)^2 T u $$ Now we apply a second-order finite difference in time as well: $$ u^h(t_{k-1})-2u^h(t_k)+u^h(t_{k+1}) = -(c \, \Delta t/h)^2 T u. $$ Rearranging this a bit gives us a simple solver iteration. $$ u^h(t_{k+1}) = 2u^h(t_k)-u^h(t_{k-1})-\xi^2 T u, \quad \xi = c \frac{\Delta t}{h} $$ The numerical wave speed $\xi$ is a measure of how fast waves move in the numerical method (in units of grid points per time step). We need to make sure that $\xi$ is sufficiently small (less than one) in order to remain stable; this is sometimes called a CFL (Courant-Friedrichs-Levy) condition.
Let's draw some pretty pictures of the wave equation applied to a Gaussian bump. The main point of this exercise is to make sure that we actually understand what we're doing well enough to make progress. So far, I've written all the other solvers in this document; assuming we get this far, let's do this one collaboratively!
N = 100
x = linspace(0,1,N+2)
u0 = exp(-100*(x-0.5).^2)
u0 -= u0[1]
plot(x,u0)
1-element Array{Any,1}: PyObject <matplotlib.lines.Line2D object at 0x119943410>
I was a bit silly about this in class -- I got the numerical method right, but forgot to make copies of the data in the right places. If I do set up the copies appropriately, we get a very pretty picture of the propagating waves.
xi = 0.4
u0 = exp(-100*(x-0.5).^2)
u0 -= u0[1]
u1 = copy(u0)
unext = copy(u1)
T = Tsparse(N)
for j = 1:100
unext[2:end-1] = 2*u1[2:end-1]-u0[2:end-1]-xi^2*(T*u1[2:end-1])
if mod(j,10) == 0
plot(x, unext)
end
u0, u1, unext = u1, unext, u0
end