I don't understand an error message I get that's related to function ode45.

The error message I get is:

Arrays have incompatible sizes for this operation.

Error in ode45 (line 299)

y5 = y + h .* (b41.*f1 + b42.*f2 + b43.*f3 + b44.*f4 );

I've checked my initial condition vector is the correct size. I'm not sure what this error indicated might be wrong in my code.

Answer 1

[Steven Lord]

More likely than not, your ODE function is returning a vector of a different size for one set of inputs than it does for other inputs. As stated on the ode45 documentation page, your function must always return a column vector and that column vector must be the same size as the second input to your function.

"The function dydt = odefun(t,y), for a scalar t and a column vector y, must return a column vector dydt of data type single or double that corresponds to f(t,y)."

At the end of your ODE function, you could check that dydt is the same size as y and enter debug mode if it is not. Add the following three lines at the end of your ODE function, right before you return from that function. [I'm using dydt as the name of the output from your ODE function and y as the second input; if you're using different names, make the appropriate substitutions below.]

if ~isequal(size(dydt), size(y))
    keyboard
end

Then run your code. If I'm right, MATLAB will stop on the line with the keyboard command. When it does, look at what your ODE function is doing for the inputs it received and determine why dydt is not the same size as y. You could save the input arguments to a MAT-file, exit the ode45 call, then call your ODE function directly with those inputs to assist in debugging if the problem only occurs after a large number of ODE time steps (that you don't want to have to execute each time.)

If for whatever reason you can't edit your function (you must use it exactly as it was provided to you by your professor or boss, as one example) you could set a conditional breakpoint on the last line of your ODE function with ~isequal(size(dydt), size(y)) as the condition.

[Sierra Cagle]

Thank you - this was the trouble shooting tip I needed to find the bug in my code!

Answer 2

[Torsten]

If you scroll through the code, you will find the line in ode45. Since the ode integrator works for allowed inputs, check the sizes of zin and zdot.

Running your code, you will notice that zdot changes size from 39x1 to 42x1 in one call to ode45. This is not allowed.

type ode45

function varargout = ode45(ode,tspan,y0,options,varargin)
%ODE45  Solve non-stiff differential equations, medium order method.
%   [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0) integrates the system of
%   differential equations y' = f(t,y) from time TSPAN(1) to TSPAN(end)
%   with initial conditions Y0. Each row in the solution array YOUT
%   corresponds to a time in the column vector TOUT. 
%     * ODEFUN is a function handle. For a scalar T and a vector Y,
%       ODEFUN(T,Y) must return a column vector corresponding to f(t,y).
%     * TSPAN is a two-element vector [T0 TFINAL] or a vector with
%       several time points [T0 T1 ... TFINAL]. If you specify more than
%       two time points, ODE45 returns interpolated solutions at the
%       requested times.
%     * YO is a column vector of initial conditions, one for each equation.
%
%   [TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) specifies integration
%   option values in the fields of a structure, OPTIONS. Create the
%   options structure with odeset.
%
%   [TOUT,YOUT,TE,YE,IE] = ODE45(ODEFUN,TSPAN,Y0,OPTIONS) produces
%   additional outputs for events. An event occurs when a specified function
%   of T and Y is equal to zero. See ODE Event Location for details.
%
%   SOL = ODE45(...) returns a solution structure instead of numeric
%   vectors. Use SOL as an input to DEVAL to evaluate the solution at
%   specific points. Use it as an input to ODEXTEND to extend the
%   integration interval.
%
%   ODE45 can solve problems M(t,y)*y' = f(t,y) with mass matrix M that is
%   nonsingular. Use ODESET to set the 'Mass' property to a function handle
%   or the value of the mass matrix. ODE15S and ODE23T can solve problems
%   with singular mass matrices.
%
%   ODE23, ODE45, ODE78, and ODE89 are all single-step solvers that use
%   explicit Runge-Kutta formulas of different orders to estimate the error
%   in each step.
%     * ODE45 is for general use.
%     * ODE23 is useful for moderately stiff problems.
%     * ODE78 and ODE89 may be more efficient than ODE45 on non-stiff problems
%       that are smooth except possibly for a few isolated discontinuities.
%     * ODE89 may be more efficient than ODE78 on very smooth problems, when 
%       integrating over long time intervals, or when tolerances are tight.
%
%   Example
%         [t,y]=ode45(@vdp1,[0 20],[2 0]);   
%         plot(t,y(:,1));
%     solves the system y' = vdp1(t,y), using the default relative error
%     tolerance 1e-3 and the default absolute tolerance of 1e-6 for each
%     component, and plots the first component of the solution. 
%   
%   Class support for inputs TSPAN, Y0, and the result of ODEFUN(T,Y):
%     float: double, single
%
%   See also ODE23, ODE78, ODE89, ODE113, ODE15S, ODE23S, ODE23T, ODE23TB,
%            ODE15I, ODESET, ODEPLOT, ODEPHAS2, ODEPHAS3, ODEPRINT, DEVAL,
%            ODEEXAMPLES, FUNCTION_HANDLE.

