# Inconsistent typing making adding arrays and adding a vector to the end of an array difficult

I’m surprised that I actually found something scientific-computing related that is easier for me to do in JavaScript than it is in Julia. I’ve written functions that use the 4/5th order Runge-Kutta-Fehlberg method to solve a variety of different ODE systems. I’ve been trying to write a more general version of such a function (general in the sense that it should be able to solve ODE systems of any size), namely using the code:

``````using LinearAlgebra;
struct solObj
t
vars
end

function f(t, vars)
x = vars;
y = vars;
array = [x*cos(y), cos(x)];
return array;
end

function rkf45(f::Function, t0::Number, tf::Number, conds, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = Float64[t0];
vars = conds;
i = 1;
while ( t[i] < tf )
dt = min(dt, tf-t[i]);
K1 = dt*f(t[i], vars[i,:]);
K2 = dt*f(t[i] + dt/4, vars[i,:] + K1/4);
K3 = dt*f(t[i] + 3*dt/8, vars[i,:] + 3*K1/32 + 9*K2/32);
K4 = dt*f(t[i] + 12*dt/13, vars[i,:] + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197);
K5 = dt*f(t[i] + dt, vars[i,:] + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104);
K6 = dt*f(t[i] + dt/2, vars[i,:] - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40);
vars1 = vars[i,:] + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5;
vars2 = vars[i,:] + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55;
R = min(abs.(vars2-vars1)/dt...);
s = (epsilon/(2*R))^(1/4);
if ( R <= epsilon )
push!(t, t[i]+dt);
push!(vars, vars1);
i += 1;
end
dt *= s;
end
solution = solObj(t, vars);
return solution;
end
``````

(in this case I’ve written a function for `f` both to give you an example of what my input function could look like and also to help me test out my function). Now, K1, K2, K3, K4, K5 and K6 should ideally all be row vectors, as should vars[i,:], vars1 and vars2. vars is an array with each row corresponding to a different time value and the columns corresponding to different dependent variables (in this case x and y). Vars is used instead of x and y directly because this function is designed to be more generalized, such that it should be able to work on an ODE system of any size.

This code looked to me like it should work fine, but I get the error:

``````ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(1),), b has dims (Base.OneTo(2),), mismatch at 1")
Stacktrace:
 promote_shape at ./indices.jl:178 [inlined]
 promote_shape at ./indices.jl:169 [inlined]
 +(::Array{Float64,1}, ::Array{Float64,1}) at ./arraymath.jl:45
 rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Float64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45.jl:21
 top-level scope at REPL:1
 run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
``````

when I execute the command `rkf45(f, 0.0, 10.0, [0, 0], 1e-6, 0.1)`. Line 21 is the one on which K1 is defined. Despite this, given the error I’m getting I’m guessing it’s actually an issue with line 22 (where K2 is defined), as there’s actually an addition there. I’ve tried resolving the issue using the `reshape()` function, but feels like playing whac-a-mole as whenever one issue is temporarily resolved with it another pops up. For example the code:

``````using LinearAlgebra;
struct solObj
t
vars
end

function f(t::Number, vars)
x = vars;
y = vars;
array = [x*cos(y), cos(x)];
return array;
end

