How can I save the values inside ODE solver?

조회 수: 11 (최근 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
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 CenterFile 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