# L-26 MCS 471 Fri 21 Oct 2022 : romberg4pi.jl

# Illustrates adaptive integration with the composite trapezoidal rule,
# recycling previous function evaluations and Romberg integration.

using Printf

"""
    function adaptrap(f::Function,a::Float64,b::Float64,n::Int64)

Returns a vector of n approximations for the definite integral of
the function f over the interval [a,b], the i-th entry in the vector
uses 2^i function evaluations.

Example:

    t = adaptrap(cos,0,pi/2,10)
"""
function adaptrap(f::Function,a::Float64,b::Float64,n::Int64)
    t = zeros(n)
    h = (b-a)        # size of subinterval
    m = 1            # number of subintervals
    t[1] = (f(a) + f(b))*h/2
    for i = 2:n
        h = h/2
        for j=0:m-1
            t[i] = t[i] + f(a+h+j*2*h)
        end;
        t[i] = t[i-1]/2 + h*t[i]
        m = 2*m
    end
    return t
end

"""
    function romberg(t::Array{Float64,1})

Applies extrapolation to the approximations in t,
returned by the composite Trapezoidal rule.

Example:

    t = adaptrap(cos,0,pi/2,6)
    et = romberg(t);
"""
function romberg(t::Array{Float64,1})
    n = length(t)
    et = zeros(n,n)
    et[:,1] = t
    for i = 2:n
        for j = 2:i
            r = 4^(j-1)
            et[i,j] = (r*et[i,j-1] - et[i-1,j-1])/(r - 1);
       end
    end
    return et
end

"""
Applies Romberg integration to approximate pi.
The function sqrt(1 - x^2) is not sufficiently many times
continously differentiable on [0, 1] for Romberg to work
directly on the integral over [0, 1] and therefore,
the interval is shortened to [0, sqrt(2)/2].
"""
function main()
    f(x) = sqrt(1-x^2) - sqrt(2.0)/2.0
    exact = (pi - 2.0)/8.0
    n = 6
    t = adaptrap(f,0.0,sqrt(2.0)/2.0,n)
    println("The composite trapezium rule :")
    N = 2
    for i=1:n
        strnbr = @sprintf("%7d", N)
        strerr = @sprintf("%.2e", abs(exact - t[i]))
        strapp = @sprintf("%.16e", t[i])
        println("$strnbr  $strapp  $strerr")
        N = 2*N
    end
    et = romberg(t)
    println("Romberg integration :")
    for i=1:n
        strerr = @sprintf("%.2e", abs(exact - et[i,i]))
        strapp = @sprintf("%.16e", et[i,i])
        println("$strapp  $strerr")
    end
    println("Errors in the table ")
    for i=1:n
        for j=1:i
            strerr = @sprintf("%.2e", abs(exact - et[i,j]))
            print("  $strerr")
        end
        println("")
    end
end

main()