![]() |
Adaptive time-stepping with high-order embedded Runge-Kutta pairs and rejection sampling provides efficient approaches for solving differential equations. While many such methods exist for solving deterministic systems, little progress has been made for stochastic variants. One challenge in developing adaptive methods for stochastic differential equations (SDEs) is the construction of embedded schemes with direct error estimates. We present a new class of embedded stochastic Runge-Kutta (SRK) methods with strong order 1.5 which have a natural embedding of strong order 1.0 methods. This allows for the derivation of an error estimate which requires no additional function evaluations. Next we derive a general method to reject the time steps without losing information about the future Brownian path termed Rejection Sampling with Memory (RSwM). This method utilizes a stack data structure to do rejection sampling, costing only a few floating point calculations. We show numerically that the methods generate statistically-correct and tolerance-controlled solutions. Lastly, we show that this form of adaptivity can be applied to systems of equations, and demonstrate that it solves a stiff biological model 12.28x faster than common fixed timestep algorithms. Our approach only requires the solution to a bridging problem and thus lends itself to natural generalizations beyond SDEs.
Citation: |
Figure 1. Outline of the adaptive SODE algorithm based on Rejection Sampling with Memory. This depicts the general schema for the Rejection Sampling with Memory algorithm (RSwM). The green arrow depicts the path taken when the step is rejected, whereas the red arrow depicts the path that is taken when the step is accepted. The blue text denotes steps which are specific to the stochastic systems in order to ensure correct sampling properties which are developed in section 4.
Figure 2.
Adaptive algorithm Kolmogorov-Smirnov Tests. Equation 31 was solved from
Figure 3.
Adaptive algorithm correctness checks. QQplots of the distribution of the Brownian path at the end time
Figure 4.
Adaptive Algorithm Error Comparison on Equation 31, Equation 33, and Equation 35. Comparison of the rejection sampling with memory algorithms on examples 1-3. Along the
Figure 5.
Adaptive Algorithm Timing Comparison. Comparison of the rejection sampling with memory algorithms. Along the
Figure 7.
Stochastic cell differentiation model solutions. (A) Timeseries of the concentration of
Table 1.
ESRK1. Table (a) shows the legend for how the numbers in in Table (b) correspond to the coefficient arrays/matrices
![]() |
Algorithm 1 RSwM1 | |
1: | Set the values |
2: | Set |
3: | Take an initial |
4: | While |
5: | Attempt a step with |
6: | Calculate E according to (9) |
7: | Update |
8: | If |
9: | Take |
10: | Calculate |
11: | Push |
12: | Update |
13: | Update |
14: | else ▷% Accept the Step |
15: | Update |
16: | if ( |
17: | Update |
18: | Take |
19: | else |
20: | Pop the top of |
21: | Update |
22: | end if |
23: | end if |
24: | end while |
Algorithm 2 RSwM3 | |
1: | Set the values |
2: | Set |
3: | Take an initial |
4: | while |
5: | Attempt a step with |
6: | Calculate E according to (9) |
7: | Update |
8: | if |
9: | Set |
10: | while |
11: | Pop the top of |
12: | if |
13: | Update |
14: | Update |
15: | else |
16: | Push |
17: | end if |
18: | end while |
19: | Set |
20: | Set |
21: | Take |
22: | Take |
23: | Pop |
24: | Pop |
25: | Update |
26: | else ▷% Accept the Step |
27: | Update |
28: | Empty |
29: | Update |
30: | Set |
31: | while |
32: | Pop the top of |
33: | if ( |
34: | Update |
35: | Push a copy of |
36: | else ▷% Final part of step from stack |
37: | Set |
38: | Let |
39: | Let |
40: | Push |
41: | Update |
42: | end if |
43: | end while |
▷% Update for last portion to step. Note zero if final part is from stack | |
44: | if ( |
45: | Let |
46: | Update |
47: | Push |
48: | end if |
49: | end if |
50: | end while |
Table 2.
Fixed timestep method fails and runtimes. The fixed timestep algorithms and ESRK+RSwM3 algorithms were used to solve the stochastic cell model on
Euler-Maruyama | Runge-Kutta Milstein | Rößler SRI | ||||
|
Fails (/10,000) | Time (s) | Fails (/10,000) | Time (s) | Fails (/10,000) | Time (s) |
|
137 | 133.35 | 131 | 211.92 | 78 | 609.27 |
|
39 | 269.09 | 26 | 428.28 | 17 | 1244.06 |
|
3 | 580.14 | 6 | 861.01 | 0 | 2491.37 |
|
1 | 1138.41 | 1 | 1727.91 | 0 | 4932.70 |
|
0 | 2286.35 | 0 | 3439.90 | 0 | 9827.16 |
|
0 | 4562.20 | 0 | 6891.35 | 0 | 19564.16 |
Table 3.
Example 1 | Example 2 | Example 3 | Cell Model | |||||
qmax | Time (s) | Error | Time (s) | Error | Time (s) | Error | Time (s) | |
|
37.00 | 2.57e-8 | 60.87 | 2.27e-7 | 67.71 | 3.42e-9 | 229.83 | |
|
34.73 | 2.82e-8 | 32.40 | 3.10e-7 | 66.68 | 3.43e-9 | 196.36 | |
|
49.14 | 3.14e-8 | 132.33 | 8.85e-7 | 65.94 | 3.44e-9 | 186.81 | |
|
39.33 | 3.59e-8 | 33.90 | 1.73e-6 | 66.33 | 3.44e-9 | 205.57 | |
|
38.22 | 3.82e-8 | 159.94 | 2.58e-6 | 68.16 | 3.44e-9 | 249.77 | |
|
82.76 | 4.41e-8 | 34.41 | 3.58e-6 | 568.22 | 3.44e-9 | 337.99 | |
|
68.16 | 9.63e-8 | 33.98 | 6.06e-6 | 87.50 | 3.22e-9 | 418.78 | |
|
48.23 | 1.01e-7 | 33.97 | 9.74e-6 | 69.78 | 3.44e-9 | 571.59 |
Algorithm 3 Initial |
|
1: | Let |
2: | Calculate |
3: | Let |
4: | if |
5: | Let |
6: | else |
7: | Let |
8: | end if |
9: | Calculate an Euler step: |
10: | Calculate new estimates: |
11: | Determine |
12: | Let |
13: | if |
14: | Let |
15: | else |
16: | Let |
17: | end if |
18: | Let |
Table 4.
SRA1. Table (a) shows the legend for how the numbers in in Table (b) correspond to the coefficient arrays/matrices
![]() |
Algorithm 4 RSwM2 | |
1: | Set the values |
2: | Set |
3: | Take an initial |
4: | while |
5: | Attempt a step with |
6: | Calculate E according to (9) |
7: | Update |
8: | if |
9: | Take |
10: | Calculate |
11: | Push |
12: | Update |
13: | Update |
14: | else ▷% Accept the Step |
15: | Update |
16: | Update |
17: | Set |
18: | while |
19: | Pop the top of |
20: | if ( |
21: | Update |
22: | else ▷% Final part of step from stack |
23: | Set |
24: | Let |
25: | Let |
26: | Push |
27: | Update |
28: | end if |
29: | end while |
▷% Update for last portion to step. Note zero if final part is from stack | |
30: | if ( |
31: | Let |
32: | Update |
33: | end if |
34: | end if |
35: | end while |
Table 5. Table of Parameter Values for the Stochastic Cell Model.
Parameter | Value | Parameter | Value | Parameter | Value | Parameter | Value |
|
3 | |
0.1 | |
1 | |
0.35 |
|
0.2 | |
0.3 | |
1 | |
0.0002 |
|
0.15 | |
0.4 | |
1 | |
0.001 |
|
0.35 | |
0.4 | |
1 | |
0.09 |
|
0.9 | |
2 | |
20 | |
0.1 |
|
0.6 | |
3.5 | |
100 | |
0.1 |
|
0.5 | |
0.9 | |
0 | |
0.9 |
|
1.8 | |
1 | |
1000 | |
1.66 |
|
0.0005 | |
0.003 | |
0.5 | |
1.1 |
|
3 | |
2 | |
0.5 | |
5 |
|
2 | |
2 | |
0.5 | |
5 |
|
2 | |
2 | |
0.5 | |
15 |
|
2 | |
2 | |
0.5 | |
5 |
|
2 | |
2 | |
0.5 | |
2 |
|
2 | |
6 | |
0.5 | |
5 |
|
1.2 | |
0.02 | |
0.01 | |
0.05 |
|
0.06 | |
1.5 | |
16 | |
16 |
|
0.5 | |
0.5 | |
0.5 | |
0.5 |
|
0.5 | |
1.0 | |
0.035 | |
0.035 |
|
0.9 | |
0.05 | |
0.05 |
[1] |
R. L. Burden and J. D. Faires,
Numerical Analysis 9th ed. , Brooks/Cole, Cengage Learning, Boston, MA, 2011.
![]() |
[2] |
K. Burrage and P. Burrage, General order conditions for stochastic runge-kutta methods for both commuting and non-commuting stochastic ordinary differential equation systems, Appl. Numer. Math., 28 (1998), 161-177.
doi: 10.1016/S0168-9274(98)00042-7.![]() ![]() ![]() |
[3] |
K. Burrage and P. M. Burrage, High strong order explicit runge-kutta methods for stochastic ordinary differential equations, Applied Numerical Mathematics, 22 (1996), 81-101.
doi: 10.1016/S0168-9274(96)00027-X.![]() ![]() ![]() |
[4] |
P. M. Burrage,
Runge-Kutta Methods for Stochastic Differential Equations Thesis, The University of Queensland Brisbane, 1999.
![]() |
[5] |
P. M. Burrage and K. Burrage, A variable stepsize implementation for stochastic differential equations, SIAM Journal on Scientific Computing, 24 (2002), 848-864.
doi: 10.1137/S1064827500376922.![]() ![]() ![]() |
[6] |
J. R. Cash and A. H. Karp, A variable order runge-kutta method for initial value problems with rapidly varying right-hand sides, ACM Transactions on Mathematical Software, 16 (1990), 201-222.
doi: 10.1145/79505.79507.![]() ![]() ![]() |
[7] |
F. Ceschino, Modification de la longueur du pas dans l'integration numerique parles methodes a pas lies, Chiffres, 4 (1961), 101-106.
![]() ![]() |
[8] |
J. R. Dormand and P. J. Prince, A family of embedded runge kutta formulae, Journal of Computational and Applied Mathematics, 6 (1980), 19-26.
doi: 10.1016/0771-050X(80)90013-3.![]() ![]() ![]() |
[9] |
J. G. Gaines and T. J. Lyons, Variable step size control in the numerical solution of stochastic differential equations, SIAM Journal on Applied Mathematics, 57 (1997), 1455-1484.
doi: 10.1137/S0036139995286515.![]() ![]() ![]() |
[10] |
A. Ghasemi and S. Zahediasl, Normality tests for statistical analysis: A guide for non-statisticians, Int J Endocrinol Metab, 10 (2012), 486-489.
doi: 10.5812/ijem.3505.![]() ![]() |
[11] |
E. Hairer, S. P. N orsett and G. Wanner,
Solving Ordinary Differential Equations I 2nd rev. edn, Springer series in computational mathematics. Springer-Verlag, Berlin, New York, 1993.
![]() ![]() |
[12] |
T. Hong, K. Watanabe, C. H. Ta, A. Villarreal-Ponce, Q. Nie and X. Dai, An ovol2-zeb1 mutual inhibitory circuit governs bidirectional and multi-step transition between epithelial and mesenchymal states PLoS Comput Biol 11 (2015), e1004569.
doi: 10.1371/journal. pcbi. 1004569.![]() ![]() |
[13] |
S. M. Iacus,
Simulation and Inference for Stochastic Differential Equations: With R Examples (Springer Series in Statistics), Springer Publishing Company, Incorporated, 2008.
![]() |
[14] |
J. Kaneko, Explicit order 1.5 runge kutta scheme for solutions of ito stochastic differential equations, S/urikaisekikenky/usho K/oky/uroku, 932 (1995), 46-60.
![]() ![]() |
[15] |
P. E. Kloeden and E. Platen,
Numerical Solution of Stochastic Differential Equations Springer Berlin Heidelberg, 1992.
doi: 10.1007/978-3-662-12616-5.![]() ![]() ![]() |
[16] |
H. Lamba, An adaptive timestepping algorithm for stochastic differential equations, Journal of Computational and Applied Mathematics, 161 (2003), 417-430.
doi: 10.1016/j.cam.2003.05.001.![]() ![]() ![]() |
[17] |
T. Liggett,
Continuous Time Markov Processes: An Introduction American Mathematical Society, 2010.
doi: 10.1090/gsm/113.![]() ![]() ![]() |
[18] |
X. Mao and C. Yuan,
Stochastic Differential Equations with Markovian Switching Imperial College Press, London, 2006.
doi: 10.1142/p473.![]() ![]() ![]() |
[19] |
N. Hofmann, T. Muller-Gronbach and K. Ritter, Optimal approximation of stochastic differential equations by adaptive step-size control, Mathematics of Computation, 69 (2000), 1017-1034.
doi: 10.1090/S0025-5718-99-01177-1.![]() ![]() ![]() |
[20] |
B. Oksendal,
Stochastic Differential Equations: An Introduction with Applications Springer, 2003.
doi: 10.1007/978-3-642-14394-6.![]() ![]() ![]() |
[21] |
W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery,
Numerical Recipes 3rd Edition: The Art of Scientific Computing, Cambridge University Press, 2007.
![]() ![]() |
[22] |
A. Rößler, Runge kutta methods for the strong approximation of solutions of stochastic differential equations, SIAM Journal on Numerical Analysis, 48 (2010), 922-952.
doi: 10.1137/09076636X.![]() ![]() ![]() |
[23] |
T. Ryden and M. Wiktorsson, On the simulation of iterated ito integrals, Stochastic Processes and their Applications, 91 (2001), 151-168.
doi: 10.1016/S0304-4149(00)00053-3.![]() ![]() ![]() |
[24] |
L. F. Shampine, Some practical runge-kutta formulas, Mathematics of Computation, 46 (1986), 135-150.
doi: 10.1090/S0025-5718-1986-0815836-3.![]() ![]() ![]() |
[25] |
M. Wiktorsson, Joint characteristic function and simultaneous simulation of iterated ito integrals for multiple independent brownian motions, The Annals of Applied Probability, 11 (2001), 470-487.
doi: 10.1214/aoap/1015345301.![]() ![]() ![]() |