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 
):
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:
. 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.
. 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
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
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
