Convergence of AB4¶
We consider the IVP \(u'=\sin[(u+t)^2]\) over \(0\le t \le 4\), with \(u(0)=-1\).
using FundamentalsNumericalComputation
ivp = ODEProblem((u,p,t)->sin((t+u)^2),-1.,(0.0,4.0))
[36mODEProblem[0m with uType [36mFloat64[0m and tType [36mFloat64[0m. In-place: [36mfalse[0m
timespan: (0.0, 4.0)
u0: -1.0
We use a solver from DifferentialEquations to construct an accurate approximation to the exact solution.
u_exact = solve(ivp,Tsit5(),reltol=1e-14,abstol=1e-14);
Now we perform a convergence study of the AB4 code.
n = @. 10*2^(0:5)
err = zeros(size(n))
for (j,n) = enumerate(n)
t,u = FNC.ab4(ivp,n)
err[j] = norm(u_exact.(t)-u,Inf)
end
pretty_table((n=n,err=err),["n","inf-norm error"],backend=:html)
| n | inf-norm error |
|---|---|
| 10 | 1.421332157914805 |
| 20 | 0.2998677236341702 |
| 40 | 0.006278094753455532 |
| 80 | 0.0005392729579345446 |
| 160 | 3.975975730263759e-5 |
| 320 | 2.64515824960343e-6 |
The method should converge as \(O(h^4)\), so a log-log scale is appropriate for the errors.
plot(n,err,m=:o,label="AB4",
xaxis=(:log10,"n"),yaxis=(:log10,"inf-norm error"),
title="Convergence of AB4",leg=:bottomleft)
plot!(n,0.1*(n/n[1]).^(-4),l=:dash,label="4th order")