How do I return an extra parameter using Matlab's ODE solvers?

I would like to return a parameter, r, from the ode solution, but this parameter is not to be integrated. However, it is important to the result and for later post-processing. I want r to be a vector matching the length of t so I can do

plot(t, r)

The tricky part for me that seems to stand out from the examples I've looked at here at the forum is that the time dependent variables and their derivatives are dependent on r which in turn r is dependent on a boolean state variable. This means it is rather difficult to calculate r based on the output results from the ode solution variables.

My code is

clear; close all; clf; clc

t_span  = [0 1];
y_0     = [0, 2];
env.g   = 9.81;

figure(101); hold on

options = odeset('AbsTol',1e-6, 'RelTol',1e-6, 'Stats','on');
[t, y] = ode45(@my_ode_fcn, t_span, y_0, options, env);

figure(1)
plot(t, y(:, 1))
xlabel('time, t [s]')
ylabel('position, x [m]')

figure(2)
plot(t, y(:, 2))
xlabel('time, t [s]')
ylabel('velocity, u [m/s]')

function dydt = my_ode_fcn(t, y, env)

    persistent state_flag
    if isempty(state_flag); state_flag=0; end

    x = y(1);
    u = y(2);
    g = env.g;

    if ~state_flag
        r = 20 + 50*t;
    else
        r = 20 + 10*t;
    end

    if (r > 30) && ~state_flag
        state_flag = 1;
        sprintf('State flag = %d', state_flag)
    end

    if (r > 27)
        xdot = 2*u;
        ydot = -g*x;
    else
        xdot = 3*u;
        ydot = -2*g*x;
    end

    figure(101)
    plot(t, r, '*')

    dydt    = zeros(length(y),1);
    dydt(1) = xdot;
    dydt(2) = ydot;

end

I've tried all the suggestions in this thread, but nothing has worked so far:

https://se.mathworks.com/matlabcentral/answers/92701-how-do-i-pass-out-extra-parameters-using-ode23-or-ode45-from-the-matlab-ode-suite

Any syggestions?

[Jan]

You are using ODE45 to integrate a non-mooth function. Matlab's ODE integrators are not specified to integrate over discontinuities. In consequence the step-size controller drives crazy. Unfortunately you might get a final value, but this is not a scientific computation any more and the accumulated rounding error of the trajectory can dominate the trajectory. ODE45 discontinuities

Use event function to handle discontinuities. You have to restart the integration after each change of the function.

Setting the stateflag inside the function to be integrated and using a persistent variable to store it is a mistake also: Remember, that the ODE integrator can reject steps, if the local rounding error exceed the tolerances. This means, that the state migt be switched in a rejected step and is invalid in the following shorter time step.

Use anonymous functions to provide parameters: Answers: Anonymous for params

Answer 1

[Jan]

If the problems mentioned in my comment are fixed, the solution is easy: Run the integration at first. Then provide the calcualted time and trajectory as input for my_ode_fcn and export the wanted value as 2nd output. Therefore the function must be vectorized - or called in a loop for each time point.

See

• https://www.mathworks.com/matlabcentral/answers/92701-how-do-i-pass-out-extra-parameters-using-ode23-or-ode45-from-the-matlab-ode-suite
• https://www.mathworks.com/matlabcentral/answers/396633-how-to-store-a-intermediate-variable-from-ode45-function-to-the-workspace-as-a-vector

[oski89]

Thank you Jan for the insanely quick reply!

I will definitely look closer into this tomorrow.