Incomplete implementation of the Gauss-Jackson Integrator

gauss-jackson.jl

"""
Backwards difference

∇fᵢ = (fᵢ - fᵢ₋₁)
"""
function ∇(power, f, i)
    if power > 1
        ∇(power-1, f, i) - ∇(power-1, f, i-1)
    else
        f[i] - f[i-1]
    end
end

"""
Inverse backwards difference

∇⁻¹fᵢ = (fᵢ + ∇⁻¹fᵢ₋₁)

We also pass an integration constant vector with one constant per power.
 """
function ∇⁻¹(power, f, i, C)
    if i < 1
        return C[power]
    end

    if power == 0
        return f[i]
    end

    ∇⁻¹(power-1, f, i) + ∇⁻¹(power, f, i-1)
end

function C(y₀, ẏ₀, Δt)
    local a = [ 14797/152064, -326911/39916800,
                55067/39916800, -2539/13305600, 317/22809600 ]
    local b = [ 0, -252769/3628800,
                68119/3628800, -1469/403200, 2497/7257600 ]

    [ ẏ₀/Δt-sum(map(k -> b[k]*f[k], 1:4)) + f[1]/2;
      y₀/Δt^2-sum(map(k -> a[k]*f[k], 1:4)) + 0 ]
end

function L(ÿ, ẏ₀, y₀, Δt)
    local table = [-1/240, -1/240, -221/60480]
    local t = length(table)
    local n = length(ÿ)

    local sum = ∇⁻¹(2, ÿ, n-1, C(y₀, ẏ₀, Δt)) + 1/12*ÿ[n]

    for i ∈ 1:n
        if i > t
            break
        end

        sum += table[i]*∇(i+1, ÿ, n)
    end

    sum
end

function corrector(ÿ, y₀, ẏ₀, Δt)
    (Δt^2) * L(ÿ)
end

function J(ÿ, ẏ₀, y₀, Δt)
    local table = [1/12, 19/240, 3/40]
    local t = length(table)
    local n = length(ÿ)

    local sum = ∇⁻¹(2, ÿ, n, C(y₀, ẏ₀, Δt)) + 1/12*ÿ[n]

    for i ∈ 1:n
        if i > t
            break
        end

        sum += table[i]*∇(i, ÿ, n)
    end

    sum
end

function predictor(ÿ, y₀, ẏ₀,  Δt)
    (Δt^2) * J(ÿ)
end

∇⁻¹(2, [-sin(0), -sin(.5), -sin(1), -sin(1.5)], 4),
[-sin(0), -sin(.5), -sin(1), -sin(1.5)]

#
# RK4 double integrator to check gauss-jackson results
#

d1(t, y, ẏ) = ẏ
d2(t, y, ẏ) = -sin(t)

function double_rk4(y, ẏ, Δt::Real, t::Real)
    kp₁ = d1(t, y, ẏ)
    kv₁ = d2(t, y, ẏ)
    kp₂ = d1(t+Δt/2, y + kp₁*Δt/2, ẏ + kv₁*Δt/2)
    kv₂ = d2(t+Δt/2, y + kp₁*Δt/2, ẏ + kv₁*Δt/2)
    kp₃ = d1(t+Δt/2, y + kp₂*Δt/2, ẏ + kv₂*Δt/2)
    kv₃ = d2(t+Δt/2, y + kp₂*Δt/2, ẏ + kv₂*Δt/2)
    kp₄ = d1(t+Δt, y + kp₃*Δt, ẏ + kv₃*Δt)
    kv₄ = d2(t+Δt, y + kp₃*Δt, ẏ + kv₃*Δt)

    y+Δt/6*(kp₁+2*kp₂+2*kp₃+kp₄), ẏ+Δt/6*(kv₁+2*kv₂+2*kv₃+kv₄)
end