First off: why are the initial values to the differential equation the initial speed (speedInitial
) and the initial acceleration (accelerationInitial
)? That means that the differential equation car
will be computing the acceleration (y(1)
) and the jerk (y(2)
), the time-derivative of the acceleration, at each time t
. That seems incorrect...I would say the initial values should be the initial position (positionInitial
) and the initial speed (speedInitial
). But, I don't know your model, I could be wrong.
Now, getting the gear
in the solution directlty: you can't, not without hacking ode45
. This is also logical; you also cannot get dy
at all times directly, can you? That's just not how ode45
is set up.
There's two ways out I see here:
Global variable
DISCLAIMER: don't use this method. It's only here to show what most people would do as a first attempt.
You can store gear
in a global variable. It's probably the least amount of coding, but also the least convenient outcome:
global ts gear ii
ii = 1;
tInit = 0;
tEnd = 5,
[t,y] = ode45(...
@(t,y) car(t,y,modelStructure), ...
[tInit tEnd], ...
[speedInitial, accelerationInitial], options);
...
function [dy] = car(t,y,modelStructure)
global ts gear ii
dy = zeros(2,1);
dy(1) = y(2);
[dy(2),gear(ii)] = getAcceleration(y(1),modelStructure);
ts(ii) = t;
ii = ii + 1;
But, due to the nature of ode45
, this will get you an array of times ts
and associated gear
which contains intermediate points and/or points that got rejected by ode45
. So, you'll have to filter for those afterwards:
ts( ~ismember(ts, t) ) = [];
I'll say it again: this is NOT the method I'd recommend. Only use global variables when testing or doing some quick-n-dirty stuff, but always very quickly shift towards other solutions. Also, the global variables grow at each (sub-)iteration of ode45
, which is an unacceptable performance penalty.
It's better to use the next method:
Post-solve call
This is also not too hard for your case, and the way I'd recommend you to go. First, modify the differential equation as below, and solve as normal:
tInit = 0;
tEnd = 5,
[t,y] = ode45(...
@(t,y) car(t,y,modelStructure), ...
[tInit tEnd], ...
[speedInitial, accelerationInitial], options);
...
function [dy, gear] = car(t,y,modelStructure)
dy = [0;0];
dy(1) = y(2);
[dy(2),gear] = getAcceleration(y(1),modelStructure);
and then after ode45
completes, do this:
gear = zeros(size(t));
for ii = 1:numel(t)
[~, gear(ii)] = car(t(ii), y(ii,:).', modelStructure);
end
That will get you all the gears the car would have at times t
.
The only drawback that I can see here is that you'll have many more function evaluations of car
than ode45
would use by itself. But this is only a real problem if each evaluation of car
takes in the order of seconds or longer, which I suspect is not the case in your setup.