Store value of variable computed in ODE Solver at the specified time steps.

I am simulating a system using the ode15s(..) solver. Code is like this:

tspan = 0:0.01:10;
options = odeset('RelTol',1e-7,'AbsTol',1e-10,'Refine',4);
[t,theta] = ode15s(@ode_simulation, tspan, theta_init, options);

function d_theta = ode_simulation(t,theta)

    u = k*(p-d)/e;
    d_theta = [theta(2) ; -m*c*[theta(2) theta(4)]'-m*g+m*u ; theta(4) ; -m*c*[theta(2) theta(4)]'-m*g+m*u];

end

The solver returns the results in t, theta variable at the time steps specified by tspan. I want to find a way to obtain the values of variable u that is being calculated inside the ode_simulation(..) function. I don't want to define a global array and keep inserting values cause it is computationally-wise and memory-wise terrible. So, after the simulation is exefuted I want to have the t,theta variables and in addition a vector u contaning the computed u values for the same time steps. How could I do that ?

Answer 1

[Star Strider]

There are not enough variables provided to run your code.

This is how I would do it:

tspan = 0:0.01:10;
options = odeset('RelTol',1e-7,'AbsTol',1e-10,'Refine',4);
[t,theta] = ode15s(@ode_simulation, tspan, theta_init,options);

function [d_theta,u] = ode_simulation(t,theta)

    u = k*(p-d)/e;
    d_theta = [theta(2) ; -m*c*[theta(2) theta(4)]'-m*g+m*u ; theta(4) ; -m*c*[theta(2) theta(4)]'-m*g+m*u];

end

for k = 1:numel(t)
    [~,u(k)] = ode_simulation(t(k),theta(k,:));
end

figure
plot(t, theta,    t, u)
grid
legend('\theta(t)', 'u(t)')

.

[Teo Protoulis]

I have edited and linked pastebins contaning the whole code. Coded your initial suggestion and got the 'Invalid expression. Check for missing or extra characters.' error.

[Star Strider]

I do not understand where the error came from, since my code worked when I simulated it (with a different differential equation). It would seem to be a typographical error that would be straightforward to find and correct.

One problem in your code is the extensive use of global variables. Instead, it would be best to assign them and then pass them as extra parameters in your functions. See Passing Extra Parameters for a full description. (I have not run it, since I do not understand it or how it all fits together.)

[Teo Protoulis]

One problem is also the fact that while the simulation runs the theta variables get quantized before the signal u is computed. If I compute the signal u after the simulation is executed by using the solutions theta of the simulation, then these solutions will be again quantized and the resulting signal u will not be the same as the one which was calculated during the initial simulation.

[Star Strider]

That should not make a difference using the loop. It uses the ‘t’ and ‘theta’ values already calculated. The only difference is that it produces the ‘u’ value at each step.

I have done this previously with other differential equation functions, and it has worked without problems.

One potentially serious problem in the original code is the extensive use of global variables. It is best to avoid them at all costs. If you are getting different results in the loop than from the original result, that may be the reason. I linked to a documentation section in my previous Comment that tells how to pass extra parameters to functions.

[Teo Protoulis]

Yes, I am already working towards getting rid of the global variables entirely. Gonna test it and comment regarding the results.

[Teo Protoulis]

The code executed perfectly as you suggested. Only thing I am not sure about is the issue I wrote previously regarding the quantization of the values. Even though it produces the 'u' values I am not sure if they are the identical ones with the first execution of the simulation.

[Star Strider]

Thank you!

The ‘u’ values should be as originally calculated, since all the calculations are the same. The ‘d_theta’ outputs will of course be the derivatives, so they will not match the ‘theta’ values from the integration, however when plotted against the ‘theta’ values, it should be obvious that they are the derivatives.

You can also check the derivatives of the ‘theta’ values by calculating them as:

d_theta_check = gradient(theta) ./ gradient(t);

They should be approximately equal to the ‘d_theta’ values when plotted or otherwise compared. That will also serve as an indirect check on the accuracy of the ‘u’ vector.

[Teo Protoulis]

And yes I tested it. Results are different. Obtaining the values of 'u' after the simulation has ended produced complex numbers for the isgnal. I coded an `if` statement inside the 'ode_simulation' function forcing it to pause if it detects complex numbers for u and it didn't pause at all. So, I think I have to find another way of taking the values of 'u'. If you have any ideas, I would really appreciate it if you shared them.

[Stephen23]

@Teo Protoulis : the solution provided by Star Strider is the recommended approach. Trying to store values while the ODE solver is running is not recommended because the ODE solver can use values not on the requested steps, can use values that are very close, and can backtrack and try new values. It is not trivial to keep track of this and return only the relevant values. It is much simpler to use the ODE solver's output values, just as Star Strider showed.

If you are getting unexpected complex values then you could start by investigating where that complexity first arises.

[Teo Protoulis]

@Stephen Cobeldick : I am aware of everything you wrote but still when the simulation runs I don't get any complex numbers at all. Not any component of the solutions includes any complex number. And when I am running again the function to store the u values I get the complex numbers. This is kind of weird considering also that I have ckecked the derivatives and they are very much alike. Anyway, thank you both for the time, learned a bunch of new things. Hope, I will discover where this complexity arises from.

[Star Strider]

Stephen Cobledick — Thank you!

Teo Protoulis — I have no idea what the problelm is or the reason you are not getting the desired result. I would check the magnitudes of the imaginary components. If they are vanishingly small (on the order of 1E-10 or less), they could simply be the result of floating-point calculation approximation errors, and not significant. In that event, you can safely ignore them, using the real() values of the vectors.

Also, note that one complex value calculated for a vector results in the entire vector being complex, even if the imaginary values for the other elements are 0.

[Teo Protoulis]

I have checked these conditions. There are more than one complex values and their magnitude can't be ignored. Anyway, I will keep checking. Thank you very much for your time and the useful information !