ode45 for Langevin equation

Question

I have a question about the use of Matlab to compute solution of stochastic differentials equations. The equations are the 2.2a,b, page 3, in this paper (PDF).

My professor suggested using ode45 with a small time step, but the results do not match with those in the article. In particular the time series and the pdf. I also have a doubt about the definition of the white noise in the function.

Here the code for the integration function:

function dVdt = R_Lang( t,V )

global  sigma lambda alpha

W1=sigma*randn(1,1); 
W2=sigma*randn(1,1);
dVdt=[alpha*V(1)+lambda*V(1)^3+1/V(1)*0.5*sigma^2+W1;
        sigma/V(1)*W2];

end

Main script:

clear variables
close all
global  sigma lambda alpha
sigma=sqrt(2*0.0028);
alpha=3.81;
lambda=-5604;

tspan=[0,10];
options = odeset('RelTol',1E-6,'AbsTol',1E-6,'MaxStep',0.05);

A0=random('norm',0,0.5,[2,1]);
[t,L]=ode45(@(t,L) R_Lang(t,L),tspan,A0,options);

If you have any suggestions I'd be grateful.

Here the new code to confront my EM method and 'sde_euler'.

lambda = -5604; 
sigma=sqrt(2*0.0028) ; 
Rzero = 0.03;    % problem parameters
phizero=-1;
dt=1e-5;
T = 0:dt:10;
N=length(T);
Xi1 = sigma*randn(1,N);         % Gaussian Noise with variance=sigma^2
 Xi2 = sigma*randn(1,N);
alpha=3.81;

Rem = zeros(1,N);                 % preallocate for efficiency
Rtemp = Rzero;
phiem = zeros(1,N);                 % preallocate for efficiency
phitemp = phizero;
for j = 1:N
     Rtemp = Rtemp + dt*(alpha*Rtemp+lambda*Rtemp^3+sigma^2/(2*Rtemp)) + sigma*Xi1(j);
   phitemp=phitemp+sigma/Rtemp*Xi2(j);
   phiem(j)=phitemp;
   Rem(j) = Rtemp;

end

f = @(t,V)[alpha*V(1)+lambda*V(1)^3+0.5*sigma^2/V(1)/2;
           0];                 % Drift function
g = @(t,V)[sigma;
           sigma/V(1)];        % Diffusion function
A0 = [0.03;0];                % 2-by-1 initial condition
opts = sdeset('RandSeed',1,'SDEType','Ito'); % Set random seed, use Ito formulation
L = sde_euler(f,g,T,A0,opts);       

plot(T,Rem,'r')
hold on
plot(T,L(:,1),'b')

Thanks again for the help !

Answer 1

ODEs and SDEs are very different and one should not use tools for ODEs, like ode45, to try to solve SDEs. Looking at the paper you linked to, they used a basic Euler-Maruyama scheme to integrate the system. This a very simple solver to implement yourself.

Before proceeding, you (and your professor!) should take some time to read up on SDEs and how to solve them numerically. I recommend this paper, which includes many Matlab examples:

Desmond J. Higham, 2001, An Algorithmic Introduction to Numerical Simulation of Stochastic Differential Equations, SIAM Rev. (Educ. Sect.), 43 525–46. http://dx.doi.org/10.1137/S0036144500378302

The URL to the Matlab files in the paper won't work; use this one. Note, that as this a 15-year old paper, some of the code related to random number generation is out of date (use rng(1) instead of randn('state',1) to seed the generator).

If you are familiar with ode45 you might look at my SDETools Matlab toolbox on GitHub. It was designed to be fast and has an interface that works very similarly to Matlab's ODE suite. Here is how you might code up your example using the Euler-Maruyma solver:

sigma = 1e-1*sqrt(2*0.0028);
lambda = -5604;
alpha = 3.81;
f = @(t,V)[alpha*V(1)+lambda*V(1)^3+0.5*sigma^2/V(1);
           0];                 % Drift function
g = @(t,V)[sigma;
           sigma/V(1)];        % Diffusion function
dt = 1e-3;                     % Time step
t = 0:dt:10;                   % Time vector
A0 = [0.03;-2];                % 2-by-1 initial condition
opts = sdeset('RandSeed',1,'SDEType','Ito'); % Set random seed, use Ito formulation
L = sde_euler(f,g,t,A0,opts);                % Integrate

figure;
subplot(211);
plot(t,L(:,2));
ylabel('\phi');
subplot(212);
plot(t,L(:,1));
ylabel('r');
xlabel('t');

I had to reduce the size of sigma or the noise was so large that it could cause the radius variable to go negative. I'm not sure if the paper discusses how they handle this singularity. You can try the 'NonNegative' option within sdeset to try to handle this or you may need to construct your own solver. I also couldn't find what integration time step the paper used. You should also consider contacting the authors of the paper directly.

UPDATE

Here's an Euler-Maruyama implementation that matches the sde_euler code above:

sigma = 1e-1*sqrt(2*0.0028);
lambda = -5604;
alpha = 3.81;
f = @(t,V)[alpha*V(1)+lambda*V(1)^3+0.5*sigma^2/V(1);
           0];                 % Drift function
g = @(t,V)[sigma;
           sigma/V(1)];        % Diffusion function
dt = 1e-3;                     % Time step
t = 0:dt:10;                   % Time vector
A0 = [0.03;-2];                % 2-by-1 initial condition

% Create and initialize state vector (L here is transposed relative to sde_euler output)
lt = length(t);
n = length(A0);
L = zeros(n,lt);
L(:,1) = A0;

% Set seed and pre-calculate Wiener increments with order matching sde_euler
rng(1);
r = sqrt(dt)*randn(lt-1,n).';

% General Euler-Maruyama integration loop
for i = 1:lt-1
    L(:,i+1) = L(:,i)+f(t(i),L(:,i))*dt+r(:,i).*g(t(i),L(:,i));
end

figure;
subplot(211);
plot(t,L(2,:));
ylabel('\phi');
subplot(212);
plot(t,L(1,:));
ylabel('r');
xlabel('t');