1
$\begingroup$

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?

$\endgroup$
  • $\begingroup$ What do you mean by "a constant source"? Is $i$ just a constant parameter or is it noise with a uniform distribution? Or is the right-hand side an Ornstein-Uhlenbeck process? $\endgroup$ – horchler Jan 27 '16 at 1:17
3
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Thank you for your help! I got the code to work, but get the same result each time I run it. I'm ultimately hoping to generate a histogram of a certain x' threshold value, and thought the white noise would allow for this variation. Is there a reason that the code would converge to the same solution every time it's run, despite the additive noise? Thanks! $\endgroup$ – Emily Jan 27 '16 at 21:18
  • $\begingroup$ Yes, a random seed is specified so the same sequence of pseudorandom variates will be produced each time. This allows you to re-run and replicate your results exactly. If you want $N$ independent runs, you either need to set the seed once at the very top of your code outside of any loop that runs the simulation $N$ times, or vectorize the integration such that all $N$ simulations are calculated simultaneously, or you could re-seed using a different value for each run (slightly less efficient). $\endgroup$ – horchler Jan 27 '16 at 21:49

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.