IRK

irk.jl

using Printf, LinearAlgebra, StaticArrays

s3o6=sqrt(3)/6
const a= SA[1/4 1/4-s3o6; 1/4+s3o6 1/4]
const b= SA[1/2,1/2]
const c= SA[1/2-s3o6,1/2+s3o6]

function euler(t,y,h)
    return y+h*fF(t,y)
end

function rk2(t,y,h)
    k1=fF(t,y)
    k2=fF(t,y+h*2/3*k1)
    return y+h*(1/4*k1+3/4*k2)
end

@inbounds function fF(t,y)
    SA[
        sin(t)*y[2]*y[3] - y[1]*y[2]*y[3],
        (y[1]*y[3])/20 - sin(t)*y[1]*y[3],
        -(y[1]*y[2])/20 + y[1]^2*y[2]
    ]
end

@inbounds function fG(t,y,h,k)
    SA[
        k[1] - sin(t + h*c[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) + (h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]),
        k[2] - sin(t + h*c[1])*(-(h*(a[1, 1]*k[1] + a[1, 2]*k[4])) - y[1])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]) - ((h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]))/20,
        k[3] + ((h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]))/20 - (h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])^2*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]),
        k[4] - sin(t + h*c[2])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) + (h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]),
        k[5] - sin(t + h*c[2])*(-(h*(a[2, 1]*k[1] + a[2, 2]*k[4])) - y[1])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]) - ((h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]))/20,
        k[6] + ((h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]))/20 - (h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])^2*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]),
    ]
end

@inbounds function fDG(t,y,h,k)
    DG= @MMatrix zeros(6,6)
DG[1, 1] = 1 + h*a[1, 1]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])
DG[1, 2] = -(sin(t + h*c[1])*h*a[1, 1]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])) + h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])
DG[1, 3] = -(sin(t + h*c[1])*h*a[1, 1]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])) + h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])
DG[1, 4] = h*a[1, 2]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])
DG[1, 5] = -(sin(t + h*c[1])*h*a[1, 2]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])) + h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])
DG[1, 6] = -(sin(t + h*c[1])*h*a[1, 2]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])) + h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])
DG[2, 1] = -(h*a[1, 1]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]))/20 + sin(t + h*c[1])*h*a[1, 1]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])
DG[2, 2] = 1
DG[2, 3] = -(sin(t + h*c[1])*h*a[1, 1]*(-(h*(a[1, 1]*k[1] + a[1, 2]*k[4])) - y[1])) - (h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1]))/20
DG[2, 4] = -(h*a[1, 2]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3]))/20 + sin(t + h*c[1])*h*a[1, 2]*(h*(a[1, 1]*k[3] + a[1, 2]*k[6]) + y[3])
DG[2, 5] = 0
DG[2, 6] = -(sin(t + h*c[1])*h*a[1, 2]*(-(h*(a[1, 1]*k[1] + a[1, 2]*k[4])) - y[1])) - (h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1]))/20
DG[3, 1] = (h*a[1, 1]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]))/20 - 2*h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])
DG[3, 2] = (h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1]))/20 - h*a[1, 1]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])^2
DG[3, 3] = 1
DG[3, 4] = (h*a[1, 2]*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2]))/20 - 2*h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])*(h*(a[1, 1]*k[2] + a[1, 2]*k[5]) + y[2])
DG[3, 5] = (h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1]))/20 - h*a[1, 2]*(h*(a[1, 1]*k[1] + a[1, 2]*k[4]) + y[1])^2
DG[3, 6] = 0
DG[4, 1] = h*a[2, 1]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])
DG[4, 2] = -(sin(t + h*c[2])*h*a[2, 1]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])) + h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])
DG[4, 3] = -(sin(t + h*c[2])*h*a[2, 1]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])) + h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])
DG[4, 4] = 1 + h*a[2, 2]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])
DG[4, 5] = -(sin(t + h*c[2])*h*a[2, 2]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])) + h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])
DG[4, 6] = -(sin(t + h*c[2])*h*a[2, 2]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])) + h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])
DG[5, 1] = -(h*a[2, 1]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]))/20 + sin(t + h*c[2])*h*a[2, 1]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])
DG[5, 2] = 0
DG[5, 3] = -(sin(t + h*c[2])*h*a[2, 1]*(-(h*(a[2, 1]*k[1] + a[2, 2]*k[4])) - y[1])) - (h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1]))/20
DG[5, 4] = -(h*a[2, 2]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3]))/20 + sin(t + h*c[2])*h*a[2, 2]*(h*(a[2, 1]*k[3] + a[2, 2]*k[6]) + y[3])
DG[5, 5] = 1
DG[5, 6] = -(sin(t + h*c[2])*h*a[2, 2]*(-(h*(a[2, 1]*k[1] + a[2, 2]*k[4])) - y[1])) - (h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1]))/20
DG[6, 1] = (h*a[2, 1]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]))/20 - 2*h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])
DG[6, 2] = (h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1]))/20 - h*a[2, 1]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])^2
DG[6, 3] = 0
DG[6, 4] = (h*a[2, 2]*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2]))/20 - 2*h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])*(h*(a[2, 1]*k[2] + a[2, 2]*k[5]) + y[2])
DG[6, 5] = (h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1]))/20 - h*a[2, 2]*(h*(a[2, 1]*k[1] + a[2, 2]*k[4]) + y[1])^2
DG[6, 6] = 1
    return DG
end

@inbounds function irk(t,y,h)
    k1=fF(t+h*c[1],rk2(t,y,h*c[1]))
    k2=fF(t+h*c[2],rk2(t,y,h*c[2]))
    Ki=[k1;k2]
    lastone=false
    for i=1:10
        Gi=fG(t,y,h,Ki)
        DGi=fDG(t,y,h,Ki)
        DK=-DGi\Gi
        Ki+=DK
        if lastone
            break
        end
        if norm(DK)<=norm(Ki)*10e-10
            lastone=true
        end
    end
    k1=Ki[1:3]
    k2=Ki[4:6]
    return y+h*(b[1]*k1+b[2]*k2)
end

function solve(t0,y0,h,n)
    yn=copy(y0)
    for j=1:n
        tn=t0+(j-1)*h
        yn=irk(tn,yn,h)
    end
    return yn
end

function main(args)
    @printf("conjulia -- Energy Conserving IRK Scheme\n\n")
    t0=0; T=100
    y0=[1.0,1,1]
    if length(args)==3
        for j=1:3
            y0[j]=parse(Float64,args[j])
        end
    end
    println("y($t0)=$y0\n")
    einit=y0'*y0
    yold=zeros(3)
    @printf("%5s %7s %14s %14s %14s %14s\n",
        "p","n","h","E(100)","dE(100)","error")
    yn=copy(y0)
    n=64
    for p=6:16
        h=T/n
        yn=solve(t0,y0,h,n)
        if p>6
            @printf("%5d %7d %14.7e %14.7e %14.7e %14.7e\n",
            p,n,h,yn'*yn,yn'*yn-einit,norm(yold-yn))
        end
        flush(stdout)
        yold=copy(yn)
        n*=2
    end
    println("\ny($T)=$yn")
end

main(ARGS)

tsec=@elapsed main(ARGS)
println("\nTotal execution time ",tsec," seconds.")