Using persistent variable to pass out extra parameters using ode45 from the MATLAB ODE suite?

0
votes

In How do I pass out extra parameters using ODE23 or ODE45 from the MATLAB ODE suite? the "MathWorks Support Team" suggests using a persistent variable to pass out extra parameters using ode45 from the MATLAB ODE suite.

However, what if an integration failed and solver will call your function for the previous time step? The ode45 solver uses an adaptive time step, and sometimes integration can fail, so the solver automatically decreases time step, and "has to go back."

Is "MathWorks Support Team" wrong in their suggestion of using a persistent variable? I see that the 'OutputFcn' method can be used to pass out extra parameters, but I do not know how to use it. Could you give me an outline or example code using the 'OutputFcn' method to detect failed steps / flag etc. to give the correct extra parameters? If I use the 'OutputFcn' method, is a global variable needed?

1
matlabodenumerical-integration
Your concerns are correct, and Output Functions are the right way to solve the problem. The documentation is in it.mathworks.com/help/matlab/ref/…. You will have to define a global variable in the Output Functions, and in the odefun. – Jommy
@Jommy Thanks. I have tried using outputFcn for many attempts but i do not have a general outline on how to use it and if how to define or which to define a global variable/s. Can you give me an outline or example code on how to achieve this? – Ka-Wa Yip

1 Answers

2
votes

Here is an example, slightly modified from the documentation.

Define the ODE function. param will be the variable to be exchanged.

function dy = rigid(t, y)

global param

dy = zeros(3,1);    % a column vector
dy(1) = param + y(2) * y(3);
dy(2) = -y(1) * y(3);
dy(3) = -0.51 * y(1) * y(2);

param = dy(1) * dy(3);

The output function must declare the variable to be exchanged as global. myout will be called after each successful step.

function status = myout(t, y, flag)

global param

% Don't ever halt the integration
status = 0;

if strcmp(flag, 'init')
    % Initialization (apparently useless).
    param = -2;

elseif isempty(flag)
     % Main code to update param.
     param = param * mean(t);

elseif strcmp(flag, 'done')
     % Nothing to do.

end

Call ode45 in this way.

% True initialization.
global param
param = -2;

odeOpts = odeset(  ...
    'RelTol', 1e-6,  ...
    'AbsTol', 1e-7,  ...
    'OutputFcn', @myout);

[time, y] = ode45(@rigid, [0 10], [0 1 1], odeOpts);

plot(  ...
    time, y(:,1), '-',   ...
    time, y(:,2), '-.',  ...
    time, y(:,3), '.');