# Convergence of RK4 solver
**URL:** https://discourse.julialang.org/t/convergence-of-rk4-solver/113770
**Category:** Numerics
**Created:** [May 3, 2024, 1:22am UTC](https://discourse.julialang.org/t/convergence-of-rk4-solver/113770 "2024-05-03T01:22:17Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@Pranav\_Vinod](https://discourse.julialang.org/u/Pranav_Vinod)
#### Post date: [May 3, 2024, 1:22am UTC](https://discourse.julialang.org/t/convergence-of-rk4-solver/113770/1 "2024-05-03T01:22:17Z")
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.
```julia
# 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
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 :
```julia
[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.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 3, 2024, 1:33am UTC](https://discourse.julialang.org/t/convergence-of-rk4-solver/113770/2 "2024-05-03T01:33:06Z")
> [@Pranav\_Vinod](#):
>
> I’m not sure what I’m missing in my code.
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:
```julia
velocity = 0.05f0
dt = Float32[100,10,1,0.1,0.01]
```
then I get results identical to DifferentialEquations.jl:
```julia
julia> error_custom_rk4
5-element Vector{Any}:
0.0f0
0.0f0
0.048553467f0
0.4050598f0
2.3095856f0
```
---
### Author:  [@Pranav\_Vinod](https://discourse.julialang.org/u/Pranav_Vinod)
#### Post date: [May 3, 2024, 2:11am UTC](https://discourse.julialang.org/t/convergence-of-rk4-solver/113770/3 "2024-05-03T02:11:48Z")
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?
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [May 3, 2024, 2:27am UTC](https://discourse.julialang.org/t/convergence-of-rk4-solver/113770/4 "2024-05-03T02:27:08Z")
Once you’re running up against floating point precision, smaller dt just means more occurrences of rounding error.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [May 3, 2024, 12:39pm UTC](https://discourse.julialang.org/t/convergence-of-rk4-solver/113770/5 "2024-05-03T12:39:55Z")
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
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
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.)
---
### Author:  [@Pranav\_Vinod](https://discourse.julialang.org/u/Pranav_Vinod)
#### Post date: [May 3, 2024, 2:40pm UTC](https://discourse.julialang.org/t/convergence-of-rk4-solver/113770/6 "2024-05-03T14:40:32Z")
Thank you to both of you! That makes sense and is really helpful!