# 2021-09-27 Transient problems¶

## Last time¶

recap/discuss methods for differentiating

intro to ordinary differential equations

accuracy and stability

## Today¶

A and L stability

transient advection

transient heat equation

```
using Plots
using LinearAlgebra
using SparseArrays
using Zygote
function vander(x, k=nothing)
if k === nothing
k = length(x)
end
V = ones(length(x), k)
for j = 2:k
V[:, j] = V[:, j-1] .* x
end
V
end
function fdstencil(source, target, k)
"kth derivative stencil from source to target"
x = source .- target
V = vander(x)
rhs = zero(x)'
rhs[k+1] = factorial(k)
rhs / V
end
function newton(residual, jacobian, u0; maxits=20)
u = u0
uhist = [copy(u)]
normhist = []
for k in 1:maxits
f = residual(u)
push!(normhist, norm(f))
J = jacobian(u)
delta_u = - J \ f
u += delta_u
push!(uhist, copy(u))
end
uhist, normhist
end
```

```
newton (generic function with 1 method)
```

# Ordinary Differential Equations¶

Given initial condition \(u_0 = u(t=0)\), find \(u(t)\) for \(t > 0\) that satisfies

Application |
\(u\) |
\(f\) |
---|---|---|

Orbital dynamics |
position, momentum |
conservation of momentum |

Chemical reactions |
concentration |
conservation of atoms |

Epidemiology |
infected/recovered population |
transmission and recovery |

Heat transfer |
temperature |
conservation of energy |

Seismology |
displacement, momentum |
conservative of momentum |

# Solving differential equations¶

## Linear equations¶

Autonomous if \(A(t) = A\) and source independent of \(t\)

Suppose \(u\) and \(a = A\) are scalars: \(u(t) = e^{at} u_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¶

```
function ode_euler(f, u0; tfinal=10., h=0.1)
u = copy(u0)
t = 0.
thist = [t]
uhist = [u0]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
u += h * f(t, u)
t = tnext
push!(thist, t)
push!(uhist, u)
end
thist, hcat(uhist...)
end
```

```
ode_euler (generic function with 1 method)
```

```
f1(t, u; k=20) = -k * (u .- cos(t))
thist, uhist = ode_euler(f1, [.5], tfinal=5, h=.102)
scatter(thist, uhist[1,:])
plot!(cos)
```

# Forward Euler on a linear system¶

```
f2(t, u) = [0 1; -1 0] * u
thist, uhist = ode_euler(f2, [0., 1], h=.1, tfinal=20)
scatter(thist, uhist')
plot!([cos, sin])
```

```
eigen([0 1; -1 0])
```

```
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
2-element Vector{ComplexF64}:
0.0 - 1.0im
0.0 + 1.0im
vectors:
2×2 Matrix{ComplexF64}:
0.707107-0.0im 0.707107+0.0im
0.0-0.707107im 0.0+0.707107im
```

# Linear Stability Analysis¶

Why did Euler converge on the first problem (with sufficiently small steps) while always diverging on the second problem?

We can understand the convergence of methods by analyzing the test problem

When does this map cause solutions to “blow up” and when is it stable?

# Stability diagrams (a property of methods)¶

```
function plot_stability(Rz, title; xlims=(-2, 2), ylims=(-2, 2))
x = LinRange(xlims[1], xlims[2], 100)
y = LinRange(ylims[1], ylims[2], 100)
heatmap(x, y, (x, y) -> abs(Rz(x + 1im*y)), c=:bwr, clims=(0, 2), aspect_ratio=:equal, title=title)
end
plot_stability(z -> 1 + z, "Forward Euler")
```

# Implicit methods¶

Recall that forward Euler is the step

**explicitly**; all the terms on the right hand side are known so the approximation \(\tilde u(h)\) is computed merely by evaluating the right hand side. Let’s consider an alternative,

**backward Euler**(or “implicit Euler”),

```
plot_stability(z -> 1 / (1 - z), "Backward Euler")
```

# Advection with centered stencil¶

```
function advect_center(t, u)
n = length(u)
dx = 2 / n
uwrap = [u[end], u..., u[1]]
f = zero(u)
for i in 1:n
f[i] = -(uwrap[i+2] - uwrap[i]) / (2*dx)
end
f
end
x = LinRange(-1, 1, 51)[1:end-1]
u0 = exp.(-9 * x .^ 2)
thist, uhist = ode_euler(advect_center, u0, h=.02, tfinal=1);
```

```
nsteps = size(uhist, 2)
@show nsteps
plot(x, uhist[:, 1:(nsteps÷8):end])
```

```
nsteps = 51
```

# Advection with upwind stencil¶

```
function advect_upwind(t, u)
n = length(u)
dx = 2 / n
uwrap = [u[end], u..., u[1]]
f = zero(u)
for i in 1:n
f[i] = -(uwrap[i+1] - uwrap[i]) / dx
end
f
end
x = LinRange(-1, 1, 51)[1:end-1]
u0 = exp.(-9 * x .^ 2)
@time thist, uhist = ode_euler(advect_upwind, u0, h=.002, tfinal=1);
```

```
0.850757 seconds (206.39 k allocations: 11.905 MiB, 99.75% compilation time)
```

```
nsteps = size(uhist, 2)
plot(x, uhist[:, 1:(nsteps÷8):end])
```

# Centered advection with Runge-Kutta¶

```
function ode_rk4(f, u0; tfinal=10., h=0.4)
u = copy(u0)
t = 0.
thist = [t]
uhist = [u0]
while t < tfinal
tnext = min(t+h, tfinal)
h = tnext - t
ydot1 = f(t, u)
ydot2 = f(t+h/2, u + h/2*ydot1)
ydot3 = f(t+h/2, u + h/2*ydot2)
ydot4 = f(t+h, u + h*ydot3)
u += h/6 * (ydot1 + 2*ydot2 + 2*ydot3 + ydot4)
t = tnext
push!(thist, t)
push!(uhist, u)
end
thist, hcat(uhist...)
end
```

```
ode_rk4 (generic function with 1 method)
```

```
x = LinRange(-1, 1, 51)[1:end-1]
u0 = exp.(-9 * x .^ 2)
@time thist, uhist = ode_rk4(advect_center, u0, h=.1, tfinal=1);
nsteps = size(uhist, 2)
plot(x, uhist[:, 1:(nsteps÷8):end])
```

```
0.168826 seconds (628.29 k allocations: 35.002 MiB, 12.00% gc time, 99.66% compilation time)
```

# Midpoint/trapezoid and \(\theta\) methods¶

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=.8
plot_stability(z -> Rz_theta(z, theta), "\$\\theta=$theta\$")
```

# A \(\theta\) solver¶

```
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]
thist, uhist = ode_theta_linear(A, [0., 1], h=.1, theta=1, tfinal=20)
scatter(thist, uhist')
plot!([cos, sin])
```

# Stiff decay to cosine¶

```
k = 10000
thist, uhist = ode_theta_linear(-k, [.2], forcing=t -> k*cos(t), tfinal=5, h=.5, theta=1)
scatter(thist, uhist[1,:])
plot!(cos)
```

# Advection as 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])
```

```
1.565558 seconds (4.13 M allocations: 235.608 MiB, 5.64% gc time, 99.91% compilation time)
```

# 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?

```
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])
```

```
0.046499 seconds (93.74 k allocations: 6.951 MiB, 98.56% compilation time)
```

# 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

# 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))
```