# L-24 MCS 471 Mon 17 Oct 2022 : richardson.jl

# Illustration of Richardson extrapolation,
# on forward differences.

using Printf
Base.show(io::IO, f::Float64) = @printf(io, "%.8e", f)

"""
    richardson(f::Function,z::Float64,h::Float64,n::Int)

returns the triangular table of numerical approximations of
the derivative of f at z, with initial h value and n steps.
"""
function richardson(f::Function,z::Float64,h::Float64,n::Int)
    R = zeros(n, n)
    w = h
    for i=1:n
        R[i, 1] = (f(z+w) - f(z))/w
        w = w/2
    end
    for j=2:n
        for i=j:n
            w = 2^(j-1) - 1
            R[i, j] = (2^(j-1)*R[i, j-1] - R[i-1, j-1])/w
        end
    end
    return R
end

"""
    Runs a test on Richardson extrapolation.
"""
function main()
    n = 4
    R = richardson(exp,1.0,0.1,n)
    println("Richardson extrapolation on exp(x) :")
    show(stdout, "text/plain", R); println("");
    E = zeros(n, n)
    for i=1:n
        for j=1:i
            E[i, j] = abs(R[i, j] - exp(1.0))
        end
    end
    println("The table with errors :")
    show(stdout, "text/plain", E); println("");
end

main()