Runge–Kutta method (SDE) explained

In mathematics of stochastic systems, the Runge–Kutta method is a technique for the approximate numerical solution of a stochastic differential equation. It is a generalisation of the Runge–Kutta method for ordinary differential equations to stochastic differential equations (SDEs). Importantly, the method does not involve knowing derivatives of the coefficient functions in the SDEs.

Most basic scheme

X satisfying the following Itō stochastic differential equation$$dX\_t\; =\; a(X\_t)\; \backslash ,\; dt\; +\; b(X\_t)\; \backslash ,\; dW\_t,$$with initial condition X _{0=x} _{0} , where W _{t} stands for the Wiener process, and suppose that we wish to solve this SDE on some interval of time [ 0 ,T] . Then the basic Runge–Kutta approximation to the true solution X is the Markov chain Y defined as follows:^[1]

• partition the interval

[ 0 ,T] into N subintervals of width \delta = T/N > 0 :

• set

Y _{0} : = x _{0} ;

• recursively compute

Y _{n} for 1 \leq n\leq N by $$Y\_\; :=\; Y\_\; +\; a(Y\_)\; \backslash delta\; +\; b(Y\_)\; \backslash Delta\; W\_\; +\; \backslash frac\; \backslash left(b(\backslash hat\_)\; -\; b(Y\_)\; \backslash right)\; \backslash left((\backslash Delta\; W\_)^\; -\; \backslash delta\; \backslash right)\; \backslash delta^,$$ where \Delta W _{n} = W _{\tau _{n}} - W _{\tau _{n}} and \hat{\Upsilon} _{n} = Y _{n} + a(Y _{n)} \delta + b(Y _{n} ) \delta ^{1/2} . The random variables \Delta W _{n} are independent and identically distributed normal random variables with expected value zero and variance \delta .

This scheme has strong order 1, meaning that the approximation error of the actual solution at a fixed time scales with the time step

\delta . It has also weak order 1, meaning that the error on the statistics of the solution scales with the time step \delta . See the references for complete and exact statements.

The functions

a and b can be time-varying without any complication. The method can be generalized to the case of several coupled equations; the principle is the same but the equations become longer.

Variation of the Improved Euler is flexible

A newer Runge—Kutta scheme also of strong order 1 straightforwardly reduces to the improved Euler scheme for deterministic ODEs.^[2] Consider the vector stochastic process

\vec X(t)\in R ^{n} that satisfies the general Ito SDE$$d\backslash vec\; X=\backslash vec\; a(t,\backslash vec\; X)\backslash ,dt+\backslash vec\; b(t,\backslash vec\; X)\backslash ,dW,$$where drift \vec a and volatility \vec b are sufficiently smooth functions of their arguments.Given time step h , and given the value estimate \vec X(t _{k+1} ) by \vec X _{k+1} for time t _{k+1} = t _{k+h} via$$\backslash begin\backslash vec\; K\_1=h\backslash vec\; a(t\_k,\backslash vec\; X\_k)+(\backslash Delta\; W\_k-S\_k\backslash sqrt\; h)\backslash vec\; b(t\_k,\backslash vec\; X\_k),\backslash \backslash \backslash vec\; K\_2=h\backslash vec\; a(t\_,\backslash vec\; X\_k+\backslash vec\; K\_1)+(\backslash Delta\; W\_k+S\_k\backslash sqrt\; h)\backslash vec\; b(t\_,\backslash vec\; X\_k+\backslash vec\; K\_1),\backslash \backslash \backslash vec\; X\_=\backslash vec\; X\_k+\backslash frac12(\backslash vec\; K\_1+\backslash vec\; K\_2),\backslash end$$

• where

\Delta W _{k=\sqrt} hZ _{k} for normal random Z _{k\sim} N( 0 , 1 ) ;

• and where

S _{k=\pm1} , each alternative chosen with probability 1 / 2 .

The above describes only one time step.Repeat this time step

(t _{m-t} _{0)/h} times in order to integrate the SDE from time t = t _{0} to t = t _{m} .

The scheme integrates Stratonovich SDEs to

O(h) provided one sets S _{k=0} throughout (instead of choosing \pm 1 ).

Higher order Runge-Kutta schemes

Higher-order schemes also exist, but become increasingly complex.Rößler developed many schemes for Ito SDEs,^{[3] [4]} whereas Komori developed schemes for Stratonovich SDEs.^{[5] [6] [7]} Rackauckas extended these schemes to allow for adaptive-time stepping via Rejection Sampling with Memory (RSwM), resulting in orders of magnitude efficiency increases in practical biological models,^[8] along with coefficient optimization for improved stability.^[9]

Notes and References

1. P. E. Kloeden and E. Platen. Numerical solution of stochastic differential equations, volume 23 of Applications of Mathematics. Springer--Verlag, 1992.
2. Roberts . A. J. . Oct 2012 . Modify the Improved Euler scheme to integrate stochastic differential equations . math.NA . 1210.0933 .
3. 10.1137/060673308. Second Order Runge–Kutta Methods for Itô Stochastic Differential Equations. SIAM Journal on Numerical Analysis. 47. 3. 1713–1738. 2009. Rößler . A. .
4. 10.1137/09076636X . Runge–Kutta Methods for the Strong Approximation of Solutions of Stochastic Differential Equations. SIAM Journal on Numerical Analysis. 48. 3. 922–952. 2010. Rößler . A. .
5. 10.1016/j.apnum.2006.02.002. Multi-colored rooted tree analysis of the weak order conditions of a stochastic Runge–Kutta family. Applied Numerical Mathematics. 57. 2. 147–165. 2007. Komori . Y. . 49220399 .
6. 10.1016/j.cam.2006.03.010. Weak order stochastic Runge–Kutta methods for commutative stochastic differential equations. Journal of Computational and Applied Mathematics. 203. 57–79. 2007. Komori . Y. . free.
7. 10.1016/j.cam.2006.06.006. Weak second-order stochastic Runge–Kutta methods for non-commutative stochastic differential equations. Journal of Computational and Applied Mathematics. 206. 158–173. 2007. Komori . Y. . free.
8. 10.3934/dcdsb.2017133. Discrete and Continuous Dynamical Systems - Series B. 22. 7. 2731–2761. 2017. Rackauckas. Christopher . Nie . Qing . Adaptive methods for stochastic differential equations via natural embeddings and rejection sampling with memory. 29527134. 5844583. free.
9. Rackauckas. Christopher . Nie . Qing . 2018 . Stability-optimized high order methods and stiffness detection for pathwise stiff stochastic differential equations . math.NA . 1804.04344.