ODE OutputFcn returning fewer points than expected

Hi,

I'm using OutputFcn to return exta information from ode45 and whilst the returned signal looks as I expected it to it's only comprised of 93 points rather than 1001.

tvec is my tspan vector and contains 1001 points for a run of 0 to 5 seconds at 200 Hz.

The ODE setup 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, tvec, lr, ...
    krx, crx, adr, adr_dbz, Fxr_max));
% Run solver
[T,x] = ode45(odeFcn, tvec, x0, options);

The intention is that Fxr is calculated during integration and its value capped under certain conditions (I'm aware this is an abuse of the ODE solver). I want Fxr out so I repeat the same equations in the output function, which contains:

function status = myOutputFcn(t, x, flag, ix, tvec, lr, krx, crx, adr, adr_dbz, Fxr_max)

persistent Fxr_out

ix = interp1(tvec, ix, t); % interpolates ix at t within tvec

switch flag
    case 'init'        
        tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr)); % anti-dive tangent wrt body displacement

        xir = (x(1)+lr.*x(3)).*tanadr; % contact patch longitudinal displacement...
        xdotir = (x(4)+lr.*x(6)).*tanadr; % ..and velocity
        
        Fxr = krx.*(-xir-x(2))+crx.*(-xdotir-x(5)); % rear axle longitudinal force
        
        if ix ~= 0 % if longitudinal acceleration is negative...
            if Fxr < Fxr_max % if Fxr exceeds limit...
                Fxr = Fxr_max; % cap Fxr to Fxr_max
            end
        end

        Fxr_out = Fxr;
        
    case ''
        tanadr = tan(deg2rad(adr_dbz.*(x(1)+lr.*x(3))+adr)); % anti-dive tangent wrt body displacement

        xir = (x(1)+lr.*x(3)).*tanadr; % contact patch longitudinal displacement...
        xdotir = (x(4)+lr.*x(6)).*tanadr; % ..and velocity
        
        Fxr = krx.*(-xir-x(2))+crx.*(-xdotir-x(5)); % rear axle longitudinal force
        
        if ix ~= 0 % if longitudinal acceleration is negative...
            if Fxr < Fxr_max % if Fxr exceeds limit...
                Fxr = Fxr_max; % cap Fxr to Fxr_max
            end
        end

        Fxr_out = [Fxr_out, Fxr];
        
    case 'done' % when it's done
        assignin('base','Fxr',Fxr_out); % get the data to the workspace. 
        
end

status = 0;

As i say, the resulting data in Fxr is correct. However, it only contains a fraction of the expected number of points. Why is this?

Regards,

Simon.

Answer 1

[Jan]

The code looks okay. Check again if tspan has really 1001 time points.

Use the debugger to check, what's going on. Set a breakpoint in the output function and see, when it is triggered.

I'd stay away from using assignin to create a variable in the base workspace. What about:

function status = myOutputFcn(t, x, flag, ix, tvec, lr, krx, crx, adr, adr_dbz, Fxr_max)
persistent Fxr_out
...
status = 0;
switch flag
    case 'init'        
        ...
    case ''
        ...
    case 'done'
        % Nothing
    case 'flush'
        status  = Fxr_out;
        Fxr_out = [];
end
end

Then calling this function after the integration to flush the internal persistent buffer is save, because the data are actively requested. Poking it by assignin into the base workspace is susceptible as all global access to variables: If there is another forgotton assignin anywhere in your code, the confusion is perfect.

Note: tand() is nicer than tand(deg2rad()).

[Simon Aldworth]

Hi Jan,

Thanks for answering and for the tips.

I set a breakpoint in the output function and was very surprised at the results. I had expected that t would be a single value representing the most recent successful integration time step. However, what I see appears to be otherwise.

The first loop through gives the following variable values and statement returns:

tvec = [0,0.005,0.01,0.015 etc] (1x1001 double)

flag = 'init'

t = [0,5]

ix = [0,0]

Fxr = 0

if ix ~= 0 is FALSE

Fxr_out = 0

It appears to me that t is just showing the tspan range and so ix interpolates accordingly. Is this presumably just a response to the initialisation loop? I don't know how the if statement reacts to ix in vector format. As it happens, Fxr and Fxr_out contain the expected result. tvec is as expected.

Moving on to the second loop, the above list changes to:

tvec = [0,0.005,0.01,0.015 etc] (1x1001 double)

flag = ''

t = 0.005

ix = 0

Fxr = 0

if ix ~= 0 is FALSE

Fxr_out = [0,0]

This is entirely as expected. However, things become problematic on the third loop through:

tvec = [0,0.005,0.01,0.015 etc] (1x1001 double)

flag = ''

t = [0.0100,0.0150,0.0200,0.0250,0.0300]

ix = [0,0,0,0,0]

Fxr = 0

if ix ~= 0 is FALSE

Fxr_out = [0,0,0]

Now t is a vector of 5 values, rather than the expected single value (I assumed it would be just the first one, 0.01). I won't post any further loops, except to show that t then starts from the LAST element value from the previous loop:

t(1) = 0.0350

With t skipping forward like that it is little wonder I only get 93 points back. But why is this happening?

Many thanks,

Simon.

[Simon Aldworth]

Hi Jan,

From the help page for odeset:

The solver calls status = myOutputFcn(t,y,[]) after each integration step for which output is requested. t contains points where output was generated during the step, and y is the numerical solution at the points in t. If t is a vector, then the ith column of y corresponds to the ith element of t.

So, it appears I may get several results for t at any one time. I assume this means that I need to format the code in my output function to deal with multiple values in t?

Regards,

Simon.

[Simon Aldworth]

Hi Jan,

I have managed to reformat my code to deal with multiple values in t:

switch flag
    case 'init'        
        tanadr = tand(adr_dbz.*(x(1,:)+lr.*x(3,:))+adr); % anti-dive tangent wrt body displacement

        etc...

I've also changed the if statement to deal with multiple values in ix and Fxr:

idx = Fxr < Fxr_max & ix~=0;
Fxr(idx) = Fxr_max;

So, thanks for your help and guidance.

Just one last question: I can't find any help documentation for using 'flush' with flag. Where did that come from and is there a help page hidden somewhere?

Thanks and regards,

Simon.

[Simon Aldworth]

@Jan Hi, can you tell me about using 'flush' with flag - is there any documentation for that?

Thanks.