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

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.

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[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

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

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[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,:])
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[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])
1 Like

I am in awe of how immensely helpful you’ve been, Rafael, thank you! While your code cuts down execution time on my system, it doesn’t quite go down to 0.03s, but ~0.18s is still pretty good, significantly faster than my JS code (which takes about 0.77s).

3 Likes

As RKF45 is meant to solve ODE systems of any size, conds should probably be generalized (i.e. not have its size fixed at 3 elements long). I’ve tried changing to conds::SVector{Float64}, but that gives an error (namely:

MethodError: no method matching rkf45(::typeof(Lorenz), ::NamedTuple{(:sigma, :beta, :rho),Tuple{Float64,Float64,Float64}}, ::Float64, ::Float64, ::SArray{Tuple{3},Float64,1,3}, ::Float64, ::Float64)

). Any ideas on how to resolve this problem?

StaticArrays is a package that makes the size of the array available to the compiler for optimizations. SVector{3, UInt} means “This vector holds 3 UInt”. There’s no need to restrict the argument to that though (it doesn’t allow for more optimizations, the performance will depend on the types of the objects passed in, not on the types of the assertions) - if you remove the type assertions from the function arguments and simply pass in objects of the type you want (SVector in this case), the function will be compiled optimized for the types that are passed in.

So if you remove the type assertion and simply make conds = @Svector [x0, y0, z0, w0], it should work fine for the 4 dimensional case.

3 Likes

Ah, so for SVectors specifying the type in function definitions does not further optimize the performance of the function?

Asserting function argument types doesn’t optimize anything, ever. It only restricts what types are allowed to be used to call the function at all and is used to disambiguate between different methods of the same function with the same arity (=number of arguments).

Note that there is an upper limit to the usefulness of SVectors - at some point, the large known size can increase compile times dramatically.

6 Likes

@time includes compilation time on the first run. Try using @btime from BenchmarkTools.jl for somewhat similar functionality instead.

2 Likes