ode45 - Plotting temporarily calculated quantity

I am using ode45 to integrate a second order differential equation describing state vector Z:

[t Z]= ode45(@diffEq,tSpan,Z0);

My function diffEq takes input Z and outputs its time derivative Zdot at time t. To do this, it needs to calculate a quantity called T:

Zdot= Z\T, for instance.

So after running ode45 I am able to plot the integration results (namely Z with respect to time). However, I am also interested in plotting T vs time but have been having some trouble.

At first I stored a time history of T using a global variable, but soon realized that these values do not correspond to those calculated at times specified in tSpan. Instead ode45 chooses its own time steps. I really need a time history of T that, like Z, matches with the tSpan I specified.

I have been trying to incorporate T into my state vector Z, but I haven't been able to make this work and am not sure how to do it. I don't think this will work because ode45 just keeps on integrating T.

I may try taking the global variable time history and trimming out the entries that don't correspond to tSpan. However, this seems like a very inelegant solution, and I'd like some help trying to come up with a better fix.

Thanks in advance.

Answer 1

[Teja Muppirala]

If you are creating T inside diffEq somewhere, then you could just as easily create T again after the solver is all done. In other works, after you run this:

[t Z]= ode45(@diffEq,tSpan,Z0);

You have t, and you have Z, so just just say T = whatever(t,Z) where "whatever" means the code you used inside the diffEq function.

For example if diffEq was like:

function Zdot = diffEq(t,Z)

T = [cos(t) sin(t); cos(t) -sin(t)]*sum(Z.^2)
Zdot =  Z\T;

Then instead of trying to get T from out of there, wait until the ODE solver is done and then use

T = [cos(t) sin(t); cos(t) -sin(t)*sum(Z.^2)

again on your output t and Z

[Kenny]

I see what you mean. That certainly works and is better than the workaround I've been doing. Thanks so much for your thoughts.

I was hoping there would be an easy way to extract T out of diffEq, but I think this is acceptable for now. Most of the "work" is already done by ode45 anyway, so this method doesn't make my program that much longer to run. However, if anyone finds a way to do this without re-running diffEq, please let me know.

Thanks.

[Richard Crozier]

One of the major annoyances of the solvers, in my opinion, is that you cannot extract internally calculated values. It could possibly be done, however, using global variables and the 'OutputFcn' parameter. This would be achieved by making T a global variable. In the OutputFcn, you could in principle discover the last calculated value of T when the solver was called at a given time step by looking at the contents of T (you may need to be more careful than this in reality) and storing this value, in yet another global array for access later or saving to disk.

[Teja Muppirala]

Well, it's not quite that simple. The time-steps that get output by ODE45 are not necessarily time steps where it integrated at (i.e., where it calculates the function dydx). It interpolates the output based on a number of integration steps. The OutputFcn, however, only gets called after an actual integration step, and so in general, the time steps where OutputFcn is called, and the time values that are output to the user will not be the same.

Even if you specify the integration time as a vector (0:0.1:10), it will give you outputs at those times, it does not necessarily evaluate dydx at those times. Because again, it may be interpolating those values. You can actually view the source for ODE45 and verify all of this.