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

26 views (last 30 days)
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:
Any syggestions?
  1 Comment
Jan
Jan on 20 Nov 2019
Edited: Jan on 20 Nov 2019
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

Sign in to comment.

Accepted Answer

Jan
Jan on 20 Nov 2019
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
  1 Comment
oski89
oski89 on 20 Nov 2019
Thank you Jan for the insanely quick reply!
I will definitely look closer into this tomorrow.

Sign in to comment.

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!