Solving second order SDE with Gaussian white noise for first time derivative in Matlab

Question

I'm having trouble solving a second order differential equation with Gaussian white noise. The equation I'm solving follows the form:

$$Ax'' + Bx' + \sin(x) = i + i_{n}$$

where $i_{n}$ is the Gaussian white noise, and $i$ is a constant source. I'm ultimately interested in $x'$, not $x$. Without the Gaussian source, I've been able to simply use an ODE solver like ode15s or ode45. However, with the random source, I don't know how to use an SDE solver to get $x'$. Can anyone with experience in solving SDEs in MATLAB lend some advice?

Answer 1

The most straightforward way to solve your SDE is with an Euler-Maruyama scheme. This is a simple and effective method for additive noise, i.e., the diffusion/noise term is not a function of the state, as appears to be the case for your example. Here is some Matlab code to solve your system:

n = 2;        % Order of system
t0 = 0;       % Initial time
dt = 1e-3;    % Fixed time step
tf = 1;       % Final time
t = t0:dt:tf; % Time vector
tlen = length(t);

A = 1;
B = 1;
c = 1;
f = @(t,x)[x(2);(-B*x(2)-sin(x(1))+c)/A]; % Drift function
ep = 1e-1;                                % Size of additive noise
g = @(t,x)[0;ep];                         % Diffusion function

x = zeros(lt,n); % Allocate output
x0 = [1;1];
x(1,:) = x0;     % Set initial condition

seed = 1;  % Seed value
rng(seed); % Always seed random number generator

% Euler-Maruyama
for i = 1:tlen-1
    x(i+1,:) = x(i,:).' + f(t(i),x(i,:))*dt + g(t(i),x(i,:)).*randn(n,1)*sqrt(dt);
end

figure;
plot(t,x);
xlabel('Time');
ylabel('State');
legend('x(t)','xdot(t)')

This is not an adaptive method like ode45, so you need to ensure that your integration step size, dt, is sufficiently small. There are many ways to optimize the above code and make it less general (e.g., simplify anonymous functions, pre-calculate normal variates, etc.). You could also try using sde_euler in my SDETools Matlab toolbox on Github, which has many options akin to those in the ODE suite. See this answer on Math.SE.

For methods specialized to second-order SDEs like yours, you could check out this paper (PDF) by Burrage, et al. For an introduction to numerically solving SDEs, 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. Also, note that random number generator seeding used in the code is now deprecated.