Euler–Maruyama method

Consider the stochastic differential equation (see Itô calculus)

d X_{t} = a (X_{t}, t) d t + b (X_{t}, t) d W_{t},

with initial condition X₀ = x₀, where W_t denotes the Wiener process, and suppose that we wish to solve this SDE on some interval of time [0, T]. Then the Euler–Maruyama approximation to the true solution X is the Markov chain Y defined as follows:

Partition the interval [0, T] into N equal subintervals of width $Δ t > 0$ :

0 = τ_{0} < τ_{1} < \dots < τ_{N} = T and Δ t = T / N;

Set Y₀ = x₀
Recursively define Y_n for 0 ≤ n ≤ N-1 by

Y_{n + 1} = Y_{n} + a (Y_{n}, τ_{n}) Δ t + b (Y_{n}, τ_{n}) Δ W_{n},

where

Δ W_{n} = W_{τ_{n + 1}} - W_{τ_{n}} .

The random variables ΔW_n are independent and identically distributed normal random variables with expected value zero and variance Δt.

Example

Numerical simulation

An area that has benefited significantly from SDEs is mathematical biology. As many biological processes are both stochastic and continuous in nature, numerical methods of solving SDEs are highly valuable in the field.

The graphic depicts a stochastic differential equation solved using the Euler-Maruyama method. The deterministic counterpart is shown in blue.

Computer implementation

The following Python code implements the Euler–Maruyama method and uses it to solve the Ornstein–Uhlenbeck process defined by

d Y_{t} = θ \cdot (μ - Y_{t}) d t + σ d W_{t}

Y_{0} = Y_{i n i t} .

The random numbers for $d W_{t}$ are generated using the NumPy mathematics package.

