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

General Usage
,
12

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