How to save a variable computed inside ODE?

Inside an ode function I perform some calculations to compute a variable. I would like to save the values of that variable in a matrix so after the ode is done I can plot it vs time. Because it would not be useful to save the values during the integration, since the function is evaluated several times for each step and a step could be rejected, I followed michio's suggestion to use an Output function.

His solution is also described here: Get variable out of ode 45

function status = myOutputFcn(time, x, flag)
% q: value of variable at each succesfull time step
% q_save: matrix of variable values I want to save for plotting after ode is done

persistent q_save
switch flag
    case 'init'
        q_save = q;
    case ''
        q_save = [q_save; q];
    case 'done' % when it's done
        assignin('base','q_save',q_save); % get the data to the workspace.
end
status = 0;
end

The problem I have with this approach is that after the ode finishes, the size of the time vector is different from the size of the variable 'q_save'. The Output function is called less times than the succesfull steps. For example, time=800x1 and q_save=150x1.

I have also tried declaring tspan in two ways

tspan = linspace(t0,tf,100)  
tspan = [t0 tf]

So, how can I save a variable computed inside ODE for every succesfull ode step?

Meaning that I want a value of q for every (time,x) pair that the ode gives me as a solution.

Answer 1

[Jan]

assignin is a bad programming technique. Producing variables remotely in another workspace is not clean. A better method is to let the integration run regularily and request the wanted values afterwards:

[t, y] = ode45(@fcn, tSpan, y0);

function dy = fcn(t, y)
p  = sin(y(2));   % <- The wanted parameter
dy = [2 * t; p / y(1)];
end

If this is the function to be integrated, it must be modified such, that it accepts vectors as input and replies p as output:

function [dy, p] = fcn(t, y)
p  = sin(y(2, :));   % <- The wanted parameter  [EDITED] y(:,2) --> y(2,:)
dy = [2 .* t; p ./ y(1, :)];                 %  [EDITED] y(:,1) --> y(1,:)
end

Then:

[t, y] = ode45(@fcn, tSpan, y0);
[~, p] = fcn(t.', y.');

This is much easier than letting the output function create variables remotely.

For some cases the output function is the best choice. Then provide the value as output instead of assignin

function status = myOutputFcn(time, x, flag)
persistent q_save
status = 0;
switch flag
    case 'init'
        q_save = q;
    case ''
        q_save = [q_save; q];
    case 'done' % when it's done
        status = q_save;
        q_save = [];
end
end

But now to the actual problem: Unfortunately the ODE integrators reply different time steps, if the output is obtained as a struct. You did not show, how you call the integrator, so the most important detail is missing. (In opposite to this, the code or the output function is not a part of the problem.)

If you want to call the integrator in the form sol = ode45(___), and use a tspan defined as vector also, you have to modify ode45 such that the Output function is called regularily as for the [t,y,te,ye,ie] = ode45(odefun,tspan,y0,options) form. I consider this a bug and sent a bug report to MathWorks some years ago.

You wrote:

"I have also tried declaring tspan in two ways"

  tspan = linspace(t0,tf,100)  
  tspan = [t0 tf]

And what did you observe as effect? How did you call which integrator?

[elena]

Thank you so much for your response.

• "And what did you observe as effect?"

The size of the steps is different (size of time) but in both cases the size of time and q_save does not mach.

• "assignin is a bad programming technique"

Actually I did not use assignin (it just was in the solution I copied) but thanks for letting me know.

• The way I call my function is:

options = odeset('Event',@Event_Function,'OutputFcn',@(t,xp,flag) myOutputFcn(t,xp,flag));
[t,xp,te,ye,ie] = ode45(@(t,xp) my_ode(t,xp),tspan,xp0,options);

Does this mean that there is a problem with that?

• I tried to put 2 outputs in the ode function like you suggest and call it like that:

[t, y] = ode45(@fcn, tSpan, y0); 
[~, p] = fcn(t(end), y(end,:)); % example for one value of t to show size. 
                                % I will put it in a loop  

function [dy, p] = fcn(t, y)
p  = sin(y(2));   % <- The wanted parameter that has size 4x1
dy = ..
end

(I did not use vectorization since I have a lot of computations where I already use vectorization and it would be difficult for me to do it like that. Also, at this step, I don't mind about performance so I can use a for loop after the ode is done)

The problem is that even though I compute a 4x1 variable in the ode, when I call it afterwards, I get a 4x4 p variable.

Why is that?

[Jan]

You have to provide y as column vector:

p = zeros(size(t));
for ip = 1:numel(p)
  [~, p(ip)] = fcn(t(ip), y(ip, :).');
  %                               ^^
end

"The size of the steps is different (size of time) but in both cases the size of time and q_save does not mach." - The actual values might contain useful information, so please post the details. Actually ode45 replies as many steps as you define in tspan, if it is a vector. This is a bad design, because it needs a special behaviour, if tspan is meant as a [1 x 2] vector, but than ode45 treats it as time span. Then there is the known bug, that ode45 ignores a tspan defined as vector, if the output is caught as a struct. But this is not te case for you.

This means, that the cause of your problem is still unknown. Maybe there is another mistake anywhere else. If you post the sizes of tspan and the replied t, this might shed light on the problem.

[elena]

Oh ok.

Thank you very much for your help!

[elena]

I post my results in case they are of any help.

If I use

tspan = [0 10];

I get time=317x1 and q_save=80x1

If I use

tspan = linspace(0,10,5000);

I get time=5000x1 and q_save=76x1

If I use

tspan = linspace(0,10,2000);

I get time=2000x1 and q_save=76x1

If I use

tspan = linspace(0,10,1000);

I get time=1000x1 and q_save=75x1

If I use

tspan = linspace(0,10,100);

I get time=100x1 and q_save=71x1

[Akash Singh]

Hey @Jan

Shouldn't it be ?

p  = sin(y(2,:));   % <- (Instead of y(:,2))
dy = [2 .* t; p ./ y(1,:)];

[Mohammad]

I have the following code and I want to get F, G , U , xdot , x1 , x2 and any other variable in the code as outputs but I dont know how. I need them to make other functions and I need their exact dimensions because they are vectors and in some cases matrices . I appreciate it if you could help me out.

function [] = osc()

tspan = 0 : .1 : 20;

x0 = [3;-2];

[t,x] = ode45(@osc2d_2,tspan,x0);

x

plot(t,x)

function xdot = osc2d_2(t,x)

xdot = f(t,x) + g(t,x) * u(t,x)

end

function F = f(t,x)

F = [x(2) ; -x(1) * (pi/2 + atan(5*x(1))) - (5*x(1).^2/(2*(1+25*x(1).^2))) + 4*x(2)]

end

function G = g(t,x)

G = [0;3]

end

function U = u(t,x)

U = -3*x(2)

end

end