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

# Illustrates adaptive integration with the composite trapezoidal rule,
# applied so previous function values are recycled.

using Printf

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

Returns a vector of at most n approximations 
for the definite integral of the function f 
over the interval [a,b], the i-th entry in the 
returned vector uses 2^i function evaluations.
Stops when the difference between two 
consecutive approximations is less than tol.

Example:

    t = adaptrap(cos,0,pi/2,10,1.e-5)
"""
function adaptrap(f::Function,
                  a::Float64,b::Float64,
                  n::Int64,tol::Float64)
    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]
        if(abs(t[i] - t[i-1])) < tol
            return t[1:i], i
        end 
        m = 2*m
    end
    return t, n
end

"""
Applies adaptive integration.
"""
function main()
    exact = 1.0
    n = 20
    tol = 1.0e-4
    t, nit = adaptrap(cos,0.0,pi/2,n,tol)
    println("The composite trapezoidal rule :")
    for i=1:length(t)
        strerr = @sprintf("%.2e", abs(exact - t[i]))
        strapp = @sprintf("%.16e", t[i])
        println("$strapp  $strerr")
    end
    print("The number of iterations : ")
    println(nit)
end

main()