# 2022-09-26 Stiffness#

## Last time#

• Stability diagrams

• Energy dissipation

• Exploring the $$\theta$$ method

• PDE as ODE

## Today#

• $$A$$- and $$L$$-stability

• Spatial, temporal, and physical dissipation

• Stiffness

using Plots
default(linewidth=3)
using LinearAlgebra
using SparseArrays

function plot_stability(Rz, title; xlims=(-3, 3), 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

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

Rz_theta(z, theta) = (1 + (1 - theta)*z) / (1 - theta*z)

Rz_theta (generic function with 1 method)


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

10-element Vector{ComplexF64}:
0.0 - 2.938926261462368im
0.0 - 2.9389262614623664im
0.0 - 1.2560848302247556e-16im
0.0 + 1.2560848302247556e-16im
0.0 + 2.9389262614623664im
0.0 + 2.938926261462368im
4.81482486096809e-34 - 4.755282581475769im
4.81482486096809e-34 + 4.755282581475769im
3.0814879110195774e-33 - 4.7552825814757655im
3.0814879110195774e-33 + 4.7552825814757655im

n = 50
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=0.5, tfinal=1.);
nsteps = size(uhist, 2)
plot(x, uhist[:, 1:(nsteps÷8):end])

  0.035987 seconds (78.89 k allocations: 6.054 MiB, 92.42% 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
eigvals(Matrix(heat_matrix(8)))

8-element Vector{Float64}:
-64.0
-54.62741699796952
-54.6274169979695
-31.999999999999996
-31.99999999999998
-9.372583002030481
-9.37258300203047
1.4210854715202004e-14

n = 100
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=.005, theta=.5, tfinal=.1);
nsteps = size(uhist, 2)
plot(x, uhist[:, 1:5])

  0.033819 seconds (82.24 k allocations: 5.521 MiB, 95.34% compilation time)


# Stability classes and the $$\theta$$ method#

## Definition: $$A$$-stability#

A method is $$A$$-stable if the stability region

$\{ z : |R(z)| \le 1 \}$
contains the entire left half plane
$\Re[z] \le 0 .$
This means that the method can take arbitrarily large time steps without becoming unstable (diverging) for any problem that is indeed physically stable.

## Definition: $$L$$-stability#

A time integrator with stability function $$R(z)$$ is $$L$$-stable if

$\lim_{z\to\infty} R(z) = 0 .$
For the $$\theta$$ method, we have
$\lim_{z\to \infty} \frac{1 + (1-\theta)z}{1 - \theta z} = \frac{1-\theta}{\theta} .$
Evidently only $$\theta=1$$ is $$L$$-stable.

# Spectrum of operators#

theta=0.5
h = .05
plot_stability(z -> Rz_theta(z, theta), "\$\\theta=$theta, h=$h\$")
scatter!(real(ev), imag(ev))

theta=0
h = .001
plot_stability(z -> Rz_theta(z, theta), "\$\\theta=$theta, h=$h\$")
ev = eigvals(Matrix(h*heat_matrix(40)))
scatter!(real(ev), imag(ev))


(17)#\begin{align} (-\kappa u_x + w u)_x &= 0 & Pe_L &= \frac{L \lvert w \rvert}{\kappa} \end{align}
function advdiff_matrix(n; kappa=1, wind=1, upwind=0.)
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+1))
diffuse = [-1, 2, -1] * kappa / dx^2
advect_upwind = [-1, 1, 0] * wind / dx
advect_center = [-1, 0, 1] * wind / 2dx
stencil = -diffuse - upwind * advect_upwind - (1 - upwind) * advect_center
append!(vals, stencil)
end
sparse(rows, cols, vals)
end

5×5 SparseMatrixCSC{Float64, Int64} with 15 stored entries:
-1.25   -0.625    ⋅       ⋅      1.875
1.875  -1.25   -0.625    ⋅       ⋅
⋅      1.875  -1.25   -0.625    ⋅
⋅       ⋅      1.875  -1.25   -0.625
-0.625    ⋅       ⋅      1.875  -1.25

n = 50
x = LinRange(-1, 1, n+1)[1:end-1]
u0 = exp.(-16 * x .^ 2)
thist, uhist = ode_theta_linear(A, u0, h=.04, theta=1)
plot(x, uhist[:, 1:5])


