Convergence of RK4 solver

Specific Domains Numerics
1

Hi,
I am working on solving some differential equations and noticed some discrepancies in the convergence of the RK4 solver. I tested the error by finding the difference between the solver solution and analytical results and computing the absolute maximum of the error vector. To check the method itself, I coded up the RK4 method from scratch and found that it seemed to converge on smaller values of dt.
Here is a MWE of the method I am using on a much simpler problem.

# ODE is \dot{u}   = velocity
# solution is u(t) = velocity*t

# using DifferentialEquations
# using DataFrames
# using Plots

u0    = Float32[0.0]    # initial conditions
tspan = (0.0f0, 1.0f4)  # time range
model_params = []       # RHS is not parameterized
velocity = 0.05     # RHS value



function RHS(u, model_params, t)
    return [velocity]
end

## Custom RK4 solver

function rk4_step(f, u, p, t, dt)
    k1 = f(u, p, t)
    k2 = f(u + dt/2 * k1, p, t + dt/2)
    k3 = f(u + dt/2 * k2, p, t + dt/2)
    k4 = f(u + dt * k3, p, t + dt)
    return u + dt/6 * (k1 + 2*k2 + 2*k3 + k4)
  end
  
function rk4_solve(f, u0, tspan, p, dt)
    t0, tf = tspan
    t = t0
    u = u0
    ts = [t]
    us = [u]
    while t < tf
        u = rk4_step(f, u, p, t, dt)
        t += dt
        push!(ts, t)
        push!(us, u)
    end
    return ts, us
end

dt = [100,10,1,0.1,0.01]
error_julia_rk4 = []
error_custom_rk4 = []
for dts in dt

    ## Solution and error from custom solver
    custom_tsteps ,custom_sol = rk4_solve(RHS, u0, tspan, model_params, dts)
    sol_phi = [u[1] for u in custom_sol]
    x = sol_phi .- custom_tsteps*velocity
    error_custom = maximum(abs.(x))
    
    ## Solution and error from julia solver
    julia_sol = solve(prob, RK4(),dt = dts, adaptive=false)
    y = julia_sol[1,:] .- julia_sol.t*velocity
    error_julia = maximum(abs.(y))
    
    push!(error_julia_rk4, error_julia)
    push!(error_custom_rk4, error_custom)
end

The output I get is the following :

julia> error_julia_rk4
5-element Vector{Any}:
 0.0
 0.0
 0.048553466796875
 0.40500488281247726
 2.3095886230468636

julia> error_custom_rk4
5-element Vector{Any}:
 0.0
 0.0
 1.2207031261368684e-5
 3.0517578125e-5
 3.662109378410605e-5

I’m not sure what I’m missing in my code. I would appreciate any help and suggestions.

Here are the package versions I’m using :

[0c46a032] DifferentialEquations v7.6.0
[91a5bcdd] Plots v1.39.0
[a93c6f00] DataFrames v1.3.6

Please let me know if I can provide more information.

2

Your custom solver is doing some of the calculations in double precision (Float64) because your dt and velocity variables are double precision, whereas DifferentialEquations.jl is doing all the calculations in single precision (Float32) because that is the precision of tspan and u0.

If I fix your code to use single precision too, by setting:

velocity = 0.05f0
dt = Float32[100,10,1,0.1,0.01]

then I get results identical to DifferentialEquations.jl:

julia> error_custom_rk4
5-element Vector{Any}:
 0.0f0
 0.0f0
 0.048553467f0
 0.4050598f0
 2.3095856f0
3 Likes
3

Thank you! That fixes the discrepancy between the custom solver and the once used by DifferentialEquations.jl .
Any ideas as to why they are not converging with the correct order with smaller dt? In that, the error seems to be growing as dt becomes smaller, but wouldn’t we expect the opposite?

4

Once you’re running up against floating point precision, smaller dt just means more occurrences of rounding error.

2 Likes
5

Right, in fact in this example your errors are all rounding error.

Your ODE is dx/dt = v where v is a constant — RK4, or even a simple Euler scheme x_k = x_{k-1} + v \Delta t (equivalent to RK4 in this case!), is exact for any \Delta t in the absence of roundoff error.

So, your error is entirely due to the roundoff error of the running sum, which accumulates linearly with the number of steps in this case because v and \Delta t are not exactly representable.

In fact, if I make v and \Delta t exactly representable (= integers \times powers of 2):

julia> velocity = 2.0f0^-5
0.03125f0

julia> dt = 2.0f0 .^ [6, 3, 0, -3, -6]
5-element Vector{Float32}:
 64.0
  8.0
  1.0
  0.125
  0.015625

then I get an error of zero independent of \Delta t:

julia> error_custom_rk4
5-element Vector{Any}:
 0.0f0
 0.0f0
 0.0f0
 0.0f0
 0.0f0

(This won’t be the case for a more complicated ODE with irrational solutions, however.)

2 Likes
6

Thank you to both of you! That makes sense and is really helpful!