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.
采纳的回答
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?
6 个评论
elena
2019-5-28
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
2019-5-28
Akash Singh
2020-4-14
Hey @Jan
Shouldn't it be ?
p = sin(y(2,:)); % <- (Instead of y(:,2))
dy = [2 .* t; p ./ y(1,:)];
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
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Translated by 
