10

My professor would like me to solve a system similar to the following: $$ dx_i=[f_i(x_1,x_2,...x_n)]dt + g_ix_idW_i$$

Where $g_i$ are positive constants that measure the amplitude of the random perturbations, and $W_i$ are random variables normally distributed.

Im not sure how I can implement this in Matlab for ode45 to solve. What is throwing me off is the $dt$ tacked on to $f_i(x_1,...)$ and $dW_i$.

Is it as simple as coding 

    function dxdt=money(t,x,a,b,c)

x1=x(1);
x2=x(2);
x3=x(3);

dx1=x1.*(x2-a)+x3 + 10*rand();
dx2=1-b*x2-x1.^2 + 10*rand;
dx3=-x1-c*x3 +10*rand();

dxdt=[dx1;dx2;dx3];

end
Rodrigo de Azevedo
  • 23,223
Demetri Pananos
  • 3,294

3 Answers3

13

You absolutely cannot use ode45 to simulate this. ODEs and SDEs are different beasts. Also note that rand in Matlab is uniformly, not normally, distributed, so you're also not even generating proper Wiener increments. Before proceeding, you 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.

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 and anonymous functions:

a = ...
b = ...
c = ...
f = @(t,x)[ x(1)*(x(2)-a)+x(3);
            1-b*x(2)-x(2)^2;
           -x(1)-c*x(3)]; % drift function
g = 10; % or @(t,x)10;    % diffusion function, constant in this case
dt = 1e-3;                % time step
t = 0:dt:1;               % time vector
x0 = [...];               % 3-by-1 initial condition
opts = sdeset('RandSeed',1);  % Set random seed
y = sde_euler(f,g,t,x0,opts); % Integrate
figure;
plot(t,y);

In your example, the drift function, $g(t,x)$, doesn't depend on the state variable, $x$. This is so-called additive noise. In this case, more complex SDE integration schemes like Milstein reduce to Euler-Maruyma. More importantly, if you do use more complex diffusion functions (e.g., multiplicative noise where $g(t,x) = g_1(t)x$ like in the equation at the beginning of your question), you need to specify what SDE formulation you're using, Itô (most common in finance) or Stratonovich. My SDETools currently defaults to Stratonovich, but you can change it via opts = sdeset('SDEType','Ito');.

horchler
  • 3,258
  • 2
  • 27
  • 41
  • EDIT: Oh, I see, f is a vector. Thank you.

    This looks promising. I notice that your example is for 1 SDE. What about a system of SDEs. Can I do something similar?

     – Demetri Pananos Apr 08 '14 at 20:50
  • Is your edited comment just for the first three sentences? This is a system of SDEs in $x$. The output, y, of sde_euler is 3-by-length(t). Or are you asking about simultaneously simulating multiple instances of this system? Sure you can do that too, you just nee to write f and g. As with ODEs, a system of first order SDEs can be generalized to pretty much anything. – horchler Apr 08 '14 at 21:41
  • @horchler Would the approach be similar with Simulink? Could the Simulink ODE4 runge-kutta solver be used with a sufficiently small fixed time step? – m_power Jun 19 '19 at 19:00
  • @horchler I would like to know: Why can't I generate a pre-defined noise vector, e.g. noise=randn(1,1000) and then call ode45 with and interpolate the noise like so: ode45(@(t_,x_)fFun(t_,x_)+gFun(t_,x_)*interp1(linspace(0,tf,1000),noise,t_),tSpan,x0)? – link Jul 06 '23 at 13:35
  • @link While you can do this, real noise is not smoothed and differentiable like your interpolated case. You could be a bit closer if you used something like Brownian Bridge Interpolation. But fundamentally, ode45 is designed to solve ODEs, not SDEs, and therefore makes certain assumptions on smoothness of the underlying solution. If you're interested in a rigorous and statistically correct solution that captures the true distribution of the noise, you'd need to derive an SDE formulation and solve it with an SDE solver. – horchler Jul 06 '23 at 18:16
  • So, if you just want to perturb a system with some arbitrary signal your approach might serve a purpose. But if you want to model and simulate the actual statistical behavior of your system with respect to defined noise you'll want to look to SDEs. – horchler Jul 06 '23 at 18:18
  • That makes sense. It also yields vastly different trajectories (with same seed) when I compare the processes using your sde MATLAB functions (which are well-documented btw). Also, using ode45 and a different number of random numbers, say randn(1,100), also Leads to vastly different results compared to randn(1,1000) which also makes sense. Thank you for the quick response – link Jul 08 '23 at 06:50
  • You’re welcome. Maybe more important is what the ensemble distribution of behavior of ~1000 runs at different seeds looks like, as well as what your noise is attempting to model. – horchler Jul 09 '23 at 21:15
2

Unfortunately, I don't think you can use ode45 to solve it, whereas, since $dW$ distributes like $N(0,\sqrt{dt})$ you must use the following discretisation: $$ x_{n+1,i}=x_{n,i}+f_i(x_{n,i})\Delta t+g_i(x_{n,i})x_{n,i}\sqrt{\Delta t}\cdot\epsilon_i, $$ $$ \epsilon_i=\mbox{randn()} $$ $$ x_{n,i}\approx x_i(t) $$

7raiden7
  • 1,904
  • I.e., Euler-Maruyama (assuming that the Itô formulation is desired). You might clarify that rand returns a normal variate from $N(0,1)$ to not confuse people with the Matlab function that return uniform variates (or maybe just $\epsilon_i \sim N(0,1)$). – horchler Apr 08 '14 at 18:21
  • Absolutely right, it actually is "randn", I will correct the answer immediately. – 7raiden7 Apr 09 '14 at 07:48
1

Does an exact solution exist for your SDE? If not, you must use a numerical method.

What 7Raiden7 suggests is called the Euler-Maruyama method and is based on the truncated Ito-Taylor expansion.

Another simple numerical method would be the Milstein Scheme, which contains additional terms from the Ito-Taylor expansion. Its implementation is easy to program in MATLAB and exhibits a higher order of convergence than the Euler-Maruyama method.

These will not solve your system, but they might get you started.

Bryan Johnson
  • 111