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!
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.