ode

ode.jl

using BenchmarkTools
using StaticArrays
function kernelrk(f, t, y, dt)
    
    k1 = f(t, y)
    k2 = f(t + 1/2*dt, y + 1/2*dt*k1)
    k3 = f(t + 3/4*dt, y + 3/4*dt*k2)
    y! = y + dt*(2/9*k1 + 1/3*k2 + 4/9*k3)
    # k4 = f(t + dt, y‘)
    # zn1 = y + dt*(7/24*k1 + 1/4*k2 + 1/3*k3 + 1/8*k4)
end

@fastmath function kernelrk2(f, t, y, dt)
    
    k1 = f(t, y)
    k2 = f(muladd(1/2, dt, t), muladd(1/2*dt, k1, y))
    k3 = f(muladd(3/4, dt, t), muladd(3/4*dt, k2, y))
    y! = muladd(dt, (2/9*k1 + 1/3*k2 + 4/9*k3), y)
    # k4 = f(t + dt, y‘)
    # zn1 = y + dt*(7/24*k1 + 1/4*k2 + 1/3*k3 + 1/8*k4)
end

function kernele(f, t, y, dt)   
    k1 = f(t, y)
    y! = y + dt*k1
    # k4 = f(t + dt, y‘)
    # zn1 = y + dt*(7/24*k1 + 1/4*k2 + 1/3*k3 + 1/8*k4)
end
f(t, x::Float64) = 1.2*x
const B = SMatrix{2,2,Float64}((1.1, 0.1, -0.1, 1.1))
f(t, x::SVector{2}) = B*x


function bench(kernelf, t0, y0, T, n = 100)
    dt = T/n
  
    t = t0
    y = y0

    for i in 1:n
        t = t + dt
        y = kernelf(f, t, y, dt)
    end
    y
end

T = 5.0


println("exp(-T*B), T = 5, B = -1.2")

println("Exp: ", exp(T*1.2))
println("Euler: ", bench(kernele, 0.0, 1.0, T, 1000))
println("RK: ", bench(kernelrk, 0.0, 1.0, T, 20))

println("Time Euler: ")
@btime  bench(kernele, 0.0, 1.0, T, 1000)
println("Time RK: ")
@btime  bench(kernelrk, 0.0, 1.0, T, 20)

y0 = ones(SVector{2, Float64})
println("Exp: ", expm(T*B)*y0)

println("Euler: ", bench(kernele, 0.0, y0, T, 1000))
println("RK: ", bench(kernelrk, 0.0, y0, T, 20))


println("Time expm: ")
@btime expm(T*B)*y0
println("Time Euler: ")
@btime  bench(kernele, 0.0, y0, T, 1000)
println("Time RK: ")
@btime  bench(kernelrk, 0.0, y0, T, 20)
@btime  bench(kernelrk, 0.0, y0, T, 1000)


using OrdinaryDiffEq

tspan = (0.0,T)
prob = ODEProblem(f, y0, tspan)
sol = solve(prob,Tsit5(), adaptive=false, dt=T/20 )