# L-30 MCS 471 Mon 31 Oct 2022 : rkexp.jl

# Experiments with Runge-Kutta methods
# on the test equation y' = y, y(0) = 1. 

using Printf

"""
    modeulerexp(n::Int64,verbose::Bool=true)

The Modified Euler's method with 
n steps on the interval [0,1] on 
the test equation y' = y, y(0) = 1.
"""
function modeulerexp(n::Int64,verbose::Bool=true)
    h = 1.0/n
    if verbose
        println("  i   x       predict      correct        exact       error")
    end
    y0 = 1.0
    y1 = 1.0
    for i=1:n
        x = i*h
        y = (1.0+h)*y0
        y1 = y0 + h*(y0 + y)/2.0
        exact = exp(x)
        if verbose
            stri = @sprintf("%3d", i)
            strx = @sprintf("%.2f", x)
            stry = @sprintf("%.6e", y)
            stry1 = @sprintf("%.6e", y1)
            strexp = @sprintf("%.6e", exact)
            strerr = @sprintf("%.2e", abs(y1-exact))
            println("$stri  $strx  $stry  $stry1  $strexp  $strerr")
        end
        y0 = y1
    end
    return y1
end

"""
    rk2exp(n::Int64,verbose::Bool=true)

The Modified Euler's method as a 2-stage Runge-Kutta
method with n steps on the interval [0,1] on 
the test equation y' = y, y(0) = 1.
"""
function rk2exp(n::Int64,verbose::Bool=true)
    h = 1.0/n
    if verbose
        println("  i   x         k2        2-stage RK      exact       error")
    end
    y0 = 1.0
    y1 = 1.0
    for i=1:n
        x = i*h
        k1 = y0
        k2 = y0 + h*k1
        y1 = y0 + (h/2)*(k1 + k2)
        exact = exp(x)
        if verbose
            stri = @sprintf("%3d", i)
            strx = @sprintf("%.2f", x)
            stry = @sprintf("%.6e", k2)
            stry1 = @sprintf("%.6e", y1)
            strexp = @sprintf("%.6e", exact)
            strerr = @sprintf("%.2e", abs(y1-exact))
            println("$stri  $strx  $stry  $stry1  $strexp  $strerr")
        end
        y0 = y1
    end
    return y1
end

"""
    rk3exp(n::Int64,verbose::Bool=true)

A 3-stage Runge-Kutta method with n steps on 
the interval [0,1] on the test equation y' = y, y(0) = 1.
"""
function rk3exp(n::Int64,verbose::Bool=true)
    h = 1.0/n
    if verbose
        println("  i   x         k3        3-stage RK      exact       error")
    end
    y0 = 1.0
    y1 = 1.0
    for i=1:n
        x = i*h
        k1 = y0
        k2 = y0 + h*k1/2
        k3 = y0 + 3*h*k2/4
        y1 = y0 + (h/9)*(2*k1 + 3*k2 + 4*k3)
        exact = exp(x)
        if verbose
            stri = @sprintf("%3d", i)
            strx = @sprintf("%.2f", x)
            stry = @sprintf("%.6e", k3)
            stry1 = @sprintf("%.6e", y1)
            strexp = @sprintf("%.6e", exact)
            strerr = @sprintf("%.2e", abs(y1-exact))
            println("$stri  $strx  $stry  $stry1  $strexp  $strerr")
        end
        y0 = y1
    end
    return y1
end

"""
    rk4exp(n::Int64,verbose::Bool=true)

A 4-stage Runge-Kutta method with n steps on 
the interval [0,1] on the test equation y' = y, y(0) = 1.
"""
function rk4exp(n::Int64,verbose::Bool=true)
    h = 1.0/n
    if verbose
        println("  i   x         k4        4-stage RK      exact       error")
    end
    y0 = 1.0
    y1 = 1.0
    for i=1:n
        x = i*h
        k1 = y0
        k2 = y0 + (h/2)*k1
        k3 = y0 + (h/2)*k2
        k4 = y0 + h*k3
        y1 = y0 + (h/6)*(k1 + 2*k2 + 2*k3 + k4)
        exact = exp(x)
        if verbose
            stri = @sprintf("%3d", i)
            strx = @sprintf("%.2f", x)
            stry = @sprintf("%.6e", k4)
            stry1 = @sprintf("%.6e", y1)
            strexp = @sprintf("%.6e", exact)
            strerr = @sprintf("%.2e", abs(y1-exact))
            println("$stri  $strx  $stry  $stry1  $strexp  $strerr")
        end
        y0 = y1
    end
    return y1
end

"""
Runs the modified Euler method
and Runge-Kutta methods with 2, 3, and 4 stages.
"""
function main()
    println("Running the modified Euler method ...")
    modeulerexp(10)
    println("Running a 2-stage Runge-Kutta method ...")
    rk2exp(10)
    println("Running a 3-stage Runge-Kutta method ...")
    rk3exp(10)
    println("Running a 4-stage Runge-Kutta method ...")
    rk4exp(10)
end

main()