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.
采纳的回答
Walter Roberson
2022-9-18
1 个投票
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.
7 个评论
Muhammad Qaisar Ali
2022-9-18
编辑：Muhammad Qaisar Ali
2022-9-18
Walter Roberson
2022-9-18
%%%%% 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
2022-9-18
编辑：Walter Roberson
2022-9-18
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
2022-9-18
%%%%% 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
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
2022-9-18
Read your code. What happens when t becomes 10 exactly ?
Muhammad Qaisar Ali
2022-9-18
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
