Comparing RK4 to Euler method

Question

I took a course in computational physics and I was given in an assigment to solve the equation of simple harmonic oscillator x_dot_dot = -x . I wrtoe the script in Julia for both methods

function RK4(x0 , F, T , dt) 

F1(f) = x -> f(x...)
F2(f) = x -> f((x .+ 0.5*dt*F1(f)(x))...)
F3(f) = x -> f((x .+ 0.5*dt*F2(f)(x))...)
F4(f) = x -> f((x .+ dt*F3(f)(x))...)

t = 0:dt:T
x = zeros(length(t), length(x0))
x[1,:] = x0
for n in 1:length(t) - 1
    act(F)  = map(f -> f(x[n,:]) , F)'
    x[n+1,:]= x[n,:] .+ (dt/6)*(act(F1.(F)) .+ 2*act(F2.(F)) + 2*act(F3.(F)) + act(F4.(F)))
end
return (x,t)
end

function euler(x0 , F, T , dt) 

F1(f) = x -> f(x...)

t = 0:dt:T
x = zeros(length(t), length(x0))
x[1,:] = x0
for n in 1:length(t)-1
    act(F)  = map(f -> f(x[n,:]) , F)'
    x[n+1,:]= x[n,:] .+ dt.* act(F1.(F))
end
return (x,t)
end

the code works and it solves the ODE for both methods by when i plot the solutions it seems that the RK4 isn't compensating the Euler method when trying them both with equal time steps. When i try to use with RK4 with say half time step than Euler, Euler always shows more accurte results and when i try to use a timestep of say dt = 0.01 for both methods the plots look almost identically the same.

The rest of the code :

using Plots
using LaTeXStrings
include("./methods.jl")


x0 = [1,0]
f1(x,v)= v
f2(x,v) = -x

plt = plot(xlabel = L"t")

xe , t = euler(x0 , [f1 f2] , 10*2π , 0.0005)
plot!(t, xe[:,1] , label = "Euler")
xr , t = RK4(x0 , [f1 f2] , 10*2π , 0.0005)
plot!(t, xr[:,1] , label = "RK4")

# Analytical solution
#plot!(t, cos.(t))

display(plt)

Thanks!

Answer 1

The reason RK4 is less accurate is because of a bug. The implementation in the question is really obfuscated. And even non-indented. There is definitely better and more efficient ways to implement this and there are certainly packages in Julia which do so.

A fixed RK4 can be:

function RK4(x0 , F, T , dt) 
    F1(f) = x -> f(x...)
    F2(f) = x -> f((x .+ 0.5*dt*[F1(g)(x) for g in F]')...)
    F3(f) = x -> f((x .+ 0.5*dt*[F2(g)(x) for g in F]')...)
    F4(f) = x -> f((x .+ dt*[F3(g)(x) for g in F]')...)

    t = 0:dt:T
    x = zeros(length(t), length(x0))
    x[1,:] = x0
    for n in 1:length(t) - 1
        act(F)  = map(f -> f(x[n,:]) , F)'
        x[n+1,:]= x[n,:] .+ (dt/6)*(act(F1.(F)) .+ 2*act(F2.(F)) + 2*act(F3.(F)) + act(F4.(F)))
    end
    return (x,t)
end

Note the definitions of F2, F3, F4 have changed. The original definition worked because the following broadcasting

julia> [1.2 1.2] .+ 0.5
1×2 Matrix{Float64}:
 1.7  1.7

worked on the F1(f)(x) output erroneously.