Richardson extrapolation¶
Suppose we have a finite difference formula, say central difference...
$$ \begin{align*} L &= \phi(h) + a_2 \, h^2 + \mathcal{O}(h^4) \\ L &= \phi(h/2) + \frac{1}{2^2} a_2 \, h^2 + \mathcal{O}(h^4) \\ => L &= \frac{2^2 \phi(h/2) - \phi(h) }{2^2 - 1} + \mathcal{O}(h^4) \end{align*} $$We thus get a new formula for $L$, and we repeat this trick (with step size $h/2^2$)
differentiation: center diff x adaptive scheme¶
In [1]:
function df_adapt(f, x; Max=20, tol=1E-8, dx0=0.20)
function center(n)
# for internal points
dx = dx0 / 2^(n-1)
return (f(x+dx)-f(x-dx))/2/dx
end
Rmat = zeros(Max, Max)
# initial result
Rmat[1, 1] = center(1)
estim = Rmat[1, 1]
for i1 in 2:Max
Rmat[i1, 1] = center(i1)
for i2 in 2:i1
Rmat[i1, i2] = (4^(i2-1)*Rmat[i1, i2-1]-Rmat[i1-1, i2-1])/(4^(i2-1)-1)
end
if abs(estim-Rmat[i1,i1])*(4^(i1-1)-1) < tol || i1 == Max
i1 == Max ? println("max iteration reached...") : nothing
break
else
estim = Rmat[i1, i1]
end
end
return estim
end
Out[1]:
df_adapt (generic function with 1 method)
In [2]:
f = x -> sin(x)
df = x -> df_adapt(f, x)
d2f = x -> df_adapt(df, x)
d3f = x -> df_adapt(d2f, x)
d4f = x -> df_adapt(d3f, x)
[df(pi), d2f(pi), d3f(pi), d4f(pi), -1.0, 0.0, 1.0, 0.0]
Out[2]:
8-element Vector{Float64}: -0.9999999999999948 -2.220446049250313e-15 0.9999999998019038 5.440092820663267e-13 -1.0 0.0 1.0 0.0
In [3]:
# Romberg Integration
function romberg(f, a, b; Max=20, tol=1E-10)
function interior(n)
# for internal points
dx = (b - a) / 2^(n-1)
return dx * sum([f(a + (2*i1-1)*dx) for i1 in 1:2^(n-2)])
end
Rmat = zeros(Max, Max)
# initial result
Rmat[1, 1] = 0.5 * (b-a) * (f(a)+f(b))
estim = Rmat[1, 1]
for i1 in 2:Max
# without duplicated calculations
Rmat[i1, 1] = 0.5 * Rmat[i1-1, 1] + interior(i1)
for i2 in 2:i1
# Rmat[i1, i2] = Rmat[i1, i2-1] + (Rmat[i1, i2-1]-Rmat[i1-1, i2-1])/(4^(i2-1)-1)
Rmat[i1, i2] = (4^(i2-1)*Rmat[i1, i2-1]-Rmat[i1-1, i2-1])/(4^(i2-1)-1)
end
if abs(estim-Rmat[i1,i1])*(4^(i1-1)-1) < tol || i1 == Max
i1 == Max ? println("max iteration reached...") : nothing
break
else
estim = Rmat[i1, i1]
end
end
return estim
end
Out[3]:
romberg (generic function with 1 method)
In [4]:
using BenchmarkTools
println([romberg(x->x*log(1+x), 0.0, 1.0), 0.25])
@benchmark ret1 = romberg(x->x*log(1+x), 0.0, 1.0), 0.25
[0.25000000000000006, 0.25]
Out[4]:
BenchmarkTools.Trial: 10000 samples with 10 evaluations. Range (min … max): 1.017 μs … 485.258 μs ┊ GC (min … max): 0.00% … 99.28% Time (median): 1.250 μs ┊ GC (median): 0.00% Time (mean ± σ): 1.451 μs ± 9.235 μs ┊ GC (mean ± σ): 12.67% ± 1.99% ▄▁█▇▅▃ ▁▂▃▃▂▁▁▁▂▂▂▂▂▃▄████████▅▆▆▅▄▃▂▂▂▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁▁ ▂ 1.02 μs Histogram: frequency by time 1.71 μs < Memory estimate: 4.69 KiB, allocs estimate: 8.