Solving linear ODE using `expv`

expv-example.jl

using OrdinaryDiffEq, ExponentialUtilities
using LinearAlgebra, SparseArrays

# Example: Use expv to solve constant-coefficient linear systems
#   u' = ϵu_xx + u, u ∈ [0, 1], t ∈ [0, 1], Dirichlet BC
# The analytical solution is
#   u(t) = exp(tA)u(0), A = L + I
# which we approximate with expv/expv_timestep
N = 100; T = 1.0
dx = 1 / (N + 1)
xs = (1:N) * dx
f0 = x -> sin(3*pi*x)^3
u0 = f0.(xs)

dd = -2 * ones(N)
du = ones(N - 1)
ϵ = 0.01
A = spdiagm(-1 => du, 0 => dd, 1 => du) * (ϵ / dx^2) + spdiagm(0 => ones(N))

## Reference solution using OrdinaryDiffEq
prob = ODEProblem((du,u,p,t) -> mul!(du, A, u), u0, (0.0, T))
sol_ref = solve(prob, Tsit5(); abstol=1e-14, reltol=1e-14)(T)

## "Exact" solution by computing the matrix exponential directly.
## This is analogous to caching ExpRK algorithms.
sol1 = exp(T * Matrix(A)) * u0 # slow
err1 = norm(sol1 - sol_ref)
println("Caching: error = ", err1) # 6.238776757814508e-14

## Approximate using Krylov expv. The signature is
##   expv(t, A, b) -> exp(t*A) * b
## Pass in keyword argument `m` (Krylov dimension) to control the accuracy
## of Krylov approximation.
sol2 = expv(T, A, u0; m=5)
err2 = norm(sol2 - sol_ref)
println("Krylov expv, m = 5: error = ", err2) # 1.1338086880671259e-13

## One way to increase accuracy is to not do expv in one batch, but instead
## subdivide the time interval. To do this, use the adaptive `expv_timestep`.
sol3 = expv_timestep(T, A, u0; tol=1e-14)
err3 = norm(sol3 - sol_ref)
println("Adaptive expv, tol = 1e-14: error = ", err3) # 1.373287832909154e-14