How to return an intermediate variable from ode45 in Matlab

What I want

Hi there.

I want to store or return an intermediate calculations in global variable at spacific time instants. For example.

%%%%%  Caller Code

x0 = [1 ; 2];
t = 0:0.1:10
[t,x] = ode45(@(t,x) odefun(t,x), t, x0);


%%%%%  ODE Function
function dx = odefun(t,x)

    if t<3
    e = x(1)-x(2) + 2;   % variable of my intrest, to be stored or returned each 0.1 sec 
    
    elseif t<6
    e = x(1)-x(2) - sin(t);
    
    elseif t<10
    e = x(1)-x(2)/x(1);
    end
    
    dx(1) = -x(2)
    dx(2) = sin(2*pi*t)
end

My requirement is that e should be calculated at exact time instants at which x is calculated. this way both would have the same length.

My Approch so far ,

%%%%%  Caller Code

global e indx;
e=0;
indx=1;

x0 = [1 ; 2];
t = 0:0.1:10
[t,x] = ode45(@(t,x) odefun(t,x), t, x0);


%%%%%  ODE Function
function dx = odefun(t,x)

    global e indx;
    
    if t<3
    e(indx) = x(1)-x(2) + 2;   % variable of my intrest, to be stored or returned each 0.1 sec exactly 
    
    elseif t<6
    e(indx) = x(1)-x(2) - sin(t);
    
    elseif t<10
    e(indx) = x(1)-x(2)/x(1);
    end
    
    dx(1) = x(2)
    dx(2) = sin(2*pi*t)
    
    indx = indx+1;
end

But using this method the length of e is extremly large, for example, x length is 200, then i have length of e more than 10000. My origional code is way more complex and simulated for longer time, which makes the process exreamply diffficult and inaccurate.

Answer 1

[Walter Roberson]

My requirement is that e should be calculated at exact time instants at which x is calculated. this way both would have the same length.

That request is based upon a misunderstanding of how the ode*() functions work.

The ode*() functions are all variable step-sized functions. At any one time, they are working with a current set of boundary conditions. They evaluate the provided function at a series of carefully chosen locations "near" the current set, evaluating with different combinations of changes to the time variable and changes to the boundary conditions. They use the calculated value to make a prediction about the value of the function at another nearby location, and then they ask for the calculated value at that location. If the prediction and the calculated value are close enough then the ode* routines "accept" a calculated location -- and increase the step size. If the prediction and the calculate value are not close enough, then the ode* routines reject the current step and start again from the same location but with a smaller step size; in areas that are rapidly changing, they might reject several proposed step sizes in a row.

The output values you are getting for e are the values at all of the locations, including the ones that were never expected to be correct but are being evaluated to figure out what the slope is doing locally -- and also including the locations that were used for the cross-check of the prediction.

In the case when you have specified a vector of more than two time points, the ode*() routines seldom evaluate the function at the given times. Instead, when they detect that a step has successfully crossed between requested times, they interpolate the output for the given time. So even if you dug through the entire list looking for the e values appropriate for those times, you might not find it, since the function was only evaluated "near" those times.

So... unless you want to do some kind of analysis of the surface at all the locations that were probed, then do not bother to record intermediate variables as you are going along. Instead, after you have completed the ode*() call and gotten the times and boundary condition outputs, run code that calculates the intermediate variables at those locations. One way to do that is to make the variable of interest the second output of your ode function; the ode*() routines will ignore the extra output, but you can save the output when you call your ode function "yourself" with the combinations of times and boundary conditions.

[Muhammad Qaisar Ali]

Hi Walter! Thank you very much for such a detailed answer. You are very much right. My method is very inaccurate and also computationaly complex. In one simulation when I stored a figure on which I had plotted e, that figure was alsmot 2 Giga Bytes.

Your answer solved all my misconsecptions, but I am a bit confused about the solution proposed by you at the last,

One way to do that is to make the variable of interest the second output of your ode function; the ode*() routines will ignore the extra output, but you can save the output when you call your ode function "yourself" with the combinations of times and boundary conditions.

If you can please edit above code, I would be thankful.

Best Regards. Qaisar

[Walter Roberson]

%%%%%  Caller Code

x0 = [1 ; 2];
tspan = 0:0.1:10
[t, x] = ode45(@(t,x) odefun(t,x), tspan, x0);
num_t = length(t);
e = zeros(num_t, 1);
for K = 1 : num_t
    [~, e(K)] = odefun(t(K), x(K,:));
end


%%%%%  ODE Function
function [dx, e] = odefun(t,x)

    if t<3
        e = x(1)-x(2) + 2;   % variable of my intrest, to be stored or returned each 0.1 sec exactly 
    
    elseif t<6
        e = x(1)-x(2) - sin(t);
    
    elseif t<10
        e = x(1)-x(2)/x(1);
    end
    
    dx(1) = x(2)
    dx(2) = sin(2*pi*t)

end

However, your tspan runs to t = 10 but you do not define e for t = 10.

[Walter Roberson]

Note:

It is a usage error to have use two different branches of a time-based if statement in the same ode45() call, unless the first and second derivatives of the left and right side of the equations match.

You need to break up your tspan so that you run multiple ode45() calls, one for 0 to 3, one for 3 to 6, one for 6 to 10.

The mathematics of ode45() is only valid if the function has continuity of first and second derivatives, but in most cases if you have an if then you have discontinuous derivatives.

[Muhammad Qaisar Ali]

Thank you very much from the core of my heart.

Regards.

Qaisar

[Muhammad Qaisar Ali]

%%%%%  Caller Code
clear all;

x0 = [1 2];
tspan = 0:0.1:10;
[t, x] = ode45(@(t,x) odefun(t,x), tspan, x0);

figure;
plot(t,x)
legend('x1','x2')
grid on;

num_t = length(t);
E = zeros(num_t, 1);
for K = 1 : num_t
    [~, E(K)] = odefun(t(K), x(K,:));
end

Output argument "e" (and possibly others) not assigned a value in the execution with "solution>odefun" function.

size(e)

figure;
plot(t,e)


%%%%%  ODE Function

function [dx, e] = odefun(t,x)
    dx = zeros(2,1);

    if t<3
        e = x(1)-x(2) + 2;   % variable of my intrest, to be stored or returned each 0.1 sec exactly 
    
    elseif t<6
        e = x(1)-x(2) - sin(t);
    
    elseif t<10
        e = x(1)-x(2)/x(1);
    end
    
    dx(1) = x(2);
    dx(2) = -2*(x(1) - 2) - (x(2) );


end

Now getting this error.

[Walter Roberson]

Read your code. What happens when t becomes 10 exactly ?

[Muhammad Qaisar Ali]

oh, sorry my bad.

got the point. when t=10, there is no true condition and hence, no thing is assigned to e.

Thank you Sir..