# Explain the mechanism and obstacles to fixing#

A = advdiff_matrix(50, kappa=.02, upwind=0)
thist, uhist = ode_theta_linear(A, u0, h=.04, theta=.5, tfinal=1.);
plot(x, uhist[:, 1:5])

A = advdiff_matrix(50, kappa=0, upwind=0)
thist, uhist = ode_theta_linear(A, u0, h=.04, theta=1, tfinal=1.)
plot(x, uhist[:, 1:5])


# Explain the mechanism#

A = advdiff_matrix(50, kappa=.02, upwind=0)
thist, uhist = ode_theta_linear(A, u0, h=.04, theta=.5, tfinal=1.);
plot(x, uhist[:, 1:5])

n = 50
x = LinRange(-1, 1, n)
u0 = exp.(-16 * x .^ 2)
thist, uhist = ode_theta_linear(A, u0, h=.04, theta=1, tfinal=1.)
plot(x, uhist[:, 1:20])


# Stability limits: CFL principle for advection#

A stable method must be capable of propagating information at least as fast as the continuous system it discretizes.

• Courant, Friedrichs, Levy (1928)

• Discrete information moves at most $$\Delta x$$ in time $$\Delta t$$ (with explicit integrator)

• Physical information moves distance $$|w| \Delta t$$

• “CFL number” or “Courant number”

$\textrm{CFL} = \frac{\Delta t |w|}{\Delta x} \le 1$

• What does this mean for implicit methods?

theta=0
n = 20
dx = 2 / n
cfl = .5
h = cfl * dx / 1

