Problem using OutputFcn to get parameters out of ODE45

Hi,

I'm using ode45 to solve a three degree of freedom mass, spring, damper problem. It's been running well until I've modified the code to extract parameters calculated during the solving process using odeset with OutputFcn. I've used the following previous questions for help:

https://uk.mathworks.com/matlabcentral/answers/305698-how-do-i-extract-extra-parameters-from-the-ode45-ordinary-differential-equation-solver-function

https://uk.mathworks.com/matlabcentral/answers/312549-get-variable-out-of-ode-45

I am now receiving an error saying that an index exceeds the number of array elements (0) and the line in question is in myOutputFcn. It is presumably x that is returning the error:

tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr));

Where my code may differ from the examples I've used is that the calculations in myOutputFcn use the current values of the displacement and velocity variables that are being solved for in the ode. x(1) is the body vertical displacement and x(3) is the body pitch displacement.

The solver setup and calling code is:

% Set up ODE function to pass to solver
odeFcn = @(t, x) Three_DOF_Half_Car_ODE(t, x, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
    kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max);
% Set up options for output function
options = odeset('OutputFcn',@(t, x, flag) myOutputFcn(t, x, flag, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
    kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max));
% Run solver
[T,x] = ode45(odeFcn, tvec, x0, options);

The ode contains, for example:

function dx = Three_DOF_Half_Car_ODE(t, x, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
    kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max)

ix = interp1(tvec, ix, t); % this is interpolating ix at t within tvec
iz = interp1(tvec, iz, t);
ip = interp1(tvec, ip, t);
Fxr_max = interp1(tvec, Fxr_max, t);

% Calculation of anti-dive tangent wrt body displacement
tanadf = tan(deg2rad(adf_dbz.*(x(1)-lf.*x(3))+adf));
tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr));

% Contact patch longitudinal displacement and velocity equations
xif = -(x(1)-lf.*x(3)).*tanadf;
xir = (x(1)+lr.*x(3)).*tanadr;
xdotif = -(x(4)-lf.*x(6)).*tanadf;
xdotir = (x(4)+lr.*x(6)).*tanadr;

% Calculate rear axle longitudinal force and limit it according to decel
% Only valid during braking, when fxr is negative
Fxr = krx.*(-xir-x(2))+crx.*(-xdotir-x(5));
          
if Fxr < Fxr_max
    Fxr = Fxr_max;
end

dx(1) = x(4);
dx(2) = x(5);
dx(3) = x(6);
dx(4) = 1./mb.*(kfz.*(-x(1)+lf.*x(3))+cfz.*(-x(4)+lf.*x(6))+krz.*(-x(1)-lr.*x(3))+crz.*(-x(4)-lr.*x(6))-...
    (kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))).*tanadf+Fxr.*tanadr+mb.*iz);
dx(5) = 1./mb.*(kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))+Fxr+mb.*ix);
dx(6) = 1./Iby.*(-kfz.*lf.*(-x(1)+lf.*x(3))-cfz.*lf.*(-x(4)+lf.*x(6))+krz.*lr.*(-x(1)-lr.*x(3))+crz.*lr.*(-x(4)-lr.*x(6))-...
    kfx.*h.*(-xif-x(2))-cfx.*h.*(-xdotif-x(5))-Fxr.*h+...
    lf.*(kfx.*(-xif-x(2))+cfx.*(-xdotif-x(5))).*tanadf-lr.*Fxr.*tanadr+Iby.*ip);

It is the variable Fxr that I am trying to extract, so the following is in myOutputFcn (including the line beginning tanadr at which the code fails):

function status = myOutputFcn(t, x, flag, ix, iz, ip, tvec, mb, Iby, lf, lr, h, ...
    kfz, cfz, kfx, cfx, adf, adf_dbz, krz, crz, krx, crx, adr, adr_dbz, Fxr_max)

persistent Fxr

Fxr_max = interp1(tvec, Fxr_max, t);

tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr));

xir = (x(1)+lr.*x(3)).*tanadr;
xdotir = (x(4)+lr.*x(6)).*tanadr;

Fxr_tmp = krx.*(-xir-x(2))+crx.*(-xdotir-x(5));
          
if Fxr_tmp < Fxr_max
    Fxr_tmp = Fxr_max;
end

switch flag
    case 'init'

        Fxr = Fxr_tmp;

    case ''
        
        Fxr = [Fxr, Fxr_tmp];

    case 'done'
        
        assignin('base','Fxr',Fxr);
        
end

status = 0;

I'd be grateful for any guidance.

Regards, Simon.

[Jan]

By the way:

if Fxr < Fxr_max
    Fxr = Fxr_max;
end

This looks like your function to be integrated is not smooth. Then the computed trajectory can be dominated by rounding errors or the integration can stop with an error. Using ODE45 for non-smooth functions replies a trajectory, but this is not a reliable result from a scientifiv point of view.

Use an event function to stop and restart the integration instead.

The same problem concerns the linear interpolation by interp1.

It is still not clear to me, what you want to achieve. You calculate the trajectory already. As far as I can see, you try export one component a second time, but why?

[Simon Aldworth]

Hi Jan,

Yes, I understand that what I'm asking of ODE45 isn't strictly correct, but I've had (sufficient) success in the past so am sticking with it for now. I'm afraid I don't know how to stop and restart the integration so that's something else to learn.

To be clear: it is the value of Fxr as calculated in the ODE at each time step (as per tvec) that I want to keep.

I don't follow what you mean by "you calculate the trajectory already". What trajectory? What calculation are you referring to?

What do you mean by "a second time"? What "one component" am I exporting a first and second time?

Thanks.