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

# Experiments with Runge-Kutta methods on the pendulum problem.
# For the gravitational constant, we take 9.807.

using Printf

"""
    eulermodpend(p::Int,n::Int,verbose::Bool=true)

applies the modified Euler's method with n steps
on the interval [0,p*2*pi] on the pendulum problem.
"""
function eulermodpend(p::Int,n::Int,verbose::Bool=true)
    h = 2*p*pi/n
    if verbose
        println("  i   t         y1          y2 ")
    end
    (rt, ry1, ry2) = (zeros(n+1), zeros(n+1), zeros(n+1))
    y1 = pi/4
    y2 = 0
    (rt[1], ry1[1], ry2[1]) = (0, y1, y2)
    for i=1:n
        t = i*h
        (y1a, y2a) = (y1 + h*y2, y2 + h*(-9.807*sin(y1)))
        (y1b, y2b) = (y1, y2)
        y1 = y1 + (h/2)*(y2a + y2b)
        y2 = y2 + (h/2)*(-9.807*sin(y1a)-9.807*sin(y1b))
        (rt[i+1], ry1[i+1], ry2[i+1]) = (t, y1, y2)
        if verbose
            stri = @sprintf("%3d", i)
            strt = @sprintf("%.2f", t)
            stry1 = @sprintf("%.6e", y1)
            stry2 = @sprintf("%.6e", y2)
            println("$stri  $strt  $stry1  $stry2")
        end
    end
    return (rt, ry1, ry2)
end

"""
    rk2pend(p::Int,n::Int,verbose::Bool=true)

applies a 2-stage Runge-Kutta method with n steps
on the interval [0,p*2*pi] on the pendulum problem.
"""
function rk2pend(p::Int,n::Int,verbose::Bool=true)
    h = 2*p*pi/n
    if verbose
        println("  i   t         y1          y2 ")
    end
    (rt, ry1, ry2) = (zeros(n+1), zeros(n+1), zeros(n+1))
    y1 = pi/4
    y2 = 0
    (rt[1], ry1[1], ry2[1]) = (0, y1, y2)
    for i=1:n
        t = i*h
        k1a = y2
        k1b = -9.807*sin(y1)
        k2a = y2 + h*k1b
        k2b = -9.807*sin(y1 + h*k1a)
        y1 = y1 + (h/2)*k1a + (h/2)*k2a
        y2 = y2 + (h/2)*k1b + (h/2)*k2b
        (rt[i+1], ry1[i+1], ry2[i+1]) = (t, y1, y2)
        if verbose
            stri = @sprintf("%3d", i)
            strt = @sprintf("%.2f", t)
            stry1 = @sprintf("%.6e", y1)
            stry2 = @sprintf("%.6e", y2)
            println("$stri  $strt  $stry1  $stry2")
        end
    end
    return (rt, ry1, ry2)
end

"""
    rk3pend(p::Int,n::Int,verbose::Bool=true)

applies a 3-stage Runge-Kutta method with n steps
on the interval [0,p*2*pi] on the pendulum problem.
"""
function rk3pend(p::Int,n::Int,verbose::Bool=true)
    h = 2*p*pi/n
    if verbose
        println("  i   t         y1          y2 ")
    end
    (rt, ry1, ry2) = (zeros(n+1), zeros(n+1), zeros(n+1))
    y1 = pi/4
    y2 = 0
    (rt[1], ry1[1], ry2[1]) = (0, y1, y2)
    for i=1:n
        t = i*h
        k1a = y2
        k1b = -9.807*sin(y1)
        k2a = y2 + h*k1b/2
        k2b = -9.807*sin(y1 + h*k1a/2)
        k3a = y2 + 3*h*k2b/4
        k3b = -9.807*sin(y1 + 3*h*k2a/4)
        y1 = y1 + (h/9)*(2*k1a + 3*k2a + 4*k3a)
        y2 = y2 + (h/9)*(2*k1b + 3*k2b + 4*k3b)
        (rt[i+1], ry1[i+1], ry2[i+1]) = (t, y1, y2)
        if verbose
            stri = @sprintf("%3d", i)
            strt = @sprintf("%.2f", t)
            stry1 = @sprintf("%.6e", y1)
            stry2 = @sprintf("%.6e", y2)
            println("$stri  $strt  $stry1  $stry2")
        end
    end
    return (rt, ry1, ry2)