plot_stability(z -> Rz_theta(z, theta), "\$\\theta=$theta, h=$h, \\mathrm{CFL}=$cfl\") ev = eigvals(Matrix(h*advdiff_matrix(n, kappa=0, upwind=1))) scatter!(real(ev), imag(ev))  # Stability limits: explicit diffusion# • The diffusion operator has stencil $\frac{\kappa}{(\Delta x)^2} \begin{bmatrix} 1 & -2 & 1 \end{bmatrix}$ • Apply this to test function $$e^{i\theta x}$$ to get the symbol (18)#\begin{align} L_{\Delta x} e^{i\theta x} &= \frac{\kappa}{(\Delta x)^2} \Big[ e^{-i\theta\Delta x} - 2 + e^{i\theta\Delta x} \Big] e^{i\theta x} \\ &= \frac{2\kappa}{(\Delta x)^2} \underbrace{\Big[ \cos (\theta \Delta x) - 1 \Big]}_{\min = -2} e^{i\theta x} \end{align} • Conclude $$\frac{-4 \kappa}{(\Delta x)^2} \le \lambda \le 0$$, need $$-2 \le h\lambda \le 0$$ for Euler theta=0 n = 20 dx = 2 / n kappa = 1 lambda_min = -4 * kappa / dx^2 h = -2 / lambda_min plot_stability(z -> Rz_theta(z, theta), "\\\theta=$theta, h=$h\$") ev = eigvals(Matrix(h*advdiff_matrix(n, kappa=kappa, wind=0))) scatter!(real(ev), imag(ev))  # Stiffness# Stiff equations are problems for which explicit methods don’t work. (Hairer and Wanner, 2002) • “stiff” dates to Curtiss and Hirschfelder (1952) k = 5 thist, uhist = ode_theta_linear(-k, [.2], forcing=t -> k*cos(t), tfinal=5, h=.5, theta=1) scatter(thist, uhist[1,:]) plot!(cos)  function mms_error(h; theta=.5, k=10) u0 = [.2] thist, uhist = ode_theta_linear(-k, u0, forcing=t -> k*cos(t), tfinal=3, h=h, theta=theta) T = thist[end] u_exact = (u0 .- k^2/(k^2+1)) * exp(-k*T) .+ k*(sin(T) + k*cos(T))/(k^2 + 1) uhist[1,end] .- u_exact end  mms_error (generic function with 1 method)  hs = .5 .^ (1:8) errors = mms_error.(hs, theta=0.5, k=1000) plot(hs, norm.(errors), marker=:circle, xscale=:log10, yscale=:log10) plot!(hs, hs, label="\$h\$", legend=:topleft) plot!(hs, hs.^2, label="\$h^2\$")  ## Discuss: is advection-diffusion stiff?# theta=0 n = 80 dx = 2 / n kappa = .05 lambda_min = -4 * kappa / dx^2 cfl = 1 h = min(-2 / lambda_min, cfl * dx) plot_stability(z -> Rz_theta(z, theta), "\$\\theta=$theta, h=$h, Pe_h = $(h/kappa)\$")
scatter!(real(ev), imag(ev))


# Cost scaling#

• Spatial discretization with error $$O((\Delta x)^p)$$

• Time discretization with error $$O((\Delta t)^q)$$

p = 2   # spatial order of accuracy
q = 2   # temporal order of accuracy
Pe = 10 # Peclet number with wind=1

dt_accuracy(dx) = dx .^ (p/q)
dt_stability(dx) = min(dx, Pe * dx ^ 2)

plot(dx, [dt_accuracy.(dx) dt_stability.(dx)],
label=["accuracy" "stability"],
xscale=:log10, yscale=:log10, legend=:bottomright,
xlabel="\$\\Delta x\$", ylabel="\$\\Delta t\$")

Cannot convert Float64 to series data for plotting

Stacktrace:
[1] error(s::String)
@ Base ./error.jl:35
[2] _prepare_series_data(x::Float64)
@ RecipesPipeline ~/.julia/packages/RecipesPipeline/OXGmH/src/series.jl:8
[3] _series_data_vector(x::Float64, plotattributes::Dict{Symbol, Any})
@ RecipesPipeline ~/.julia/packages/RecipesPipeline/OXGmH/src/series.jl:35
[4] macro expansion
@ ~/.julia/packages/RecipesPipeline/OXGmH/src/series.jl:134 [inlined]
[5] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, #unused#::Type{RecipesPipeline.SliceIt}, x::Any, y::Any, z::Any)
@ RecipesPipeline ~/.julia/packages/RecipesBase/qpxEX/src/RecipesBase.jl:289
[6] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
@ RecipesPipeline ~/.julia/packages/RecipesPipeline/OXGmH/src/user_recipe.jl:36
[7] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
@ RecipesPipeline ~/.julia/packages/RecipesPipeline/OXGmH/src/RecipesPipeline.jl:70
[8] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
@ Plots ~/.julia/packages/Plots/lW9ll/src/plot.jl:209
[9] #plot#145
@ ~/.julia/packages/Plots/lW9ll/src/plot.jl:91 [inlined]
[10] top-level scope
@ In[231]:8
[11] eval
@ ./boot.jl:368 [inlined]
[12] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

n = 2 .^ (1:10); dx = 1 ./ n
function error((dx, dt))
dt <= dt_stability(dx) ? dx ^ p + dt ^ q : 10
end
plot(n, [error.(zip(dx, dt_accuracy.(dx))) error.(zip(dx, dt_stability.(dx))) dx.^p],
label=["accuracy" "stability" "spatial only"],
marker=:circle, xscale=:log10, yscale=:log10,
xlabel="number of grid points", ylabel="error")


# Work-precision (error vs cost)#

p = 1
q = 2
Pe = 10
dt = [dt_accuracy.(dx) dt_stability.(dx) dt_accuracy.(dx)]
c = [n./dt[:,1] n./dt[:,2] (n.^1.5)./dt[:,3]] # O(n^1.5 solve cost)
e = [error.(zip(dx, dt[:,1])) error.(zip(dx, dt[:,2])) dx.^p+dt[:,3].^q]
plot(c, e, xscale=:log10, yscale=:log10,
label=["explicit_accuracy" "explicit_stable" "implicit"], xlabel="cost", ylabel="error")

MethodError: no method matching +(::Float64, ::Vector{Float64})
Closest candidates are:
+(::Any, ::Any, ::Any, ::Any...) at operators.jl:591
+(::T, ::T) where T<:Union{Float16, Float32, Float64} at float.jl:383
+(::Union{Float16, Float32, Float64}, ::BigFloat) at mpfr.jl:414
...

Stacktrace:
[1] top-level scope
@ In[230]:6
[2] eval
@ ./boot.jl:368 [inlined]
[3] include_string(mapexpr::typeof(REPL.softscope), mod::Module, code::String, filename::String)

• Order of accuracy matters; usually balance spatial order $$p$$ and temporal order $$q$$