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

fusion809 · October 2, 2020, 8:15pm

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[1];
    y = vars[2];
    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:
 [1] promote_shape at ./indices.jl:178 [inlined]
 [2] promote_shape at ./indices.jl:169 [inlined]
 [3] +(::Array{Float64,1}, ::Array{Float64,1}) at ./arraymath.jl:45
 [4] rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Float64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45.jl:21
 [5] top-level scope at REPL[114]:1
 [6] 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[1];
    y = vars[2];
    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:
 [1] rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Int64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45-reshape.jl:34
 [2] top-level scope at REPL[6]:1
 [3] 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:
 [1] rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Int64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45-reshape.jl:34
 [2] top-level scope at REPL[3]:1
 [3] 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?

Oscar_Smith · October 2, 2020, 8:49pm

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

Henrique_Becker · October 2, 2020, 10:47pm

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.

DNF · October 2, 2020, 11:24pm

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.

rafael.guerra · October 2, 2020, 11:43pm

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)

fusion809 · October 3, 2020, 6:45am

rafael.guerra:

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

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[1]
   dx = vars[2]
  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.

fusion809 · October 3, 2020, 7:55am

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[1];
    y     = vars[2];
    z     = vars[3];
    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:
 [1] top-level scope at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Lorenz.jl:38
 [2] include(::String) at ./client.jl:457
 [3] top-level scope at REPL[1]:1
 [4] 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.

fusion809 · October 3, 2020, 8:10am

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?

fusion809 · October 3, 2020, 9:25am

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.

fusion809 · October 3, 2020, 9:29am

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.

mikkoku · October 3, 2020, 9:36am

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.

fusion809 · October 3, 2020, 9:54am

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[1];
    y     = vars[2];
    z     = vars[3];
    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 [2]
Stacktrace:
 [1] getindex at ./array.jl:809 [inlined]
 [2] f(::paramObj, ::Float64, ::Array{SArray{Tuple{3},Float64,1,3},1}) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:12
 [3] 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
 [4] top-level scope at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:41
 [5] include(::String) at ./client.jl:457
 [6] top-level scope at REPL[1]:1
 [7] 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

fusion809 · October 3, 2020, 9:56am

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[1];
    y     = vars[2];
    z     = vars[3];
    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 [2]
Stacktrace:
 [1] getindex at ./array.jl:809 [inlined]
 [2] f(::paramObj, ::Float64, ::Array{Array{Float64,2},1}) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:10
 [3] RKF45(::typeof(f), ::paramObj, ::Float64, ::Float64, ::Array{Float64,2}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/RKF45.jl:44
 [4] top-level scope at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/Generalized/Lorenz.jl:39
 [5] include(::String) at ./client.jl:457
 [6] top-level scope at REPL[1]:1
 [7] 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

DNF · October 3, 2020, 9:58am

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.

mikkoku · October 3, 2020, 10:03am

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.

rafael.guerra · October 3, 2020, 11:06pm

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[1];
    y     = vars[2];
    z     = vars[3];
    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,:])

DNF · October 4, 2020, 5:14am

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

mikkoku · October 4, 2020, 6:02am

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

fusion809 · October 4, 2020, 9:58am

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.

rafael.guerra · October 4, 2020, 11:29am

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[1], vars[2], vars[3];
    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])