Fuzzy Logic for Verifying Numerical Order of a Method

Specific Domains Numerics
,
1

I often program numerical methods for solving initial value problems, and would like to test their order of convergence. However, I find that it is difficult to have automated convergence tests, due to the pre-asymptotic region and the roundoff-dominated region.

For example, consider the code below, which solves the Dahlquist test equation using the Trapezoidal rule (order 2) with various step sizes.

From the graph, we can clearly see that the order of the numerical method is 2, since the slope of the graph on a log10-log10 scale is 2.

One option for checking the order of convergence would be to fit a line to the results and check the slope, but depending on the width of the pre-asymptotic and roundoff-dominated regions, the slope may differ drastically from the the slopes in the asymptotic regime. It gets even worse if the numerical method is unstable for large stepsizes, since then the error blows up in the pre-asymptotic region, drastically changing the slope of the fitted line.

Another option is to compute the pairwise order of convergence. I.e., the order of convergence for each pair of stepsizes dt and 0.5*dt. If I expect the order of convergence to be 2, and I see a lot of pairwise orders close to 2, then the order of convergence probably really is 2. But again, if the pre-asymptotic and roundoff-dominated regions are wide, then there may be more “bad” pairwise orders (near zero or ridiculously large) than “good” ones. For very high-order methods especially, the asymptotic region will be small because the slop of the line is so steep.

Are there any packages which implement “fuzzy logic” for verifying the the numerical order of a method? E.g., if I provide a function f(dt) which outputs the error when using stepsize dt, I can tell with high-probability whether the order of f is N?

# Solve the dahlquist initial value problem
#   u̇ = λu,    u(0) = u₀,     0 ≤ t ≤ T
# Using the Trapezoidal rule for various stepsizes Δt.

using CairoMakie

const λ = 1im
const T = 100.0
const u0 = 1.0+0.0im
const uf_analytic = u0*exp(λ*T)

function solve_ivp(nsteps::Integer)
    Δt = T / nsteps 
    u = u0 
    for n in 1:nsteps
        u = (1+0.5*Δt*λ)*u / (1-0.5*Δt*λ)
    end
    return u
end

results = Tuple{Int64, ComplexF64}[]
for nsteps in 2 .^ (0:30)
    uf = solve_ivp(nsteps)
    push!(results, (nsteps, uf))
end

nsteps_vec = first.(results)
Δts_vec = T ./ nsteps_vec
ufs_vec = last.(results)
errors_vec = abs.(ufs_vec .- uf_analytic)

fig = Figure()
ax = Axis(fig[1, 1],
    xscale = log10,
    yscale = log10,
    xlabel = "Δt",
    ylabel = "Final State Error",
)
lines!(ax, Δts_vec, errors_vec)
fig

errorsaturation
errorsaturation1200×900 36.4 KB