# Passing time step of ODE solver to odefunction for stochastic simulation

24 vues (au cours des 30 derniers jours)
soeren graabaek le 30 Mar 2022
Hello,
I need help passing the step length of the ODE solve to my ODE function at simulation time. I have found two other questions dicussing the issue, but I have been unable to solve the problem.
My ODE function need the actual step length at runtime to simulate brownian motion in my dynamic system. The usual approach is to use fixed step length solves, however, my simulation would benefit from using e.g. a ODE45 solve.
------------------------------------------------------------------------------ MATWORKS ------------------------------------------------------------------------------
Step length discussion:
Stochastic Differential Equation MATLAB
Fixed step length solves
An SDE and ODE45 explanation
------------------------------------------------------------------------------ APPENDIX ------------------------------------------------------------------------------
[...]
[Source: "Simulation and Inference for Stochastic Differential Equations" by Stefano M. Iacus]
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Réponse acceptée

Riccardo Scorretti le 30 Mar 2022
Hi Soeren,
in my opinion the question is ill-posed, because (up to my knowledge) methods like ode45 are designed to solve deterministic equations. In particular, they adapt the time step so as to achive the required accuracy, and I am seriously worried about the effect of the the algorithm used to compute the optimal time step on the final result. I think that you should not use them.
However, let's try to do if tor the fun. Assume, for simplicity, that you want to solve the simplified stocastic ODE (that is, with and ):
where . When ode45 or any other similar ODE-solver is used, it is mandatory to provide the derivative of X, which depends on the derivative of W. Notice that the derivative of W should write:
Here we come to the question you asked, which is to get . This should be possible by diverting the ODE option OutputFcn. According to Matlab documentation, "The ODE solver calls the output function after each successful time step". Hence, this function is aware of the last successfull time step, which can be stored (and returned) by using a persistent variable.
However, I would like to stress that in my opinion this is ot the right approach. I think that the right approach is to desing a method like the one reported in the text provided in your own question.
I implemented both methods in the code hereafter. At a first glance, the two methods seem to provide similar results; of course, the two solutions are bound to be point-to-point different, because we are solving a stochastic equation. A more detailed study of the statistical properties of the obtained results should be performed.
Moreover, one sees that in this case ode45 is much slower than the hand-written forward Euler algorithm, hence I really don't see any reason to make use of it.
I hope it helps.
sigma = 1; % A more complex solution can be found
T_end = 1;
opts = odeset('OutputFcn', @myOutputFn);
myOutputFn(0, [], 'reset');
tic
[t, y] = ode45(@(t,y) fun(t,y,sigma), [0, T_end], 0, opts);
toc
Elapsed time is 24.448954 seconds.
figure ; plot(t, y) ; hold on
tic
nb = 10000;
t = linspace(0, T_end, nb) ; dt = T_end/(nb-1);
y = zeros(nb, 1);
w = zeros(nb, 1);
y(1) = 0;
w(1) = 0;
for k = 2 : nb
% The following lines are equivalent
% You can use this:
w(k) = w(k-1) + sqrt(dt)*randn(1);
y(k) = y(k-1) + sigma*(w(k)-w(k-1));
% Or this:
% y(k) = y(k-1) + sigma*sqrt(dt)*randn(1);
end
toc
Elapsed time is 0.016692 seconds.
plot(t, y);
legend('ODE45', 'Euler');
return
function dydt = fun(t, y, sigma)
delta_t = myOutputFn(t, [], 'delta_t');
assert(delta_t >= 0); % Of course, this test can be removed
if delta_t == 0
w = 0;
else
w = randn(1)/sqrt(delta_t);
end
dydt = sigma*w;
end
function status = myOutputFn(t, y, flag)
persistent t0
if isempty(t0) , t0 = 0 ; end
switch(flag)
case 'reset'
t0 = t(1);
case 'init'
t0 = t(1);
case 'delta_t'
status = t-t0;
case ''
t0 = t(1);
status = 0;
case 'done'
% Nothing to do
otherwise
error('*** This cannot happen!!! ***');
end
end
##### 5 commentairesAfficher 3 commentaires plus anciensMasquer 3 commentaires plus anciens
Torsten le 31 Mar 2022
Thank you, Riccardo.
soeren graabaek le 1 Avr 2022
Hi Torsten and Riccardo,

Connectez-vous pour commenter.

### Plus de réponses (1)

Torsten le 30 Mar 2022
My ODE function need the actual step length at runtime to simulate brownian motion in my dynamic system. The usual approach is to use fixed step length solves, however, my simulation would benefit from using e.g. a ODE45 solve.
Because MATLAB's ODE solvers adapt stepsize during integration in order to satisfy an error tolerance, they are only suited to be used for deterministic ODEs. The step size control becomes mad if you introduce stochastic inputs to the ODEs. You will have to stick to Euler explicit with fixed stepsize in time.
##### 0 commentairesAfficher -2 commentaires plus anciensMasquer -2 commentaires plus anciens

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by