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

2

I haven’t done it yet but these days it pretty doable to bring multi-modal LLM in the automation loop to decide the range to do regression. Do you mind to introduce LLM to your tool?

3

I guess that world work fine, but the end result is that there is some hard-coded range for the regression in the test, which makes the test more brittle to changes in the initial value problem being used (e.g., if I change lambda in the dahlquist test equation, I have to adjust the range in the test).

So I hoped maybe somebody had a more robust solution. If not, LLM is acceptable.

4

One option is to use Richardson extrapolation of error(Δt)/error(Δt/2) as \Delta t \to 0^+. The Richardson.jl package tries to automatically detect when things become roundoff-limited at small \Delta t — this is possible because roundoff errors are noisy and no longer polynomial. The key point here is that you know the exact error should be asymptotically polynomial in \Delta t, which allows one to bring a lot of analytical tools to bear.

For example, here is some code to automatically determine the order n of accuracy of a forward finite-difference approximation \frac{f(x+\Delta x) - f(x)}{\Delta x} = f'(x) + O(\Delta x^n), which should be n=1, by testing it with f(x) = \sin(x) to compute f'(1) = \cos(1). If we compute the error directly, it is only approximately linear for a limited range of \Delta x values:

image
image580×458 17.7 KB

However, we can use Richardson extrapolation to estimate \lim_{\Delta x \to 0^+} \frac{\text{error}(\Delta x)}{\text{error}(\Delta x/2)} = 2^n = 2 starting at \Delta x = 0.1:

julia> using Richardson

julia> fd_error(x, Δx) = (sin(x+Δx)-sin(x))/Δx - cos(x)
fd_error (generic function with 1 method)

julia> extrapolate(0.1) do Δx
           @show Δx
           fd_error(1, Δx) / fd_error(1, Δx/2)
       end
Δx = 0.1
Δx = 0.0125
Δx = 0.0015625
Δx = 0.0001953125
(1.9999999972699891, 5.984715389928397e-9)

julia> log2(ans[1]) # the order of accuracy
0.9999999980307134

It correctly returns the limit of \approx 2 to 9 digits (≈ the default tolerance \sqrt{\varepsilon}) after evaluating at only four \Delta x values, and log2 of this gives the order of accuracy n \approx 1 to 8–9 digits.

It can also still work to start at a value that is vastly too large, like \Delta x = 10:

julia> extrapolate(10) do Δx
           @show Δx
           fd_error(1, Δx) / fd_error(1, Δx/2)
       end
Δx = 10.0
Δx = 1.25
Δx = 0.15625
Δx = 0.01953125
Δx = 0.00244140625
Δx = 0.00030517578125
(2.000000000619581, 1.3960816502844864e-8)

Though sometimes you may want to set breaktol=Inf and play with other parameters. For example, starting at the ridiculous \Delta x = 1000, it can still give us 11–12 digits:

julia> extrapolate(1000, contract=0.5, breaktol=Inf) do Δx
           @show Δx
           fd_error(1, Δx) / fd_error(1, Δx/2)
       end
Δx = 1000.0
Δx = 500.0
Δx = 250.0
Δx = 125.0
Δx = 62.5
Δx = 31.25
Δx = 15.625
Δx = 7.8125
Δx = 3.90625
Δx = 1.953125
Δx = 0.9765625
Δx = 0.48828125
Δx = 0.244140625
Δx = 0.1220703125
Δx = 0.06103515625
Δx = 0.030517578125
(2.000000000008055, 5.1660897781857784e-11)

(The contract=0.5 gives it more data points than the default contract=0.125, which may be important for high-order methods where the asymptotic region is small.)

Note also that Richardson.extrapolate gives you an error estimate (its second output), so you can tell that something is going terribly wrong if this estimate is large.

PS. I did a bit of hunting around, and apparently they use a similar Richardson-extrapolation metric in the CFD community, which they call a “grid-convergence index (GCI)” (Celik et al., 2008). However, this “GCI” is more primitive in that one chooses only three \Delta x values, non-adaptively, possibly because CFD calculations are so costly.

5

Thanks, I’ll try that package out! It looks like it does pretty much exactly what I want.

It’s not clear to me whether Richardson Extrapolation in particular enables the package’s ability to find the region of convergence, or if there is just good code for identifying the region of convergence alongside the code which implements Richardson Extrapolation. I’ll check it out!

6

The way it works is that (a) Richardson extrapolation provides an error estimate (by comparing subsequent orders of extrapolation) and (b) the underlying “Neville tableau” algorithm allows you to cheaply get extrapolants from every consecutive subsequence. So, it picks the subsequence with the lowest error estimate; this automatically identifies the region where the extrapolation best converges.

(I keep meaning to write up some notes on this, because a lot of the descriptions of Richardson extrapolation do not talk about its ability to be adaptive in this way.)

7

Thanks for the explanation, that is very helpful!