9

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
  • 20,693
  • 5
  • 43
  • 100
Demetri Pananos
  • 2,962
  • 2
  • 23
  • 35

3 Answers3

11

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,108
  • 2
  • 24
  • 39
  • 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
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,744
  • 10
  • 9
  • 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
  • 2