How can I save the values inside ODE solver?

12 次查看（过去 30 天）
Qiming
Qiming 2012-12-24
评论： Matteo Giovanni Perego 2018-1-30
Hi:
I have a question using ODE45. I generate a "random" value inside the my ODE function. I also need this random value in my main function. How can I save these random values? I cannot regenerate those random values after calling ODE45 since they will change.
Thanks for any help.
Paler

请先登录，再回答此问题。

回答（3 个）

Jan
Jan 2012-12-24
If the random value influences the integration, ODE45 is driven outside its specifications: The function to be integrated must be smooth, otherwise the stepsize control fails and ridiculously small stepsizes will at first slow down the calculations dramatically and at second the accumulated roundoff errors will dominate the result of the integration. Therefore these commands are massive DONTs in the integrand: IF, RAND, NOW (or any other time related functions), and even ABS is "evil", when the zero limit is crossed.
If you want to integrate a function, which has random parts, use a fixed step solver.
  1 个评论
Matteo Giovanni Perego
Matteo Giovanni Perego 2018-1-30
There's actually a solution to this problem: generating the random parts of the function only when we get a successful integration step.
Here's an example:
% Main:
y0 = 5; tspan = [0 100]; b = 10;
OutputFcn(0,0,'clear');
opts = odeset('OutputFcn',@OutputFcn,'refine',1);
[ts,ys] = ode45(@(t,y)odefun(t,y,b),tspan,y0,opts);
[~,rands] = OutputFcn(0,0,'allout');
figure(1)
subplot(2,1,1)
plot(ts,ys,'o-')
ylabel('y')
xlabel('x')
subplot(2,1,2)
plot(ts,rands,'ro-')
ylabel('random part')
xlabel('x')
% Functions:
function dy_dt = odefun(t,y,b)
% ODE:
[~,k] = OutputFcn(t,y,'pass2ode');
dy_dt = b*k;
end
function [status,out_rand] = OutputFcn(t,y,flag)
%OUTPUT FUNCTION:
persistent random
% Initialize random:
if isempty(random)
random = rand;
end
% Random number generation and output:
if isempty(flag)
% Successful integration step! Generate a new value:
random = [random; rand];
elseif strcmp(flag,'pass2ode')
out_rand = random(end);
elseif strcmp(flag,'allout')
out_rand = random;
elseif strcmp(flag,'clear')
clear random
end
% Always keep integrating:
status = 0;
end
Also, we need to set the 'refine' property to 1 since interpolated values between our integration steps would be complete nonsense due to the randomness of the increments.
Hope it helps!

请先登录，再进行评论。

Laura Proctor
Laura Proctor 2012-12-24
编辑：Laura Proctor 2012-12-24
Does the random value need to be generated inside your ODE function? In other words, can you pass it in as an input variable?
For example, suppose that your ODE function will now have the declaration:
function y = odefun(t,randValue)
Then, in your main function, you can generate the random value and create an anonymous function handle to embed it into the function so that you can call your ODE solver with this anonymous function handle:
randValue = rand;
afh = @(t) odefun(t,randValue);
[t,y] = ode45(afh,tspan,y0);
And now you have access to the random value inside your main function as the variable randValue.
  3 个评论
显示 1更早的评论
Laura Proctor
Laura Proctor 2012-12-24
That's an interesting problem... you could use a global variable. But, I hate to recommend it, because global variables are generally not the best solution. I hope that someone else can provide a more elegant solution!
Laura Proctor
Laura Proctor 2012-12-24
Another idea is to use assignin to assign the value to a variable in your main function.

请先登录，再进行评论。

Image Analyst
Image Analyst 2012-12-24
I don't know anything about ODE45, but if your question is about how to share variables among other functions (because normally variables are just local to the function in which they live), you can see the FAQ for several suggestions of how to share variables. http://matlab.wikia.com/wiki/FAQ#How_can_I_share_data_between_callback_functions_in_my_GUI.28s.29.3F

类别

在 Help Center 和 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