%   ODE45 is an implementation of the explicit Runge-Kutta (4,5) pair of
%   Dormand and Prince called variously RK5(4)7FM, DOPRI5, DP(4,5) and DP54.
%   It uses a "free" interpolant of order 4 communicated privately by
%   Dormand and Prince.  Local extrapolation is done.

%   Details are to be found in The MATLAB ODE Suite, L. F. Shampine and
%   M. W. Reichelt, SIAM Journal on Scientific Computing, 18-1, 1997.

%   Mark W. Reichelt and Lawrence F. Shampine, 6-14-94
%   Copyright 1984-2023 The MathWorks, Inc.

solver_name = 'ode45';

% Check inputs
if nargin < 4
    options = [];
    if nargin < 3
        y0 = [];
        if nargin < 2
            tspan = [];
        end
    end
end

% Stats
nsteps  = 0;
nfailed = 0;
nfevals = 0;

[ode, odeIsFuncHandle, odeTreatAsMFile] = packageAsFuncHandle(ode);

% Output
output_sol = odeIsFuncHandle && nargout == 1; % sol = odeXX(...)
output_ty  = ~output_sol && nargout > 0;      % [t,y,...] = odeXX(...)
% There might be no output requested...

sol = []; f3d = [];
if output_sol
    sol.solver = solver_name;
    sol.extdata.odefun = ode;
    sol.extdata.options = options;
    sol.extdata.varargin = varargin;
end

% Handle solver arguments
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, ...
    options, threshold, rtol, normcontrol, normy, hmax, htry, htspan, dataType] = ...
    odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);
nfevals = nfevals + 1;

% Handle the output
if nargout > 0
    outputFcn = odeget(options,'OutputFcn',[]);
else
    outputFcn = odeget(options,'OutputFcn',@odeplot);
end
outputArgs = {};
if isempty(outputFcn)
    haveOutputFcn = false;
else
    haveOutputFcn = true;
    outputs = odeget(options,'OutputSel',1:neq);
    if isa(outputFcn,'function_handle')
        % With MATLAB 6 syntax pass additional input arguments to outputFcn.
        outputArgs = varargin;
    end
end
refine = max(1,odeget(options,'Refine',4));
if ntspan > 2
    outputAt = 1;          % output only at tspan points
elseif refine <= 1
    outputAt = 2;          % computed points, no refinement
else
    outputAt = 3;          % computed points, with refinement
    S = (1:refine-1) / refine;
end
printstats = strcmp(odeget(options,'Stats','off'),'on');

% Handle the event function
[haveEventFcn,eventFcn,eventArgs,valt,teout,yeout,ieout] = ...
    odeevents(odeIsFuncHandle,ode,t0,y0,options,varargin);

% Handle the mass matrix
[Mtype, M, Mfun] =  odemass(odeIsFuncHandle,ode,t0,y0,options,varargin);
if Mtype > 0  % non-trivial mass matrix
    Msingular = odeget(options,'MassSingular','no');
    if strcmp(Msingular,'maybe')
        warning(message('MATLAB:ode45:MassSingularAssumedNo'));
    elseif strcmp(Msingular,'yes')
        error(message('MATLAB:ode45:MassSingularYes'));
    end
    % Incorporate the mass matrix into ode and odeArgs.
    [ode,odeArgs] = odemassexplicit(odeIsFuncHandle,Mtype,ode,odeArgs,Mfun,M);
    f0 =ode(t0,y0,odeArgs{:});
    nfevals = nfevals + 1;
end

% Non-negative solution components
idxNonNegative = odeget(options,'NonNegative',[]);
nonNegative = ~isempty(idxNonNegative);
if nonNegative  % modify the derivative function
    [ode,thresholdNonNegative] = odenonnegative(ode,y0,threshold,idxNonNegative);
    f0 = ode(t0,y0,odeArgs{:});
    nfevals = nfevals + 1;
