Is there a way to output intermediate calculations of an ode solver in MATLAB?

Question

I have a system of stiff ODEs that I solve with ode15s. The time derivative is of the type:

yt = A*y - u.*y;

Where yt, y and u are column vectors and A is a sparse matrix. The computation of u has to be carried every evaluated time and it's computationally intensive (u is non-linear and it's the cause of the stiffness). I don't want to pre-compute u becase I would like to keep the freedom of the solver to choose any desired times later. What I'm interested in, however, is not on the solution y(t), but on u(t).*y(t), which is already calculated to compute yt. If I get the solution as a struct sol that then I can evaluate with deval, I can easily compute y(t), but u(t) would have to be recalculated. I was wondering if there is a way to somewhat store the result of u.*y while using the solver, so then I could call it as easily as I call deval(y,t). The closest alternative I have found, is to rely on deval also providing the time derivative, so then do:

[y,yt] = deval(sol,t); %t is a row vector
uy     = A*y - yt;     %uy = u.*y

Which is not terribly inefficient, but wastes some computation in operations that were already computed.

Any idea on how to do this?

Answer 1

As far as I have researched this exactly same problem I had a few months ago, this is not possible as per built-in functions. There are some workarounds to save the intermediate calculations (like using persistent or global variables), but they create more problems.

There is a lot of discussion on MATLAB's forum. This answer may be particularly useful. Other comments may be found here or here.

The problem to save intermediate calculations with MATLAB built-in ode functions is that they are adaptive step size and the function to be integrated is called several times per step. If you pass t as a vector with more than 2 elements, MATLAB will return the solution at the time steps provided. That does not mean, however, it calculates the solution on those time steps - it just interpolates from their built-in solution.

The best option is what you mentioned, to calculate the result you want after solving the ode, which means recalculating u. If you properly vectorize your functions, this step should not be time problematic.

Answer 2

You could "hack" the ODE solvers and include the control variables u in the state vector. Then your ODE system becomes a (trivial, explicit) DAE system

y' = f(t,y,u)
0  = u - g(t,x)

The only "trick" now is that you also need to pass a mass matrix in the options that is a diagonal matrix with entries 1 for the differential equations and 0 for the algebraic equations. Check the documentation if the solver you use supports DAE, all implicit methods should provide this.

There may be some efficiency loss as the new equations will also influence the step size control, but I would not expect a large difference.

In the result, the values of the control will be stored in the last components of the state vector, so extracting them is a simple slicing of that vector.