<span id="mwdg"></span><span class="linenos" data-line="1" id="mwdw"></span><span class="c1" id="mweA"># -*- coding: utf-8 -*-</span> <span class="linenos" data-line="2" id="mweQ"></span><span class="kn" id="mweg">import</span> <span class="nn" id="mwew">numpy</span> <span class="k" id="mwfA">as</span> <span class="nn" id="mwfQ">np</span> <span class="linenos" data-line="3" id="mwfg"></span><span class="kn" id="mwfw">import</span> <span class="nn" id="mwgA">matplotlib.pyplot</span> <span class="k" id="mwgQ">as</span> <span class="nn" id="mwgg">plt</span> <span class="linenos" data-line="4" id="mwgw"></span> <span class="linenos" data-line="5" id="mwhA"></span> <span class="linenos" data-line="6" id="mwhQ"></span><span class="k" id="mwhg">class</span> <span class="nc" id="mwhw">Model</span><span class="p" id="mwiA">:</span> <span class="linenos" data-line="7" id="mwiQ"></span><span class="w" id="mwig"> </span><span class="sd" id="mwiw">"""Stochastic model constants."""</span> <span class="linenos" data-line="8" id="mwjA"></span> <span class="n" id="mwjQ">THETA</span> <span class="o" id="mwjg">=</span> <span class="mf" id="mwjw">0.7</span> <span class="linenos" data-line="9" id="mwkA"></span> <span class="n" id="mwkQ">MU</span> <span class="o" id="mwkg">=</span> <span class="mf" id="mwkw">1.5</span> <span class="linenos" data-line="10" id="mwlA"></span> <span class="n" id="mwlQ">SIGMA</span> <span class="o" id="mwlg">=</span> <span class="mf" id="mwlw">0.06</span> <span class="linenos" data-line="11" id="mwmA"></span> <span class="linenos" data-line="12" id="mwmQ"></span> <span class="linenos" data-line="13" id="mwmg"></span><span class="k" id="mwmw">def</span> <span class="nf" id="mwnA">mu</span><span class="p" id="mwnQ">(</span><span class="n" id="mwng">y</span><span class="p" id="mwnw">:</span> <span class="nb" id="mwoA">float</span><span class="p" id="mwoQ">,</span> <span class="n" id="mwog">_t</span><span class="p" id="mwow">:</span> <span class="nb" id="mwpA">float</span><span class="p" id="mwpQ">)</span> <span class="o" id="mwpg">-></span> <span class="nb" id="mwpw">float</span><span class="p" id="mwqA">:</span> <span class="linenos" data-line="14" id="mwqQ"></span><span class="w" id="mwqg"> </span><span class="sd" id="mwqw">"""Implement the Ornstein–Uhlenbeck mu."""</span> <span class="linenos" data-line="15" id="mwrA"></span> <span class="k" id="mwrQ">return</span> <span class="n" id="mwrg">Model</span><span class="o" id="mwrw">.</span><span class="n" id="mwsA">THETA</span> <span class="o" id="mwsQ">*</span> <span class="p" id="mwsg">(</span><span class="n" id="mwsw">Model</span><span class="o" id="mwtA">.</span><span class="n" id="mwtQ">MU</span> <span class="o" id="mwtg">-</span> <span class="n" id="mwtw">y</span><span class="p" id="mwuA">)</span> <span class="linenos" data-line="16" id="mwuQ"></span> <span class="linenos" data-line="17" id="mwug"></span> <span class="linenos" data-line="18" id="mwuw"></span><span class="k" id="mwvA">def</span> <span class="nf" id="mwvQ">sigma</span><span class="p" id="mwvg">(</span><span class="n" id="mwvw">_y</span><span class="p" id="mwwA">:</span> <span class="nb" id="mwwQ">float</span><span class="p" id="mwwg">,</span> <span class="n" id="mwww">_t</span><span class="p" id="mwxA">:</span> <span class="nb" id="mwxQ">float</span><span class="p" id="mwxg">)</span> <span class="o" id="mwxw">-></span> <span class="nb" id="mwyA">float</span><span class="p" id="mwyQ">:</span> <span class="linenos" data-line="19" id="mwyg"></span><span class="w" id="mwyw"> </span><span class="sd" id="mwzA">"""Implement the Ornstein–Uhlenbeck sigma."""</span> <span class="linenos" data-line="20" id="mwzQ"></span> <span class="k" id="mwzg">return</span> <span class="n" id="mwzw">Model</span><span class="o" id="mw0A">.</span><span class="n" id="mw0Q">SIGMA</span> <span class="linenos" data-line="21" id="mw0g"></span> <span class="linenos" data-line="22" id="mw0w"></span> <span class="linenos" data-line="23" id="mw1A"></span><span class="k" id="mw1Q">def</span> <span class="nf" id="mw1g">dW</span><span class="p" id="mw1w">(</span><span class="n" id="mw2A">delta_t</span><span class="p" id="mw2Q">:</span> <span class="nb" id="mw2g">float</span><span class="p" id="mw2w">)</span> <span class="o" id="mw3A">-></span> <span class="nb" id="mw3Q">float</span><span class="p" id="mw3g">:</span> <span class="linenos" data-line="24" id="mw3w"></span><span class="w" id="mw4A"> </span><span class="sd" id="mw4Q">"""Sample a random number at each call."""</span> <span class="linenos" data-line="25" id="mw4g"></span> <span class="k" id="mw4w">return</span> <span class="n" id="mw5A">np</span><span class="o" id="mw5Q">.</span><span class="n" id="mw5g">random</span><span class="o" id="mw5w">.</span><span class="n" id="mw6A">normal</span><span class="p" id="mw6Q">(</span><span class="n" id="mw6g">loc</span><span class="o" id="mw6w">=</span><span class="mf" id="mw7A">0.0</span><span class="p" id="mw7Q">,</span> <span class="n" id="mw7g">scale</span><span class="o" id="mw7w">=</span><span class="n" id="mw8A">np</span><span class="o" id="mw8Q">.</span><span class="kp" id="mw8g">sqrt</span><span class="p" id="mw8w">(</span><span class="n" id="mw9A">delta_t</span><span class="p" id="mw9Q">))</span> <span class="linenos" data-line="26" id="mw9g"></span> <span class="linenos" data-line="27" id="mw9w"></span> <span class="linenos" data-line="28" id="mw-A"></span><span class="k" id="mw-Q">def</span> <span class="nf" id="mw-g">run_simulation</span><span class="p" id="mw-w">():</span> <span class="linenos" data-line="29" id="mw_A"></span><span class="w" id="mw_Q"> </span><span class="sd" id="mw_g">""" Return the result of one full simulation."""</span> <span class="linenos" data-line="30" id="mw_w"></span> <span class="n" id="mwAQA">T_INIT</span> <span class="o" id="mwAQE">=</span> <span class="mi" id="mwAQI">3</span> <span class="linenos" data-line="31" id="mwAQM"></span> <span class="n" id="mwAQQ">T_END</span> <span class="o" id="mwAQU">=</span> <span class="mi" id="mwAQY">7</span> <span class="linenos" data-line="32" id="mwAQc"></span> <span class="n" id="mwAQg">N</span> <span class="o" id="mwAQk">=</span> <span class="mi" id="mwAQo">1000</span> <span class="c1" id="mwAQs"># Compute at 1000 grid points</span> <span class="linenos" data-line="33" id="mwAQw"></span> <span class="n" id="mwAQ0">DT</span> <span class="o" id="mwAQ4">=</span> <span class="nb" id="mwAQ8">float</span><span class="p" id="mwARA">(</span><span class="n" id="mwARE">T_END</span> <span class="o" id="mwARI">-</span> <span class="n" id="mwARM">T_INIT</span><span class="p" id="mwARQ">)</span> <span class="o" id="mwARU">/</span> <span class="n" id="mwARY">N</span> <span class="linenos" data-line="34" id="mwARc"></span> <span class="n" id="mwARg">TS</span> <span class="o" id="mwARk">=</span> <span class="n" id="mwARo">np</span><span class="o" id="mwARs">.</span><span class="kp" id="mwARw">arange</span><span class="p" id="mwAR0">(</span><span class="n" id="mwAR4">T_INIT</span><span class="p" id="mwAR8">,</span> <span class="n" id="mwASA">T_END</span> <span class="o" id="mwASE">+</span> <span class="n" id="mwASI">DT</span><span class="p" id="mwASM">,</span> <span class="n" id="mwASQ">DT</span><span class="p" id="mwASU">)</span> <span class="linenos" data-line="35" id="mwASY"></span> <span class="k" id="mwASc">assert</span> <span class="n" id="mwASg">TS</span><span class="o" id="mwASk">.</span><span class="kp" id="mwASo">size</span> <span class="o" id="mwASs">==</span> <span class="n" id="mwASw">N</span> <span class="o" id="mwAS0">+</span> <span class="mi" id="mwAS4">1</span> <span class="linenos" data-line="36" id="mwAS8"></span> <span class="linenos" data-line="37" id="mwATA"></span> <span class="n" id="mwATE">Y_INIT</span> <span class="o" id="mwATI">=</span> <span class="mi" id="mwATM">0</span> <span class="linenos" data-line="38" id="mwATQ"></span> <span class="linenos" data-line="39" id="mwATU"></span> <span class="n" id="mwATY">ys</span> <span class="o" id="mwATc">=</span> <span class="n" id="mwATg">np</span><span class="o" id="mwATk">.</span><span class="kp" id="mwATo">zeros</span><span class="p" id="mwATs">(</span><span class="n" id="mwATw">TS</span><span class="o" id="mwAT0">.</span><span class="kp" id="mwAT4">size</span><span class="p" id="mwAT8">)</span> <span class="linenos" data-line="40" id="mwAUA"></span> <span class="n" id="mwAUE">ys</span><span class="p" id="mwAUI">[</span><span class="mi" id="mwAUM">0</span><span class="p" id="mwAUQ">]</span> <span class="o" id="mwAUU">=</span> <span class="n" id="mwAUY">Y_INIT</span> <span class="linenos" data-line="41" id="mwAUc"></span> <span class="k" id="mwAUg">for</span> <span class="n" id="mwAUk">i</span> <span class="ow" id="mwAUo">in</span> <span class="nb" id="mwAUs">range</span><span class="p" id="mwAUw">(</span><span class="mi" id="mwAU0">1</span><span class="p" id="mwAU4">,</span> <span class="n" id="mwAU8">TS</span><span class="o" id="mwAVA">.</span><span class="kp" id="mwAVE">size</span><span class="p" id="mwAVI">):</span> <span class="linenos" data-line="42" id="mwAVM"></span> <span class="n" id="mwAVQ">t</span> <span class="o" id="mwAVU">=</span> <span class="n" id="mwAVY">T_INIT</span> <span class="o" id="mwAVc">+</span> <span class="p" id="mwAVg">(</span><span class="n" id="mwAVk">i</span> <span class="o" id="mwAVo">-</span> <span class="mi" id="mwAVs">1</span><span class="p" id="mwAVw">)</span> <span class="o" id="mwAV0">*</span> <span class="n" id="mwAV4">DT</span> <span class="linenos" data-line="43" id="mwAV8"></span> <span class="n" id="mwAWA">y</span> <span class="o" id="mwAWE">=</span> <span class="n" id="mwAWI">ys</span><span class="p" id="mwAWM">[</span><span class="n" id="mwAWQ">i</span> <span class="o" id="mwAWU">-</span> <span class="mi" id="mwAWY">1</span><span class="p" id="mwAWc">]</span> <span class="linenos" data-line="44" id="mwAWg"></span> <span class="n" id="mwAWk">ys</span><span class="p" id="mwAWo">[</span><span class="n" id="mwAWs">i</span><span class="p" id="mwAWw">]</span> <span class="o" id="mwAW0">=</span> <span class="n" id="mwAW4">y</span> <span class="o" id="mwAW8">+</span> <span class="n" id="mwAXA">mu</span><span class="p" id="mwAXE">(</span><span class="n" id="mwAXI">y</span><span class="p" id="mwAXM">,</span> <span class="n" id="mwAXQ">t</span><span class="p" id="mwAXU">)</span> <span class="o" id="mwAXY">*</span> <span class="n" id="mwAXc">DT</span> <span class="o" id="mwAXg">+</span> <span class="n" id="mwAXk">sigma</span><span class="p" id="mwAXo">(</span><span class="n" id="mwAXs">y</span><span class="p" id="mwAXw">,</span> <span class="n" id="mwAX0">t</span><span class="p" id="mwAX4">)</span> <span class="o" id="mwAX8">*</span> <span class="n" id="mwAYA">dW</span><span class="p" id="mwAYE">(</span><span class="n" id="mwAYI">DT</span><span class="p" id="mwAYM">)</span> <span class="linenos" data-line="45" id="mwAYQ"></span> <span class="linenos" data-line="46" id="mwAYU"></span> <span class="k" id="mwAYY">return</span> <span class="n" id="mwAYc">TS</span><span class="p" id="mwAYg">,</span> <span class="n" id="mwAYk">ys</span> <span class="linenos" data-line="47" id="mwAYo"></span> <span class="linenos" data-line="48" id="mwAYs"></span> <span class="linenos" data-line="49" id="mwAYw"></span><span class="k" id="mwAY0">def</span> <span class="nf" id="mwAY4">plot_simulations</span><span class="p" id="mwAY8">(</span><span class="n" id="mwAZA">num_sims</span><span class="p" id="mwAZE">:</span> <span class="nb" id="mwAZI">int</span><span class="p" id="mwAZM">):</span> <span class="linenos" data-line="50" id="mwAZQ"></span><span class="w" id="mwAZU"> </span><span class="sd" id="mwAZY">""" Plot several simulations in one image."""</span> <span class="linenos" data-line="51" id="mwAZc"></span> <span class="k" id="mwAZg">for</span> <span class="n" id="mwAZk">_</span> <span class="ow" id="mwAZo">in</span> <span class="nb" id="mwAZs">range</span><span class="p" id="mwAZw">(</span><span class="n" id="mwAZ0">num_sims</span><span class="p" id="mwAZ4">):</span> <span class="linenos" data-line="52" id="mwAZ8"></span> <span class="n" id="mwAaA">plt</span><span class="o" id="mwAaE">.</span><span class="n" id="mwAaI">plot</span><span class="p" id="mwAaM">(</span><span class="o" id="mwAaQ">*</span><span class="n" id="mwAaU">run_simulation</span><span class="p" id="mwAaY">())</span> <span class="linenos" data-line="53" id="mwAac"></span> <span class="linenos" data-line="54" id="mwAag"></span> <span class="n" id="mwAak">plt</span><span class="o" id="mwAao">.</span><span class="n" id="mwAas">xlabel</span><span class="p" id="mwAaw">(</span><span class="s2" id="mwAa0">"time"</span><span class="p" id="mwAa4">)</span> <span class="linenos" data-line="55" id="mwAa8"></span> <span class="n" id="mwAbA">plt</span><span class="o" id="mwAbE">.</span><span class="n" id="mwAbI">ylabel</span><span class="p" id="mwAbM">(</span><span class="s2" id="mwAbQ">"y"</span><span class="p" id="mwAbU">)</span> <span class="linenos" data-line="56" id="mwAbY"></span> <span class="n" id="mwAbc">plt</span><span class="o" id="mwAbg">.</span><span class="n" id="mwAbk">show</span><span class="p" id="mwAbo">()</span> <span class="linenos" data-line="57" id="mwAbs"></span> <span class="linenos" data-line="58" id="mwAbw"></span> <span class="linenos" data-line="59" id="mwAb0"></span><span class="k" id="mwAb4">if</span> <span class="vm" id="mwAb8">__name__</span> <span class="o" id="mwAcA">==</span> <span class="s2" id="mwAcE">"__main__"</span><span class="p" id="mwAcI">:</span> <span class="linenos" data-line="60" id="mwAcM"></span> <span class="n" id="mwAcQ">NUM_SIMS</span> <span class="o" id="mwAcU">=</span> <span class="mi" id="mwAcY">5</span> <span class="linenos" data-line="61" id="mwAcc"></span> <span class="n" id="mwAcg">plot_simulations</span><span class="p" id="mwAck">(</span><span class="n" id="mwAco">NUM_SIMS</span><span class="p" id="mwAcs">)</span>