end

t = t0;
y = y0;

% Allocate memory if we're generating output.
prototype = zeros(dataType);
nout = 0;
tout = []; yout = [];
if nargout > 0
    if output_sol
        chunk = min(max(100,50*refine), refine+floor((2^11)/neq));
        tout = zeros(1, chunk, 'like', prototype);
        yout = zeros(neq, chunk, 'like', prototype);
        f3d  = zeros(neq, 7, chunk, 'like', prototype);
    else
        if ntspan > 2                         % output only at tspan points
            tout = zeros(1, ntspan, 'like', prototype);
            yout = zeros(neq, ntspan, 'like', prototype);
        else                                  % alloc in chunks
            chunk = min(max(100,50*refine), refine+floor((2^13)/neq));
            tout = zeros(1, chunk, 'like', prototype);
            yout = zeros(neq, chunk, 'like', prototype);
        end
    end
    nout = 1;
    tout(nout) = t;
    yout(:,nout) = y;
end

% Initialize method parameters.
pow = 1/5;
A = [1/5, 3/10, 4/5, 8/9, 1, 1]; % Still used by restarting criteria
% B = [
%     1/5         3/40    44/45   19372/6561      9017/3168       35/384
%     0           9/40    -56/15  -25360/2187     -355/33         0
%     0           0       32/9    64448/6561      46732/5247      500/1113
%     0           0       0       -212/729        49/176          125/192
%     0           0       0       0               -5103/18656     -2187/6784
%     0           0       0       0               0               11/84
%     0           0       0       0               0               0
%     ];
% E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40];

% Same values as above extracted as scalars (1 and 0 values ommitted)
a2=cast(1/5, 'like', prototype);
a3=cast(3/10, 'like', prototype);
a4=cast(4/5, 'like', prototype);
a5=cast(8/9, 'like', prototype);

b11=cast(1/5, 'like', prototype);
b21=cast(3/40, 'like', prototype);
b31=cast(44/45, 'like', prototype);
b41=cast(19372/6561, 'like', prototype);
b51=cast(9017/3168, 'like', prototype);
b61=cast(35/384, 'like', prototype);
b22=cast(9/40, 'like', prototype);
b32=cast(-56/15, 'like', prototype);
b42=cast(-25360/2187, 'like', prototype);
b52=cast(-355/33, 'like', prototype);
b33=cast(32/9, 'like', prototype);
b43=cast(64448/6561, 'like', prototype);
b53=cast(46732/5247, 'like', prototype);
b63=cast(500/1113, 'like', prototype);
b44=cast(-212/729, 'like', prototype);
b54=cast(49/176, 'like', prototype);
b64=cast(125/192, 'like', prototype);
b55=cast(-5103/18656, 'like', prototype);
b65=cast(-2187/6784, 'like', prototype);
b66=cast(11/84, 'like', prototype);

e1=cast(71/57600, 'like', prototype);
e3=cast(-71/16695, 'like', prototype);
e4=cast(71/1920, 'like', prototype);
e5=cast(-17253/339200, 'like', prototype);
e6=cast(22/525, 'like', prototype);
e7=cast(-1/40, 'like', prototype);

hmin = 16*eps(t);
if isempty(htry)
    % Compute an initial step size h using y'(t).
    absh = min(hmax, htspan);
    if normcontrol
        rh = (norm(f0) / max(normy,threshold)) / (0.8 * rtol^pow);
    else
        rh = norm(f0 ./ max(abs(y),threshold),inf) / (0.8 * rtol^pow);
    end
    if absh * rh > 1
        absh = 1 / rh;
    end
    absh = max(absh, hmin);
else
    absh = min(hmax, max(hmin, htry));
end
f1 = f0;

% Initialize the output function.
if haveOutputFcn
    feval(outputFcn,[t tfinal],y(outputs),'init',outputArgs{:});
end

if ~isempty(odeArgs)
    ode = @(t,y) ode(t,y,odeArgs{:});
end

% THE MAIN LOOP