end

"""
    rk4pend(p::Int,n::Int,verbose::Bool=true)

applies a 4-stage Runge-Kutta method with n steps
on the interval [0,p*2*pi] on the pendulum problem.
"""
function rk4pend(p::Int,n::Int,verbose::Bool=true)
    h = 2*p*pi/n
    if verbose
        println("  i   t         y1          y2 ")
    end
    (rt, ry1, ry2) = (zeros(n+1), zeros(n+1), zeros(n+1))
    y1 = pi/4
    y2 = 0
    (rt[1], ry1[1], ry2[1]) = (0, y1, y2)
    for i=1:n
        t = i*h
        k1a = y2
        k1b = -9.807*sin(y1)
        k2a = y2 + (h/2)*k1b
        k2b = -9.807*sin(y1 + (h/2)*k1a)
        k3a = y2 + (h/2)*k2b
        k3b = -9.807*sin(y1 + (h/2)*k2a)
        k4a = y2 + h*k3b
        k4b = -9.807*sin(y1 + h*k3a)
        y1 = y1 + (h/6)*(k1a + 2*k2a + 2*k3a + k4a)
        y2 = y2 + (h/6)*(k1b + 2*k2b + 2*k3b + k4b)
        (rt[i+1], ry1[i+1], ry2[i+1]) = (t, y1, y2)
        if verbose
            stri = @sprintf("%3d", i)
            strt = @sprintf("%.2f", t)
            stry1 = @sprintf("%.6e", y1)
            stry2 = @sprintf("%.6e", y2)
            println("$stri  $strt  $stry1  $stry2")
        end
    end
    return (rt, ry1, ry2)
end

"""
Runs the modified Euler method and Runge-Kutta methods.
"""
function main()
    (p, n) = (1, 24)
    println("Running the modified Euler method ...")
    t, y1, y2 = eulermodpend(p,p*n,false)
    for k=1:p+1
        idx = (k-1)*n + 1
        stri = @sprintf("%5d", idx)
        strt = @sprintf("%5.2f", t[idx])
        stry1 = @sprintf("%13.6e", y1[idx])
        stry2 = @sprintf("%13.6e", y2[idx])
        serr = @sprintf("%.2e", abs(y1[idx] - y1[1]))
        println("$stri  $strt  $stry1  $stry2  $serr")
    end
    println("Running a 2-stage Runge-Kutta method ...")
    t, y1, y2 = rk2pend(p,p*n,false)
    for k=1:p+1
        idx = (k-1)*n + 1
        stri = @sprintf("%5d", idx)
        strt = @sprintf("%5.2f", t[idx])
        stry1 = @sprintf("%13.6e", y1[idx])
        stry2 = @sprintf("%13.6e", y2[idx])
        serr = @sprintf("%.2e", abs(y1[idx] - y1[1]))
        println("$stri  $strt  $stry1  $stry2  $serr")
    end
    println("Running a 3-stage Runge-Kutta method ...")
    t, y1, y2 = rk3pend(p,p*n,false)
    for k=1:p+1
        idx = (k-1)*n + 1
        stri = @sprintf("%5d", idx)
        strt = @sprintf("%5.2f", t[idx])
        stry1 = @sprintf("%13.6e", y1[idx])
        stry2 = @sprintf("%13.6e", y2[idx])
        serr = @sprintf("%.2e", abs(y1[idx] - y1[1]))
        println("$stri  $strt  $stry1  $stry2  $serr")
    end
    println("Running a 4-stage Runge-Kutta method ...")
    t, y1, y2 = rk4pend(p,p*n,false)
    for k=1:p+1
        idx = (k-1)*n + 1
        stri = @sprintf("%5d", idx)
        strt = @sprintf("%5.2f", t[idx])
        stry1 = @sprintf("%13.6e", y1[idx])
        stry2 = @sprintf("%13.6e", y2[idx])
        serr = @sprintf("%.2e", abs(y1[idx] - y1[1]))
        println("$stri  $strt  $stry1  $stry2  $serr")
    end
end

main()