Euler–Maruyama method

<h2 id="example">Example</h2>
<h3>Numerical simulation</h3>

<p>An area that has benefited significantly from SDEs is <a href="/facts/Mathematical_biology/pV9KCBSD">mathematical biology</a>. As many biological processes are both stochastic and continuous in nature, numerical methods of solving SDEs are highly valuable in the field.
</p><p>The graphic depicts a stochastic differential equation solved using the Euler-Maruyama method. The deterministic counterpart is shown in blue.
</p>
<h3>Computer implementation</h3>
<p>The following <a href="/facts/Python_(programming_language)/YbuGqofa">Python</a> code implements the Euler–Maruyama method and uses it to solve the <a href="/facts/Ornstein%E2%80%93Uhlenbeck_process/oFaDcovm">Ornstein–Uhlenbeck process</a> defined by
</p>

d
        
          Y
          
            t
          
        
        =
        θ
        ⋅
        (
        μ
        −
        
          Y
          
            t
          
        
        )
        
        
          
            d
          
        
        t
        +
        σ
        
        
          
            d
          
        
        
          W
          
            t
          
        
      
    
    {\displaystyle dY_{t}=\theta \cdot (\mu -Y_{t})\,{\mathrm {d} }t+\sigma \,{\mathrm {d} }W_{t}}

Y
          
            0
          
        
        =
        
          Y
          
            
              i
              n
              i
              t
            
          
        
        .
      
    
    {\displaystyle Y_{0}=Y_{\mathrm {init} }.}

<p>The random numbers for 
  
    
      
        
          
            d
          
        
        
          W
          
            t
          
        
      
    
    {\displaystyle {\mathrm {d} }W_{t}}
  
 are generated using the <a href="/facts/NumPy/T6FnhWWD">NumPy</a> mathematics package.
</p>
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt

class Model:
    """Stochastic model constants."""
    THETA = 0.7
    MU = 1.5
    SIGMA = 0.06

def mu(y: float, _t: float) -> float:
    """Implement the Ornstein–Uhlenbeck mu."""
    return Model.THETA * (Model.MU - y)

def sigma(_y: float, _t: float) -> float:
    """Implement the Ornstein–Uhlenbeck sigma."""
    return Model.SIGMA

def dW(delta_t: float) -> float:
    """Sample a random number at each call."""
    return np.random.normal(loc=0.0, scale=np.sqrt(delta_t))

def run_simulation():
    """ Return the result of one full simulation."""
    T_INIT = 3
    T_END = 7
    N = 1000  # Compute at 1000 grid points
    DT = float(T_END - T_INIT) / N
    TS = np.arange(T_INIT, T_END + DT, DT)
    assert TS.size == N + 1

Y_INIT = 0

ys = np.zeros(TS.size)
    ys[0] = Y_INIT
    for i in range(1, TS.size):
        t = T_INIT + (i - 1) * DT
        y = ys[i - 1]
        ys[i] = y + mu(y, t) * DT + sigma(y, t) * dW(DT)

return TS, ys

def plot_simulations(num_sims: int):
    """ Plot several simulations in one image."""
    for _ in range(num_sims):
        plt.plot(*run_simulation())

plt.xlabel("time")
    plt.ylabel("y")
    plt.show()

if __name__ == "__main__":
    NUM_SIMS = 5
    plot_simulations(NUM_SIMS)

<p>The following is simply the translation of the above code into the <a href="/facts/MATLAB/qPjLISCk">MATLAB (R2019b)</a> programming language:
</p>
%% Initialization and Utility
close all;
clear all;

numSims = 5;            % display five runs
tBounds = [3 7];        % The bounds of t
N       = 1000;         % Compute at 1000 grid points
dt      = (tBounds(2) - tBounds(1)) / N ;
y_init  = 0;            % Initial y condition

% Initialize the probability distribution for our
% random variable with mean 0 and stdev of sqrt(dt)
pd = makedist('Normal',0,sqrt(dt));

c = [0.7, 1.5, 0.06];   % Theta, Mu, and Sigma, respectively

ts = linspace(tBounds(1), tBounds(2), N); % From t0-->t1 with N points
ys = zeros(1,N);     % 1xN Matrix of zeros

ys(1) = y_init;
%% Computing the Process
for j = 1:numSims
    for i = 2:numel(ts)
        t = tBounds(1) + (i-1) .* dt;
        y = ys(i-1);
        mu      = c(1) .* (c(2) - y);
        sigma   = c(3);
        dW      = random(pd);
        
        ys(i) = y + mu .* dt + sigma .* dW;
    end
    figure;
    hold on;
    plot(ts, ys, 'o')
end

<h2 id="see-also">See also</h2>
<ul><li><a href="/facts/Milstein_method/YlnTHt5J">Milstein method</a></li>
<li><a href="/facts/Runge%E2%80%93Kutta_method_(SDE)/zm8G0I5N">Runge–Kutta method (SDE)</a></li>
<li>Leimkuhler–Matthews method</li></ul>

<h2 id="references">References</h2>

<ol>
<li id="fn:1"><p>Kloeden, P.E. & Platen, E. (1992). Numerical Solution of Stochastic Differential Equations. Springer, Berlin. ISBN 3-540-54062-8. <a href="3-540-54062-8" target="_blank">3-540-54062-8</a> <a href="#fnref:1" class="footnote-back-ref">↩</a></p></li>
</ol>

Euler–Maruyama method open-in-new

Euler–Maruyama method