Passing time step from ODE solver (ode23tb) to ODE function
显示 更早的评论
I am currently utilizing MATLAB's ODE solver ode23tb to advance a system of equations with the structure:
[t, X] = ode23tb(@myfun, [t1 t2], X0)
With the ODE system defined as:
function dxdt = myfun(t, x)
% Some operations to figure out the time rate of x based on current values of x, stored ultimately in a variable 'blah'
dxdt = blah;
end
Here is the catch. The evaluation of 'blah' depends on all previous steps taken. (My governing equation involves fractional derivative, which fundamentally has a convolution integral ... without getting into too much detail.) Therefore I need to find an array of the successful steps taken by ode23tb.
So far, the only way I have been able to make this work is to create a modified version of ode23tb and globalize tout. (I also had to place a copy of the private function folder in the local directory.) By checking tout inside myfun, I can get a list of the steps that ode23tb stored as time history of the ODE advancement.
However, declaring global variable is obviously a bad practice. I also need to do parameter sweep eventually and run myfun multiple times with different input parameters. Yet, the global declaration creates a mess that can only be cleaned with the horrible clear all command -- which makes it challenging to run parameter sweep with a loop structure.
What would be the proper way to retrieve the tout information from ode23tb without generating a monster comprised of global and clear all?
Thanks in advance for the useful feedbacks!
采纳的回答
Jan
2021-2-27
You can use the OutputFcn to store the successful steps in a persitent variable, which can be forwarded to the function to bbe integrated. This is ways faster than writing to disk.
function status = myOutputFcn(t,y,flag)
persistent T_local, T_count
switch char(flag)
case ''
T_count = T_count + 1;
if T_count > numel(T_local)
T_local(T_count + 1000) = 0;
end
T_local(t_count) = t;
status = 0;
case {'init', 'done'}
T_local = zeros(1, 1000);
T_count = 0;
status = 0;
case 'get'
status = T_local(1:T_count);
end
end
Now the funciton to be integrated cann call:
formerT = myOutputFcn([], [], 'get')
更多回答（1 个）
Steven Lord
2021-2-27
2 个投票
The evaluation of 'blah' depends on all previous steps taken.
That sounds like you don't have an ordinary differential equation but a delay differential equation. In that case look at dde23.
类别
在 帮助中心 和 File Exchange 中查找有关 Historical Contests 的更多信息
Translated by 