function rkf45(f::Function, t0::Number, tf::Number, conds::Array, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = [t0];
N = length(conds);
vars = reshape(conds, 1, N);
i = 1;
while ( t[i] < tf )
dt = min(dt, tf-t[i]);
K1 = reshape(dt*f(t[i], reshape(vars[i,:], 1, N)), 1, N);
K2 = reshape(dt*f(t[i] + dt/4, reshape(vars[i,:], 1, N) + K1/4), 1, N);
K3 = reshape(dt*f(t[i] + 3*dt/8, reshape(vars[i,:], 1, N) + 3*K1/32 + 9*K2/32), 1, N);
K4 = reshape(dt*f(t[i] + 12*dt/13, reshape(vars[i,:], 1, N) + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197), 1, N);
K5 = reshape(dt*f(t[i] + dt, reshape(vars[i,:], 1, N) + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104), 1, N);
K6 = reshape(dt*f(t[i] + dt/2, reshape(vars[i,:], 1, N) - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40), 1, N);
vars1 = reshape(vars[i,:], 1, N) + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5;
vars2 = reshape(vars[i,:], 1, N) + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55;
R = min(abs.(vars2-vars1)/dt...);
s = (epsilon/(2*R))^(1/4);
if ( R <= epsilon )
push!(t, t[i]+dt);
push!(vars, vars1);
i += 1;
end
dt *= s;
end
solution = solObj(t, vars);
return solution;
end
``````

when I run `rkf45(f, 0.0, 10.0, [0, 0], 1e-6, 0.1)` yields:

``````ERROR: MethodError: no method matching push!(::Array{Int64,2}, ::Array{Float64,2})

You might have used a 2d row vector where a 1d column vector was required.
Note the difference between 1d column vector [1,2,3] and 2d row vector [1 2 3].
You can convert to a column vector with the vec() function.
Closest candidates are:
push!(::Any, ::Any, ::Any) at abstractarray.jl:2252
push!(::Any, ::Any, ::Any, ::Any...) at abstractarray.jl:2253
push!(::Base.Nowhere, ::Any) at array.jl:1347
...
Stacktrace:
 rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Int64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45-reshape.jl:34
 top-level scope at REPL:1
 run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
``````

I have tried using `push!(vars, vec(vars1)')` instead like this error suggests may be best (as instead of a column vector, I want it to be a row vector) but that just gives the error:

``````ERROR: MethodError: no method matching push!(::Array{Int64,2}, ::Adjoint{Float64,Array{Float64,1}})

You might have used a 2d row vector where a 1d column vector was required.
Note the difference between 1d column vector [1,2,3] and 2d row vector [1 2 3].
You can convert to a column vector with the vec() function.
Closest candidates are:
push!(::Any, ::Any, ::Any) at abstractarray.jl:2252
push!(::Any, ::Any, ::Any, ::Any...) at abstractarray.jl:2253
push!(::Base.Nowhere, ::Any) at array.jl:1347
...
Stacktrace:
 rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Int64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45-reshape.jl:34
 top-level scope at REPL:1
 run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
``````

Any ideas of a way around this conundrum?

The easiest solution to this is to use Vectors not Row vectors. The only time to should be using Row vectors is when you are doing linear algebra like inner or outer projects. Also this is a textbook use case of StaticArrays

5 Likes

Also, are you sure you want to use `push!` and not `hcat` or `vcat`? I am having a hard time imagining how exactly pushing elements into a multidimensional array would work.

1 Like

Why are you reshaping vectors into 1xN matrices all over the place? That creates a horrible mess.

Aside from that, the error message explains the other problem, `push!` is for appending single elements to vectors, not for concatenating arrays.

2 Likes

The following edited code runs, but I have no idea if it is correct nor efficient:

``````function rkf45(f::Function, t0::Number, tf::Number, conds, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = Float64[t0]
vars = conds
i = 1
while ( t[i] < tf )
varsi =  vars[i,:]
dt = minimum((dt, tf-t[i]))
K1 = dt*f(t[i], varsi)
K2 = dt*f(t[i] + dt/4, varsi + K1/4)
K3 = dt*f(t[i] + 3*dt/8, varsi + 3*K1/32 + 9*K2/32)
K4 = dt*f(t[i] + 12*dt/13, varsi + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197)
K5 = dt*f(t[i] + dt, varsi + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104)
K6 = dt*f(t[i] + dt/2, varsi - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40)
vars1 = varsi + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5
vars2 = varsi + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55
R = minimum(abs.(vars2-vars1)/dt)
s = (epsilon/(2*R))^(1/4)   # can R = 0 depending on initial conditions?
if ( R <= epsilon )
t = vcat(t, t[i]+dt)
vars = vcat(vars, vars1')
i += 1
end
dt *= s
end
solution = solObj(t, vars)
return solution;
end
``````

To call the function, the initial conditions should be a row vector :

``````sol = rkf45(f, 0.0, 10, [0.5 0.5], 1e-6, 0.1)
``````
1 Like

Thanks Rafael! Your code works perfectly. My example f wasn’t exactly the best, I tried with the new f (coming from the problem of the simple pendulum):

``````function f(t, vars)
g = 9.81
l = 1.0
x = vars
dx = vars
return [dx, -g/l*cos(x)]
end
``````

and it worked perfectly. In reality R should not be zero, unless the ODE is very simple as R measures the error in our approximation of vars.

I spoke too soon. It works for the 2D case, but for a system of 3 ODE it fails. Namely the code:

``````using PyPlot;
include("RKF45.jl");

function f(params, t, vars)
sigma = params.sigma;
beta  = params.beta;
rho   = params.rho;
x     = vars;
y     = vars;
z     = vars;
dx    = sigma*(y-x);
dy    = x*(rho-z) - y;
dz    = x*y - beta*z;
return [dx, dy, dz];
end

struct paramObj
sigma
beta
rho
end

# Problem parameters
sigma = 10;
beta = 8/3;
rho = 28;
params = paramObj(sigma, beta, rho);

t0 = 0.0;
tf = 60.0;
x0 = 1.0;
y0 = 1.0;
z0 = 1.0;
conds = [x0 y0 z0];
epsilon = 1e-9;
dtInitial = 0.1;

solution = RKF45(f, params, t0, tf, conds, epsilon, dtInitial);
x = solution.vars[:,1];
y = solution.vars[:,2];
z = solution.vars[:,3];
PyPlot.figure(1)
PyPlot.clf();
PyPlot.plot3D(x, y, z);
``````

(where RKF45.jl is where I stored rkf45, just renamed to RKF45) gives the error:

``````ERROR: LoadError: MethodError: no method matching RKF45(::typeof(f), ::paramObj, ::Float64, ::Float64, ::Array{Float64,2}, ::Float64, ::Float64)

You might have used a 2d row vector where a 1d column vector was required.
Note the difference between 1d column vector [1,2,3] and 2d row vector [1 2 3].
You can convert to a column vector with the vec() function.
Closest candidates are:
RKF45(::Function, ::Any, ::Number, ::Number, ::Array{T,1} where T, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45.jl:34
Stacktrace:
 top-level scope at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Lorenz.jl:38
 include(::String) at ./client.jl:457
 top-level scope at REPL:1
 run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
in expression starting at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Lorenz.jl:38
``````

Line 34 of RKF45.jl is:

``````function RKF45(f::Function, params, t0::Number, tf::Number, conds::Vector, epsilon::Float64, dtInitial::Float64)
``````

EDIT: Oops, sorry! I just discovered my error. I added `::Vector` after conds in the function definition and that caused the error.

That code does seem orders of magnitude less efficient than my JavaScript version of RKF45 (based on how fast the above 3D example executes versus how quickly this webpage solves the same problem). But, I will admit I don’t get why. It looks like it should be as efficient, if not more efficient than my JS code. Here’s my JS version of RKF45 if it helps: https://github.com/fusion809/fusion809.github.io/blob/dev/_libs/common/RKF45.js. The code of that webpage is here: https://github.com/fusion809/fusion809.github.io/tree/dev/Lorenz. Any one have any ideas?

Found the problem. Even though 0 should never be obtained for R, it sometimes is for this code but not for my JS code (oddly enough). So, guarding against division by zero fixes the issue.

Oh and I figured out why 0 is sometimes obtained for R. I was meant to be maxing R (i.e. using `R = maximum(abs.(vars2-vars1)/dt);` not `R = minimum((abs.(vars2-vars1)/dt);`. Thanks for your patience and help everyone, you’ve all been immensely helpful, one reason why I think the Julia community is one of the most helpful open-source software communities I’ve ever come across.

1 Like

To collects vars, you could use a Vector of Vectors:

``````vars = [conds]
``````

and later

``````push!(vars, vars1)
``````

If you really want a matrix later you can do

``````reduce(hcat, vars)
``````

Using SVector(0.0, 0.0) (from StaticArrays) as the initial value and modifying the rest of the code to return SVectors when given SVectors, should help getting rid of a lot of small allocations, that I assume you are getting.

3 Likes

Thanks for the alternative approach. Tried it, and I just get errors, so I’m guessing there’s something I’m missing. Here’s the modified RKF45.jl file:

``````# Solution object
struct solObj
t
vars
end

function RKF45(f::Function, params, t0::Number, tf::Number, conds, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = Float64[t0];
vars = [conds];
i = 1;
while ( t[i] < tf )
varsi =  vars[i,:];
dt = minimum((dt, tf-t[i]));

# Interpolants
K1 = dt*f(params, t[i], varsi);
K2 = dt*f(params, t[i] + dt/4, varsi + K1/4);
K3 = dt*f(params, t[i] + 3*dt/8, varsi + 3*K1/32 + 9*K2/32);
K4 = dt*f(params, t[i] + 12*dt/13, varsi + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197);
K5 = dt*f(params, t[i] + dt, varsi + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104);
K6 = dt*f(params, t[i] + dt/2, varsi - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40);

# 4th/5th order approximations
vars1 = varsi + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5;
vars2 = varsi + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55;

# Adjust step size if needed
R = maximum(abs.(vars2-vars1)/dt);
s = (epsilon/(2*R))^(1/4);
if ( R <= epsilon )
t = vcat(t, t[i]+dt);
vars = push!(vars, vars1);
i += 1;
end
dt *= s;
end
vars = reduce(hcat, vars);
solution = solObj(t, vars);
return solution;
end
``````

(I did omit a docstring that was messing with the code formatting of Markdown). My code that executes this function is:

``````using PyPlot;
using StaticArrays;
include("RKF45.jl");

function f(params, t, vars)
sigma = params.sigma;
beta  = params.beta;
rho   = params.rho;
x     = vars;
y     = vars;
z     = vars;
dx    = sigma*(y-x);
dy    = x*(rho-z) - y;
dz    = x*y - beta*z;
return SVector(dx, dy, dz);
end

struct paramObj
sigma
beta
rho
end

# Problem parameters
sigma = 10;
beta = 8/3;
rho = 28;
params = paramObj(sigma, beta, rho);

t0 = 0.0;
tf = 60.0;
x0 = 1.0;
y0 = 1.0;
z0 = 1.0;
conds = SVector(x0, y0, z0);
epsilon = 1e-9;
dtInitial = 0.1;

solution = RKF45(f, params, t0, tf, conds, epsilon, dtInitial);
x = solution.vars[:,1];
y = solution.vars[:,2];
z = solution.vars[:,3];
PyPlot.figure(1)
PyPlot.clf();
PyPlot.plot3D(x, y, z);
``````

and it returns the error:

``````ERROR: LoadError: BoundsError: attempt to access 1-element Array{SArray{Tuple{3},Float64,1,3},1} at index 
Stacktrace:
 getindex at ./array.jl:809 [inlined]
 f(::paramObj, ::Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:12
 RKF45(::typeof(f), ::paramObj, ::Float64, ::Float64, ::SArray{Tuple{3},Float64,1,3}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/RKF45.jl:44
 top-level scope at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:41
 include(::String) at ./client.jl:457
 top-level scope at REPL:1
 run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
in expression starting at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:41
``````

I did try without SVector, namely using the code:

``````using PyPlot;
include("RKF45.jl");

function f(params, t, vars)
sigma = params.sigma;
beta  = params.beta;
rho   = params.rho;
x     = vars;
y     = vars;
z     = vars;
dx    = sigma*(y-x);
dy    = x*(rho-z) - y;
dz    = x*y - beta*z;
return [dx, dy, dz];
end

struct paramObj
sigma
beta
rho
end

# Problem parameters
sigma = 10;
beta = 8/3;
rho = 28;
params = paramObj(sigma, beta, rho);

t0 = 0.0;
tf = 60.0;
x0 = 1.0;
y0 = 1.0;
z0 = 1.0;
conds = [x0 y0 z0];
epsilon = 1e-9;
dtInitial = 0.1;

solution = RKF45(f, params, t0, tf, conds, epsilon, dtInitial);
x = solution.vars[:,1];
y = solution.vars[:,2];
z = solution.vars[:,3];
PyPlot.figure(1)
PyPlot.clf();
PyPlot.plot3D(x, y, z);
``````

and it gives the error:

``````ERROR: LoadError: BoundsError: attempt to access 1-element Array{Array{Float64,2},1} at index 
Stacktrace:
 getindex at ./array.jl:809 [inlined]
 f(::paramObj, ::Float64, ::Array{Array{Float64,2},1}) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:10
 RKF45(::typeof(f), ::paramObj, ::Float64, ::Float64, ::Array{Float64,2}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/RKF45.jl:44
 top-level scope at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:39
 include(::String) at ./client.jl:457
 top-level scope at REPL:1
 run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
in expression starting at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:39
``````

Two general performance tips:

When you define structs, the fields should have concrete types, if you don’t specify they will be of type `Any`, which is very bad for performance.

Also, prefer to work (slicing and iteration for example) along array columns not rows, since Julia arrays are column major.

4 Likes

Use `vars[i]` instead of `vars[i,:]` when vars is a vector. For the SVector approach you need to rewrite `f` and maybe something else.

4 Likes

Trying to follow the advice from the experts (@mikkoku, @DNF, …) , the edited code below seems to run faster. What do you think?

``````using LinearAlgebra, Plots;

struct solObj
t::Array{Float64}
vars::Array{Float64}
end

struct paramObj
sigma::Float64
beta::Float64
rho::Float64
end

function f(params, t, vars)
sigma = params.sigma;
beta  = params.beta;
rho   = params.rho;
x     = vars;
y     = vars;
z     = vars;
dx    = sigma*(y-x);
dy    = x*(rho-z) - y;
dz    = x*y - beta*z;
return [dx, dy, dz];
end

function rkf45(f::Function, params, t0::Float64, tf::Float64, conds::Array{Float64}, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = Float64[t0]
vars = [conds]
i = 1
ti = t0
while ( ti < tf )
varsi =  vars[i]
dt = minimum((dt, tf-ti))
K1 = dt*f(params, ti, varsi)
K2 = dt*f(params, ti + dt/4, varsi + K1/4)
K3 = dt*f(params, ti + 3*dt/8, varsi + 3*K1/32 + 9*K2/32)
K4 = dt*f(params, ti + 12*dt/13, varsi + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197)
K5 = dt*f(params, ti + dt, varsi + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104)
K6 = dt*f(params, ti + dt/2, varsi - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40)
vars1 = varsi + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5
vars2 = varsi + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55
R = maximum(abs.(vars2 - vars1))/dt
s = (epsilon/(2*R))^0.25
if (R <= epsilon)
push!(t, ti+dt)
push!(vars, vars1)
i += 1
ti = t[i]
end
dt *= s
end
vars = reduce(hcat, vars)
return solObj(t, vars);
end

# Problem parameters
sigma = 10;
beta = 8/3;
rho = 28;
params = paramObj(sigma, beta, rho);
t0 = 0.0;
tf = 60.0;
x0 = 1.0;
y0 = 1.0;
z0 = 1.0;
conds = [x0;y0;z0];
epsilon = 1e-9;
dtInitial = 0.1;

sol = rkf45(f, params, t0, tf, conds, epsilon, dtInitial)

Plots.plot(sol.vars[1,:],sol.vars[2,:])
``````
1 Like

`Array{Float64}` is still an abstract type, since you didn’t specify dimensionality. Try `Array{Float64, N}` for whichever `N` is appropriate.

4 Likes

I think the bigger performance benefit would come from using SVector. The `f` function part that @fusion809 already tried looks good to me. (I read it too hastily first time.)

2 Likes

Thanks everyone. I really appreciate your help. I now have this code, after grabbing Rafael’s code and making the suggested improvements:

``````# Solution object
struct solObj
t::Array{Float64,1}
vars::Array{Float64,2}
end

function RKF45(f::Function, params, t0::Float64, tf::Float64, conds::Array{Float64}, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = Float64[t0]
vars = [conds]
i = 1
ti = t0
while ( ti < tf )
varsi =  vars[i]
dt = minimum((dt, tf-ti))
K1 = dt*f(params, ti, varsi)
K2 = dt*f(params, ti + dt/4, varsi + K1/4)
K3 = dt*f(params, ti + 3*dt/8, varsi + 3*K1/32 + 9*K2/32)
K4 = dt*f(params, ti + 12*dt/13, varsi + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197)
K5 = dt*f(params, ti + dt, varsi + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104)
K6 = dt*f(params, ti + dt/2, varsi - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40)
vars1 = varsi + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5
vars2 = varsi + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55
R = maximum(abs.(vars2 - vars1))/dt
s = (epsilon/(2*R))^0.25
if (R <= epsilon)
push!(t, ti+dt)
push!(vars, vars1)
i += 1
ti = t[i]
end
dt *= s
end
vars = transpose(reduce(hcat, vars));
return solObj(t, vars);
end
``````

and it executes a lot faster (e.g. for the Lorenz example it took about 12-14 seconds to execute before, now it takes as little as 0.5 seconds, it seems a little faster than my JS code now). I transposed vars because it feels more intuitive to work with the columns being separate variable values (x, y, z).

I haven’t implemented SVector as I just find it confusing, I tried it before and it failed and I haven’t the foggiest how to fix the error I got (which I reported earlier in this post). Besides this code is fast enough for my purposes.

We can get another 10x performance increase by using StaticArrays as recommended by @mikkoku and @Sukera. The edited code below will run in ~14 ms (@btime). Note that for this problem I did not observe any time improvement by defining in struct, t or vars as ::Array{Float64,N} versus ::Array{Float64} as recommended by @DNF, but do not understand the implications of this in general case.

``````using LinearAlgebra, Plots, StaticArrays
using BenchmarkTools

struct solObj
t::Array{Float64}
vars::Array{Float64}
end

struct paramObj
sigma::Float64
beta::Float64
rho::Float64
end

function f(params, t, vars)
sigma, beta, rho = params.sigma, params.beta, params.rho;
x, y, z = vars, vars, vars;
dx    = sigma*(y-x);
dy    = x*(rho-z) - y;
dz    = x*y - beta*z;
return  @SVector [dx, dy, dz];
end

function rkf45(f::Function, params, t0::Float64, tf::Float64, conds, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = Float64[t0]
vars = [conds]
i = 1
ti = t0
while ( ti < tf )
varsi =  vars[i]
dt = minimum((dt, tf-ti))
K1 = dt*f(params, ti, varsi)
K2 = dt*f(params, ti + dt/4, varsi + K1/4)
K3 = dt*f(params, ti + 3*dt/8, varsi + 3*K1/32 + 9*K2/32)
K4 = dt*f(params, ti + 12*dt/13, varsi + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197)
K5 = dt*f(params, ti + dt, varsi + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104)
K6 = dt*f(params, ti + dt/2, varsi - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40)
vars1 = varsi + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5
vars2 = varsi + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55
R = maximum(abs.(vars2 - vars1))/dt
s = (epsilon/(2*R))^0.25
if (R <= epsilon)
push!(t, ti+dt)
push!(vars, vars1)
i += 1
ti = t[i]
end
dt *= s
end
vars = transpose(reduce(hcat, vars))   # Optional transpose
return solObj(t, vars);
end

# Problem parameters
sigma = 10;
beta = 8/3;
rho = 28;
params = paramObj(sigma, beta, rho);
t0, tf = 0.0, 60.0;
x0, y0, z0 = 1.0, 1.0, 1.0;
conds = @SVector [x0,y0,z0];
epsilon = 1e-9;
dtInitial = 0.1;

sol = rkf45(f, params, t0, tf, conds, epsilon, dtInitial)
# @btime  rkf45(f, params, t0, tf, conds, epsilon, dtInitial)

Plots.plot(sol.vars[:,1],sol.vars[:,2])
``````
1 Like