matlab ode45 retrieving parameters

Question

I'm experimenting with ode45 in Matlab. I've learned how to pass parameters to the ode function but I still have a question. Let's suppose that I want to compute the trajectory (speed profile) of a Car and I have a function, e.g. `getAcceleration`, that gives me the acceleration of the car but also the right gear: `[acceleration, gear] = getAcceleration(speed,modelStructure)` where `modelStructure` represents the model of the car. The ode function would be: ``` function [dy] = car(t,y,modelStructure) dy = zeros(2,1); dy(1) = y(2); [dy(2),gear] = getAcceleration(y(1),modelStructure); ``` Then I call the Ode45 integrator in this way: ``` tInit = 0; tEnd = 5, [t,y] = ode45(@car,[tInit tEnd], [speedInitial,accelerationInitial],options,modelStructure); ``` The problem is: how do I get the vector storing gears? Should I have something like `[t,y,gear]=ode45(....)` or should `gear` be within the `y` vector? I've been working on my code and using the events function I'm now able to get the ca…

Accepted Answer

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.

matlab ode45 retrieving parameters

Solution

Problem