done = false;
while ~done

    % By default, hmin is a small number such that t+hmin is only slightly
    % different than t.  It might be 0 if t is 0.
    hmin = 16*eps(t);
    absh = min(hmax, max(hmin, absh));    % couldn't limit absh until new hmin
    h = tdir * absh;

    % Stretch the step if within 10% of tfinal-t.
    if 1.1*absh >= abs(tfinal - t)
        h = tfinal - t;
        absh = abs(h);
        done = true;
    end

    % LOOP FOR ADVANCING ONE STEP.
    nofailed = true;                      % no failed attempts
    while true
        y2 = y + h .* (b11.*f1 );
        t2 = t + h .* a2;
        f2 = ode(t2, y2);

        y3 = y + h .* (b21.*f1 + b22.*f2 );
        t3 = t + h .* a3;
        f3 = ode(t3, y3);

        y4 = y + h .* (b31.*f1 + b32.*f2 + b33.*f3 );
        t4 = t + h .* a4;
        f4 = ode(t4, y4);

        y5 = y + h .* (b41.*f1 + b42.*f2 + b43.*f3 + b44.*f4 );
        t5 = t + h .* a5;
        f5 = ode(t5, y5);

        y6 = y + h .* (b51.*f1 + b52.*f2 + b53.*f3 + b54.*f4 + b55.*f5 );
        t6 = t + h;
        f6 = ode(t6, y6);

        tnew = t + h;
        if done
            tnew = tfinal;   % Hit end point exactly.
        end
        h = tnew - t;      % Purify h.

        ynew = y + h.* ( b61.*f1 + b63.*f3 + b64.*f4 + b65.*f5 + b66.*f6 );
        f7 = ode(tnew, ynew);

        nfevals = nfevals + 6;

        % Estimate the error.
        NNrejectStep = false;
        fE = f1*e1 + f3*e3 + f4*e4 + f5*e5 + f6*e6 + f7*e7;
        if normcontrol
            normynew = norm(ynew);
            errwt = max(max(normy,normynew),threshold);
            err = absh * (norm(fE) / errwt);
            if nonNegative && (err <= rtol) && any(ynew(idxNonNegative)<0)
                errNN = norm( max(0,-ynew(idxNonNegative)) ) / errwt ;
                if errNN > rtol
                    err = errNN;
                    NNrejectStep = true;
                end
            end
        else
            err = absh * norm((fE) ./ max(max(abs(y),abs(ynew)),threshold),inf);
            if nonNegative && (err <= rtol) && any(ynew(idxNonNegative)<0)
                errNN = norm( max(0,-ynew(idxNonNegative)) ./ thresholdNonNegative, inf);
                if errNN > rtol
                    err = errNN;
                    NNrejectStep = true;
                end
            end
        end

        % Accept the solution only if the weighted error is no more than the
        % tolerance rtol.  Estimate an h that will yield an error of rtol on
        % the next step or the next try at taking this step, as the case may be,
        % and use 0.8 of this value to avoid failures.
        if err > rtol                       % Failed step
            nfailed = nfailed + 1;
            if absh <= hmin
                warning(message('MATLAB:ode45:IntegrationTolNotMet', sprintf( '%e', t ), sprintf( '%e', hmin )));
                solver_output = odefinalize(solver_name, sol,...
                    outputFcn, outputArgs,...
                    printstats, [nsteps, nfailed, nfevals],...
                    nout, tout, yout,...
                    haveEventFcn, teout, yeout, ieout,...
                    {f3d,idxNonNegative}, ...
                    t);
                if nargout > 0
                    varargout = solver_output;
                end
                return;
            end

            if nofailed
                nofailed = false;
                if NNrejectStep
                    absh = max(hmin, 0.5*absh);
                else
                    absh = max(hmin, absh * max(0.1, 0.8*(rtol/err)^pow));
                end
            else
                absh = max(hmin, 0.5 * absh);
            end
            h = tdir * absh;
            done = false;

        else                                % Successful step

            NNreset_f7 = false;
            if nonNegative && any(ynew(idxNonNegative)<0)
                ynew(idxNonNegative) = max(ynew(idxNonNegative),0);
                if normcontrol
                    normynew = norm(ynew);
                end
                NNreset_f7 = true;
            end

            break;

        end
    end
    nsteps = nsteps + 1;

    if haveEventFcn
        f = [f1 f2 f3 f4 f5 f6 f7];
        [te,ye,ie,valt,stop] = ...
            odezero(@ntrp45,eventFcn,eventArgs,valt,t,y,tnew,ynew,t0,h,f,idxNonNegative);
        if ~isempty(te)
            if output_sol || (nargout > 2)
                teout = [teout, te]; %#ok<AGROW>
                yeout = [yeout, ye]; %#ok<AGROW>
                ieout = [ieout, ie]; %#ok<AGROW>
            end
            if stop                 % Stop on a terminal event.
                % Adjust the interpolation data to [t te(end)].

                % Update the derivatives using the interpolating polynomial.
                taux = t + (te(end) - t)*A;
                [~,f(:,2:7)] = ntrp45(taux,t,y,[],[],h,f,idxNonNegative);
                f2 = f(:,2); f3 = f(:,3); f4 = f(:,4); f5 = f(:,5); f6 = f(:,6); f7 = f(:,7);

                tnew = te(end);
                ynew = ye(:,end);
                h = tnew - t;
                done = true;
            end
        end
    end

    if output_sol
        nout = nout + 1;
        if nout > length(tout)
            tout = [tout, zeros(1, chunk, 'like', prototype)];  %#ok<AGROW> % requires chunk >= refine
            yout = [yout, zeros(neq, chunk, 'like', prototype)]; %#ok<AGROW>
            f3d  = cat(3,f3d,zeros(neq, 7, chunk, 'like', prototype));
        end
        tout(nout) = tnew; %#ok<AGROW>
        yout(:,nout) = ynew; %#ok<AGROW>
        f3d(:,:,nout) = [f1 f2 f3 f4 f5 f6 f7]; %#ok<AGROW>
    end

    if output_ty || haveOutputFcn
        switch outputAt
        case 2      % computed points, no refinement
            nout_new = 1;
            tout_new = tnew;
            yout_new = ynew;
        case 3      % computed points, with refinement
            tref = t + (tnew-t)*S;
            nout_new = refine;
            tout_new = [tref, tnew];
            yntrp45 = ntrp45split(tref,t,y,h,f1,f3,f4,f5,f6,f7,idxNonNegative);
            yout_new = [yntrp45, ynew];
        case 1      % output only at tspan points
            nout_new =  0;
            tout_new = [];
            yout_new = [];
            while next <= ntspan
                if tdir * (tnew - tspan(next)) < 0
                    if haveEventFcn && stop     % output tstop,ystop
                        nout_new = nout_new + 1;
                        tout_new = [tout_new, tnew]; %#ok<AGROW>
                        yout_new = [yout_new, ynew]; %#ok<AGROW>
                    end
                    break;
                end
                nout_new = nout_new + 1;
                tout_new = [tout_new, tspan(next)]; %#ok<AGROW>
                if tspan(next) == tnew
                    yout_new = [yout_new, ynew]; %#ok<AGROW>
                else
                    yntrp45 = ntrp45split(tspan(next),t,y,h,f1,f3,f4,f5,f6,f7,idxNonNegative);
                    yout_new = [yout_new, yntrp45]; %#ok<AGROW>
                end
                next = next + 1;
            end
        end

        if nout_new > 0
            if output_ty
                oldnout = nout;
                nout = nout + nout_new;
                if nout > length(tout)
                    tout = [tout, zeros(1, chunk, 'like', tout)]; %#ok<AGROW> requires chunk >= refine
                    yout = [yout, zeros(neq, chunk, 'like', yout)]; %#ok<AGROW>
                end
                idx = nout;
                if nout_new ~= 1
                    idx = oldnout+1:nout;
                end
                tout(idx) = tout_new; %#ok<AGROW>
                yout(:,idx) = yout_new; %#ok<AGROW>
            end
            if haveOutputFcn
                stop = feval(outputFcn,tout_new,yout_new(outputs,:),'',outputArgs{:});
                if stop
                    done = true;
                end
            end
        end
    end

    if done
        break
    end

    % If there were no failures compute a new h.
    if nofailed
        % Note that absh may shrink by 0.8, and that err may be 0.
        temp = 1.25*(err/rtol)^pow;
        if temp > 0.2
            absh = absh / temp;
        else
            absh = 5.0*absh;
        end
    end

    % Advance the integration one step.
    t = tnew;
    y = ynew;
    if normcontrol
        normy = normynew;
    end
    if NNreset_f7
        % Used f7 for unperturbed solution to interpolate.
        % Now reset f7 to move along constraint.
        f7 = ode(tnew, ynew);
        nfevals = nfevals + 1;
    end
    f1 = f7;  % Already have f(tnew,ynew)

end % end of while loop

interp_data = {f3d,idxNonNegative};
solver_output = odefinalize(solver_name, sol,...
    outputFcn, outputArgs,...
    printstats, [nsteps, nfailed, nfevals],...
    nout, tout, yout,...
    haveEventFcn, teout, yeout, ieout, interp_data);

if nargout > 0
    varargout = solver_output;
end

[Sierra Cagle]

Thank you for this help!