The following is simply the translation of the above code into the MATLAB (R2019b) programming language:

<span id="mwAdQ"></span><span class="linenos" data-line="1" id="mwAdU"></span><span class="c" id="mwAdY">%% Initialization and Utility</span> <span class="linenos" data-line="2" id="mwAdc"></span><span class="nb" id="mwAdg">close</span><span class="w" id="mwAdk"> </span><span class="nb" id="mwAdo">all</span><span class="p" id="mwAds">;</span> <span class="linenos" data-line="3" id="mwAdw"></span><span class="n" id="mwAd0">clear</span><span class="w" id="mwAd4"> </span><span class="s" id="mwAd8">all</span><span class="p" id="mwAeA">;</span> <span class="linenos" data-line="4" id="mwAeE"></span> <span class="linenos" data-line="5" id="mwAeI"></span><span class="n" id="mwAeM">numSims</span><span class="w" id="mwAeQ"> </span><span class="p" id="mwAeU">=</span><span class="w" id="mwAeY"> </span><span class="mi" id="mwAec">5</span><span class="p" id="mwAeg">;</span><span class="w" id="mwAek"> </span><span class="c" id="mwAeo">% display five runs</span> <span class="linenos" data-line="6" id="mwAes"></span><span class="n" id="mwAew">tBounds</span><span class="w" id="mwAe0"> </span><span class="p" id="mwAe4">=</span><span class="w" id="mwAe8"> </span><span class="p" id="mwAfA">[</span><span class="mi" id="mwAfE">3</span><span class="w" id="mwAfI"> </span><span class="mi" id="mwAfM">7</span><span class="p" id="mwAfQ">];</span><span class="w" id="mwAfU"> </span><span class="c" id="mwAfY">% The bounds of t</span> <span class="linenos" data-line="7" id="mwAfc"></span><span class="n" id="mwAfg">N</span><span class="w" id="mwAfk"> </span><span class="p" id="mwAfo">=</span><span class="w" id="mwAfs"> </span><span class="mi" id="mwAfw">1000</span><span class="p" id="mwAf0">;</span><span class="w" id="mwAf4"> </span><span class="c" id="mwAf8">% Compute at 1000 grid points</span> <span class="linenos" data-line="8" id="mwAgA"></span><span class="n" id="mwAgE">dt</span><span class="w" id="mwAgI"> </span><span class="p" id="mwAgM">=</span><span class="w" id="mwAgQ"> </span><span class="p" id="mwAgU">(</span><span class="n" id="mwAgY">tBounds</span><span class="p" id="mwAgc">(</span><span class="mi" id="mwAgg">2</span><span class="p" id="mwAgk">)</span><span class="w" id="mwAgo"> </span><span class="o" id="mwAgs">-</span><span class="w" id="mwAgw"> </span><span class="n" id="mwAg0">tBounds</span><span class="p" id="mwAg4">(</span><span class="mi" id="mwAg8">1</span><span class="p" id="mwAhA">))</span><span class="w" id="mwAhE"> </span><span class="o" id="mwAhI">/</span><span class="w" id="mwAhM"> </span><span class="n" id="mwAhQ">N</span><span class="w" id="mwAhU"> </span><span class="p" id="mwAhY">;</span> <span class="linenos" data-line="9" id="mwAhc"></span><span class="n" id="mwAhg">y_init</span><span class="w" id="mwAhk"> </span><span class="p" id="mwAho">=</span><span class="w" id="mwAhs"> </span><span class="mi" id="mwAhw">0</span><span class="p" id="mwAh0">;</span><span class="w" id="mwAh4"> </span><span class="c" id="mwAh8">% Initial y condition </span> <span class="linenos" data-line="10" id="mwAiA"></span> <span class="linenos" data-line="11" id="mwAiE"></span> <span class="linenos" data-line="12" id="mwAiI"></span><span class="c" id="mwAiM">% Initialize the probability distribution for our</span> <span class="linenos" data-line="13" id="mwAiQ"></span><span class="c" id="mwAiU">% random variable with mean 0 and stdev of sqrt(dt)</span> <span class="linenos" data-line="14" id="mwAiY"></span><span class="n" id="mwAic">pd</span><span class="w" id="mwAig"> </span><span class="p" id="mwAik">=</span><span class="w" id="mwAio"> </span><span class="n" id="mwAis">makedist</span><span class="p" id="mwAiw">(</span><span class="s" id="mwAi0">'Normal'</span><span class="p" id="mwAi4">,</span><span class="mi" id="mwAi8">0</span><span class="p" id="mwAjA">,</span><span class="nb" id="mwAjE">sqrt</span><span class="p" id="mwAjI">(</span><span class="n" id="mwAjM">dt</span><span class="p" id="mwAjQ">));</span> <span class="linenos" data-line="15" id="mwAjU"></span> <span class="linenos" data-line="16" id="mwAjY"></span><span class="n" id="mwAjc">c</span><span class="w" id="mwAjg"> </span><span class="p" id="mwAjk">=</span><span class="w" id="mwAjo"> </span><span class="p" id="mwAjs">[</span><span class="mf" id="mwAjw">0.7</span><span class="p" id="mwAj0">,</span><span class="w" id="mwAj4"> </span><span class="mf" id="mwAj8">1.5</span><span class="p" id="mwAkA">,</span><span class="w" id="mwAkE"> </span><span class="mf" id="mwAkI">0.06</span><span class="p" id="mwAkM">];</span><span class="w" id="mwAkQ"> </span><span class="c" id="mwAkU">% Theta, Mu, and Sigma, respectively</span> <span class="linenos" data-line="17" id="mwAkY"></span> <span class="linenos" data-line="18" id="mwAkc"></span><span class="n" id="mwAkg">ts</span><span class="w" id="mwAkk"> </span><span class="p" id="mwAko">=</span><span class="w" id="mwAks"> </span><span class="nb" id="mwAkw">linspace</span><span class="p" id="mwAk0">(</span><span class="n" id="mwAk4">tBounds</span><span class="p" id="mwAk8">(</span><span class="mi" id="mwAlA">1</span><span class="p" id="mwAlE">),</span><span class="w" id="mwAlI"> </span><span class="n" id="mwAlM">tBounds</span><span class="p" id="mwAlQ">(</span><span class="mi" id="mwAlU">2</span><span class="p" id="mwAlY">),</span><span class="w" id="mwAlc"> </span><span class="n" id="mwAlg">N</span><span class="p" id="mwAlk">);</span><span class="w" id="mwAlo"> </span><span class="c" id="mwAls">% From t0-->t1 with N points</span> <span class="linenos" data-line="19" id="mwAlw"></span><span class="n" id="mwAl0">ys</span><span class="w" id="mwAl4"> </span><span class="p" id="mwAl8">=</span><span class="w" id="mwAmA"> </span><span class="nb" id="mwAmE">zeros</span><span class="p" id="mwAmI">(</span><span class="mi" id="mwAmM">1</span><span class="p" id="mwAmQ">,</span><span class="n" id="mwAmU">N</span><span class="p" id="mwAmY">);</span><span class="w" id="mwAmc"> </span><span class="c" id="mwAmg">% 1xN Matrix of zeros</span> <span class="linenos" data-line="20" id="mwAmk"></span> <span class="linenos" data-line="21" id="mwAmo"></span><span class="n" id="mwAms">ys</span><span class="p" id="mwAmw">(</span><span class="mi" id="mwAm0">1</span><span class="p" id="mwAm4">)</span><span class="w" id="mwAm8"> </span><span class="p" id="mwAnA">=</span><span class="w" id="mwAnE"> </span><span class="n" id="mwAnI">y_init</span><span class="p" id="mwAnM">;</span> <span class="linenos" data-line="22" id="mwAnQ"></span><span class="c" id="mwAnU">%% Computing the Process</span> <span class="linenos" data-line="23" id="mwAnY"></span><span class="k" id="mwAnc">for</span><span class="w" id="mwAng"> </span><span class="nb" id="mwAnk">j</span><span class="w" id="mwAno"> </span><span class="p" id="mwAns">=</span><span class="w" id="mwAnw"> </span><span class="mi" id="mwAn0">1</span><span class="p" id="mwAn4">:</span><span class="n" id="mwAn8">numSims</span> <span class="linenos" data-line="24" id="mwAoA"></span><span class="w" id="mwAoE"> </span><span class="k" id="mwAoI">for</span><span class="w" id="mwAoM"> </span><span class="nb" id="mwAoQ">i</span><span class="w" id="mwAoU"> </span><span class="p" id="mwAoY">=</span><span class="w" id="mwAoc"> </span><span class="mi" id="mwAog">2</span><span class="p" id="mwAok">:</span><span class="nb" id="mwAoo">numel</span><span class="p" id="mwAos">(</span><span class="n" id="mwAow">ts</span><span class="p" id="mwAo0">)</span> <span class="linenos" data-line="25" id="mwAo4"></span><span class="w" id="mwAo8"> </span><span class="n" id="mwApA">t</span><span class="w" id="mwApE"> </span><span class="p" id="mwApI">=</span><span class="w" id="mwApM"> </span><span class="n" id="mwApQ">tBounds</span><span class="p" id="mwApU">(</span><span class="mi" id="mwApY">1</span><span class="p" id="mwApc">)</span><span class="w" id="mwApg"> </span><span class="o" id="mwApk">+</span><span class="w" id="mwApo"> </span><span class="p" id="mwAps">(</span><span class="nb" id="mwApw">i</span><span class="o" id="mwAp0">-</span><span class="mi" id="mwAp4">1</span><span class="p" id="mwAp8">)</span><span class="w" id="mwAqA"> </span><span class="o" id="mwAqE">.*</span><span class="w" id="mwAqI"> </span><span class="n" id="mwAqM">dt</span><span class="p" id="mwAqQ">;</span> <span class="linenos" data-line="26" id="mwAqU"></span><span class="w" id="mwAqY"> </span><span class="n" id="mwAqc">y</span><span class="w" id="mwAqg"> </span><span class="p" id="mwAqk">=</span><span class="w" id="mwAqo"> </span><span class="n" id="mwAqs">ys</span><span class="p" id="mwAqw">(</span><span class="nb" id="mwAq0">i</span><span class="o" id="mwAq4">-</span><span class="mi" id="mwAq8">1</span><span class="p" id="mwArA">);</span> <span class="linenos" data-line="27" id="mwArE"></span><span class="w" id="mwArI"> </span><span class="n" id="mwArM">mu</span><span class="w" id="mwArQ"> </span><span class="p" id="mwArU">=</span><span class="w" id="mwArY"> </span><span class="n" id="mwArc">c</span><span class="p" id="mwArg">(</span><span class="mi" id="mwArk">1</span><span class="p" id="mwAro">)</span><span class="w" id="mwArs"> </span><span class="o" id="mwArw">.*</span><span class="w" id="mwAr0"> </span><span class="p" id="mwAr4">(</span><span class="n" id="mwAr8">c</span><span class="p" id="mwAsA">(</span><span class="mi" id="mwAsE">2</span><span class="p" id="mwAsI">)</span><span class="w" id="mwAsM"> </span><span class="o" id="mwAsQ">-</span><span class="w" id="mwAsU"> </span><span class="n" id="mwAsY">y</span><span class="p" id="mwAsc">);</span> <span class="linenos" data-line="28" id="mwAsg"></span><span class="w" id="mwAsk"> </span><span class="n" id="mwAso">sigma</span><span class="w" id="mwAss"> </span><span class="p" id="mwAsw">=</span><span class="w" id="mwAs0"> </span><span class="n" id="mwAs4">c</span><span class="p" id="mwAs8">(</span><span class="mi" id="mwAtA">3</span><span class="p" id="mwAtE">);</span> <span class="linenos" data-line="29" id="mwAtI"></span><span class="w" id="mwAtM"> </span><span class="n" id="mwAtQ">dW</span><span class="w" id="mwAtU"> </span><span class="p" id="mwAtY">=</span><span class="w" id="mwAtc"> </span><span class="n" id="mwAtg">random</span><span class="p" id="mwAtk">(</span><span class="n" id="mwAto">pd</span><span class="p" id="mwAts">);</span> <span class="linenos" data-line="30" id="mwAtw"></span><span class="w" id="mwAt0"> </span> <span class="linenos" data-line="31" id="mwAt4"></span><span class="w" id="mwAt8"> </span><span class="n" id="mwAuA">ys</span><span class="p" id="mwAuE">(</span><span class="nb" id="mwAuI">i</span><span class="p" id="mwAuM">)</span><span class="w" id="mwAuQ"> </span><span class="p" id="mwAuU">=</span><span class="w" id="mwAuY"> </span><span class="n" id="mwAuc">y</span><span class="w" id="mwAug"> </span><span class="o" id="mwAuk">+</span><span class="w" id="mwAuo"> </span><span class="n" id="mwAus">mu</span><span class="w" id="mwAuw"> </span><span class="o" id="mwAu0">.*</span><span class="w" id="mwAu4"> </span><span class="n" id="mwAu8">dt</span><span class="w" id="mwAvA"> </span><span class="o" id="mwAvE">+</span><span class="w" id="mwAvI"> </span><span class="n" id="mwAvM">sigma</span><span class="w" id="mwAvQ"> </span><span class="o" id="mwAvU">.*</span><span class="w" id="mwAvY"> </span><span class="n" id="mwAvc">dW</span><span class="p" id="mwAvg">;</span> <span class="linenos" data-line="32" id="mwAvk"></span><span class="w" id="mwAvo"> </span><span class="k" id="mwAvs">end</span> <span class="linenos" data-line="33" id="mwAvw"></span><span class="w" id="mwAv0"> </span><span class="nb" id="mwAv4">figure</span><span class="p" id="mwAv8">;</span> <span class="linenos" data-line="34" id="mwAwA"></span><span class="w" id="mwAwE"> </span><span class="n" id="mwAwI">hold</span><span class="w" id="mwAwM"> </span><span class="s" id="mwAwQ">on</span><span class="p" id="mwAwU">;</span> <span class="linenos" data-line="35" id="mwAwY"></span><span class="w" id="mwAwc"> </span><span class="nb" id="mwAwg">plot</span><span class="p" id="mwAwk">(</span><span class="n" id="mwAwo">ts</span><span class="p" id="mwAws">,</span><span class="w" id="mwAww"> </span><span class="n" id="mwAw0">ys</span><span class="p" id="mwAw4">,</span><span class="w" id="mwAw8"> </span><span class="s" id="mwAxA">'o'</span><span class="p" id="mwAxE">)</span> <span class="linenos" data-line="36" id="mwAxI"></span><span class="k" id="mwAxM">end</span>

Example

Numerical simulation

Computer implementation

More in Computing and communication