2022-04-20 Stiff equations
Contents
2022-04-20 Stiff equations¶
Last time¶
Notes on integration
Ordinary differential equations (ODE)
Explicit methods for solving ODE
Stability
Today¶
Implicit and explicit methods
Stiff equations
PDE examples
using LinearAlgebra
using Plots
default(linewidth=4, legendfontsize=12)
struct RKTable
A::Matrix
b::Vector
c::Vector
function RKTable(A, b)
s = length(b)
A = reshape(A, s, s)
c = vec(sum(A, dims=2))
new(A, b, c)
end
end
function rk_stability(z, rk)
s = length(rk.b)
1 + z * rk.b' * ((I - z*rk.A) \ ones(s))
end
rk4 = RKTable([0 0 0 0; .5 0 0 0; 0 .5 0 0; 0 0 1 0], [1, 2, 2, 1] / 6)
heun = RKTable([0 0; 1 0], [.5, .5])
Rz_theta(z, theta) = (1 + (1 - theta)*z) / (1 - theta*z)
function ode_rk_explicit(f, u0; tfinal=1, h=0.1, table=rk4)
u = copy(u0)
t = 0.
n, s = length(u), length(table.c)
fY = zeros(n, s)
thist = [t]
uhist = [u0]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
for i in 1:s
ti = t + h * table.c[i]
Yi = u + h * sum(fY[:,1:i-1] * table.A[i,1:i-1], dims=2)
fY[:,i] = f(ti, Yi)
end
u += h * fY * table.b
t = tnext
push!(thist, t)
push!(uhist, u)
end
thist, hcat(uhist...)
end
ode_rk_explicit (generic function with 1 method)
Ordinary Differential Equations¶
Given initial condition \(y_0 = y(t=0)\), find \(y(t)\) for \(t > 0\) that satisfies
Application |
\(y\) |
\(f\) |
---|---|---|
Orbital dynamics |
position, momentum |
conservation of momentum |
Chemical reactions |
concentration |
conservation of atoms |
Epidemiology |
infected/recovered population |
transmission and recovery |
\(y\) can be a scalar or a vector
Solving differential equations¶
Linear equations¶
Autonomous if \(A(t) = A\) and source independent of \(t\)
Suppose \(y\) and \(a = A\) are scalars: \(y(t) = e^{at} y_0\)
Can do the same for systems¶
What does it mean to exponentiate a matrix?¶
Taylor series!
and there are many practical ways to compute it.
Question¶
Suppose that the diagonalization \(A = X \Lambda X^{-1}\) exists and derive a finite expression for the matrix exponential using the scalar exp
function.
Forward Euler Method¶
The simplest method for solving \(y'(t) = f(t,y)\) is to use numerical differentiation to write
function ode_euler(f, y0; tfinal=10., h=0.1)
y = copy(y0)
t = 0.
thist = [t]
yhist = [y0]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
y += h * f(t, y)
t = tnext
push!(thist, t)
push!(yhist, y)
end
thist, hcat(yhist...)
end
ode_euler (generic function with 1 method)
f1(t, y; k=5) = -k * (y .- cos(t))
thist, yhist = ode_euler(f1, [1.], tfinal=10, h=.2)
scatter(thist, yhist[1,:])
plot!(cos)
Forward Euler on a linear system¶
f2(t, y) = [0 1; -1 0] * y
thist, yhist = ode_euler(f2, [0., 1], h=.5, tfinal=10)
scatter(thist, yhist')
plot!([cos, sin])
eigvals([0 1; -1 0])
2-element Vector{ComplexF64}:
0.0 - 1.0im
0.0 + 1.0im
Runge-Kutta 4¶
thist, yhist = ode_rk_explicit(f2, [0., 1], h=.5, tfinal=50)
scatter(thist, yhist')
plot!([cos, sin], size=(800, 500))
Apparently it is possible to integrate this system using large time steps.
This method evaluates \(f(y)\) four times per stepso the cost is about equal when the step size \(h\) is 4x larger than forward Euler.
Linear Stability Analysis¶
Why did Euler diverge (even if slowly) while RK4 solved this problem accurately? And why do both methods diverge if the step size is too large? We can understand the convergence of methods by analyzing the test problem
function plot_stability(Rz, method; xlim=(-3, 2), ylim=(-1.5, 1.5))
x = xlim[1]:.02:xlim[2]
y = ylim[1]:.02:ylim[2]
plot(title="Stability: $method", aspect_ratio=:equal, xlim=xlim, ylim=ylim)
heatmap!(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2))
contour!(x, y, (x, y) -> abs(Rz(x + 1im*y)), color=:black, linewidth=2, levels=[1.])
plot!(x->0, color=:black, linewidth=1, label=:none)
plot!([0, 0], [ylim...], color=:black, linewidth=1, label=:none)
end
plot_stability (generic function with 1 method)
plot_stability(z -> 1 + z, "Forward Eulor")
Stability for RK4¶
plot_stability(z -> rk_stability(4z, rk4), "RK4")
plot_stability(z -> rk_stability(2z, heun), "Heun's method")
Implicit methods¶
Recall that forward Euler is the step
plot_stability(z -> 1/(1-z), "Backward Euler")
Computing with implicit methods¶
Linear solve for linear problem
Nonlinear (often Newton) solve for nonlinear problem
Need Jacobian or finite differencing
plot_stability(z -> Rz_theta(z, .5), "Midpoint method")
The \(\theta\) method¶
Forward and backward Euler are bookends of the family known as \(\theta\) methods.
which, for linear problems, is solved as
\(\theta=0\) is explicit Euler, \(\theta=1\) is implicit Euler, and \(\theta=1/2\) are the midpoint or trapezoid rules (equivalent for linear problems). The stability function is
Rz_theta(z, theta) = (1 + (1-theta)*z) / (1 - theta*z)
theta = .6
plot_stability(z -> Rz_theta(z, theta), "\$\\theta=$theta\$",
xlim=(-20, 1))
\(\theta\) method for the oscillator¶
function ode_theta_linear(A, u0; forcing=zero, tfinal=1, h=0.1, theta=.5)
u = copy(u0)
t = 0.
thist = [t]
uhist = [u0]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
rhs = (I + h*(1-theta)*A) * u .+ h*forcing(t+h*theta)
u = (I - h*theta*A) \ rhs
t = tnext
push!(thist, t)
push!(uhist, u)
end
thist, hcat(uhist...)
end
ode_theta_linear (generic function with 1 method)
# Test on oscillator
A = [0 1; -1 0]
theta = 1
thist, uhist = ode_theta_linear(A, [0., 1], h=.2, theta=theta, tfinal=20)
scatter(thist, uhist')
plot!([cos, sin])
\(\theta\) method for the cosine decay¶
k = 50
theta = .5
thist, uhist = ode_theta_linear(-k, [.2], forcing=t -> k*cos(t), tfinal=5, h=.5, theta=theta)
scatter(thist, uhist[1,:], title="\$\\theta = $theta\$")
plot!(cos, size=(800, 500))
Stability classes and the \(\theta\) method¶
Definition: \(A\)-stability¶
A method is \(A\)-stable if the stability region
Definition: \(L\)-stability¶
A time integrator with stability function \(R(z)\) is \(L\)-stable if
Heat equation as linear ODE¶
How do different \(\theta \in [0, 1]\) compare in terms of stability?
Are there artifacts even when the solution is stable?
using SparseArrays
function heat_matrix(n)
dx = 2 / n
rows = [1]
cols = [1]
vals = [0.]
wrap(j) = (j + n - 1) % n + 1
for i in 1:n
append!(rows, [i, i, i])
append!(cols, wrap.([i-1, i, i+1]))
append!(vals, [1, -2, 1] ./ dx^2)
end
sparse(rows, cols, vals)
end
heat_matrix(5)
5×5 SparseMatrixCSC{Float64, Int64} with 15 stored entries:
-12.5 6.25 ⋅ ⋅ 6.25
6.25 -12.5 6.25 ⋅ ⋅
⋅ 6.25 -12.5 6.25 ⋅
⋅ ⋅ 6.25 -12.5 6.25
6.25 ⋅ ⋅ 6.25 -12.5
n = 200
A = heat_matrix(n)
x = LinRange(-1, 1, n+1)[1:end-1]
u0 = exp.(-200 * x .^ 2)
@time thist, uhist = ode_theta_linear(A, u0, h=.1, theta=1, tfinal=1);
nsteps = size(uhist, 2)
plot(x, uhist[:, 1:5])
1.398294 seconds (4.36 M allocations: 229.282 MiB, 5.36% gc time, 99.92% compilation time)
Advection as a linear ODE¶
function advect_matrix(n; upwind=false)
dx = 2 / n
rows = [1]
cols = [1]
vals = [0.]
wrap(j) = (j + n - 1) % n + 1
for i in 1:n
append!(rows, [i, i])
if upwind
append!(cols, wrap.([i-1, i]))
append!(vals, [1., -1] ./ dx)
else
append!(cols, wrap.([i-1, i+1]))
append!(vals, [1., -1] ./ 2dx)
end
end
sparse(rows, cols, vals)
end
advect_matrix(5)
5×5 SparseMatrixCSC{Float64, Int64} with 11 stored entries:
0.0 -1.25 ⋅ ⋅ 1.25
1.25 ⋅ -1.25 ⋅ ⋅
⋅ 1.25 ⋅ -1.25 ⋅
⋅ ⋅ 1.25 ⋅ -1.25
-1.25 ⋅ ⋅ 1.25 ⋅
n = 50
A = advect_matrix(n, upwind=false)
x = LinRange(-1, 1, n+1)[1:end-1]
u0 = exp.(-9 * x .^ 2)
@time thist, uhist = ode_theta_linear(A, u0, h=.04, theta=1, tfinal=1.);
nsteps = size(uhist, 2)
plot(x, uhist[:, 1:(nsteps÷8):end])
0.127159 seconds (179.07 k allocations: 11.513 MiB, 13.49% gc time, 83.18% compilation time)
Spectrum of operators¶
theta=.5
h = .1
plot_stability(z -> Rz_theta(z, theta), "\$\\theta=$theta, h=$h\$")
ev = eigvals(Matrix(h*advect_matrix(20, upwind=true)))
scatter!(real(ev), imag(ev))
theta=.5
h = .2 / 4
plot_stability(z -> Rz_theta(z, theta), "\$\\theta=$theta, h=$h\$")
ev = eigvals(Matrix(h*heat_matrix(20)))
scatter!(real(ev), imag(ev))