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

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

df_adapt (generic function with 1 method)

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]

8-element Vector{Float64}:
 -0.9999999999999948
 -2.220446049250313e-15
  0.9999999998019038
  5.440092820663267e-13
 -1.0
  0.0
  1.0
  0.0

# 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

romberg (generic function with 1 method)

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]

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.