differential_equations/tableau/
dorman_prince.rs

1//! Dormand-Prince Runge-Kutta methods with dense output interpolation.
2
3use crate::{tableau::ButcherTableau, traits::Real};
4
5impl<T: Real> ButcherTableau<T, 7> {
6    /// Dormand-Prince 5(4) Tableau with dense output interpolation.
7    ///
8    /// # Overview
9    /// This provides a 7-stage Runge-Kutta method with:
10    /// - Primary order: 5
11    /// - Embedded order: 4 (for error estimation)
12    /// - Number of stages: 7 primary + 0 additional stages for interpolation
13    /// - Built-in dense output of order 4
14    ///
15    /// # Efficiency
16    /// The DOPRI5 method is popular due to its efficient balance between accuracy and
17    /// computational cost. It is particularly good for non-stiff problems.
18    ///
19    /// # Interpolation
20    /// - The method provides a 4th-order interpolant using the existing 7 stages
21    /// - The interpolant has continuous first derivatives
22    /// - The interpolant uses the values of the stages to construct a polynomial
23    ///   that allows evaluation at any point within the integration step
24    ///
25    /// # Notes
26    /// - The method was developed by Dormand and Prince in 1980
27    /// - It is one of the most widely used Runge-Kutta methods and is implemented
28    ///   in many software packages for solving ODEs
29    /// - The DOPRI5 method is a member of the Dormand-Prince family of embedded
30    ///   Runge-Kutta methods
31    ///
32    /// # References
33    /// - Dormand, J. R. & Prince, P. J. (1980), "A family of embedded Runge-Kutta formulae",
34    ///   Journal of Computational and Applied Mathematics, 6(1), pp. 19-26
35    /// - Hairer, E., Nørsett, S. P. & Wanner, G. (1993), "Solving Ordinary Differential Equations I:
36    ///   Nonstiff Problems", Springer Series in Computational Mathematics, Vol. 8, Springer-Verlag
37    pub fn dopri5() -> Self {
38        let mut c = [0.0; 7];
39        let mut a = [[0.0; 7]; 7];
40        let mut b = [0.0; 7];
41        let mut bi4 = [[0.0; 7]; 7];
42        let mut er = [0.0; 7];
43
44        c[0] = 0.0;
45        c[1] = 0.2;
46        c[2] = 0.3;
47        c[3] = 0.8;
48        c[4] = 8.0 / 9.0;
49        c[5] = 1.0;
50        c[6] = 1.0;
51
52        a[1][0] = 0.2;
53
54        a[2][0] = 3.0 / 40.0;
55        a[2][1] = 9.0 / 40.0;
56
57        a[3][0] = 44.0 / 45.0;
58        a[3][1] = -56.0 / 15.0;
59        a[3][2] = 32.0 / 9.0;
60
61        a[4][0] = 19372.0 / 6561.0;
62        a[4][1] = -25360.0 / 2187.0;
63        a[4][2] = 64448.0 / 6561.0;
64        a[4][3] = -212.0 / 729.0;
65
66        a[5][0] = 9017.0 / 3168.0;
67        a[5][1] = -355.0 / 33.0;
68        a[5][2] = 46732.0 / 5247.0;
69        a[5][3] = 49.0 / 176.0;
70        a[5][4] = -5103.0 / 18656.0;
71
72        a[6][0] = 35.0 / 384.0;
73        a[6][1] = 0.0;
74        a[6][2] = 500.0 / 1113.0;
75        a[6][3] = 125.0 / 192.0;
76        a[6][4] = -2187.0 / 6784.0;
77        a[6][5] = 11.0 / 84.0;
78
79        b[0] = 35.0 / 384.0;
80        b[1] = 0.0;
81        b[2] = 500.0 / 1113.0;
82        b[3] = 125.0 / 192.0;
83        b[4] = -2187.0 / 6784.0;
84        b[5] = 11.0 / 84.0;
85        b[6] = 0.0;
86
87        er[0] = 71.0 / 57600.0;
88        er[1] = 0.0;
89        er[2] = -71.0 / 16695.0;
90        er[3] = 71.0 / 1920.0;
91        er[4] = -17253.0 / 339200.0;
92        er[5] = 22.0 / 525.0;
93        er[6] = -1.0 / 40.0;
94
95        bi4[0][0] = -12715105075.0 / 11282082432.0;
96        bi4[0][1] = 0.0;
97        bi4[0][2] = 87487479700.0 / 32700410799.0;
98        bi4[0][3] = -10690763975.0 / 1880347072.0;
99        bi4[0][4] = 701980252875.0 / 199316789632.0;
100        bi4[0][5] = -1453857185.0 / 822651844.0;
101        bi4[0][6] = 69997945.0 / 29380423.0;
102
103        let c = c.map(|x| T::from_f64(x).unwrap());
104        let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
105        let b = b.map(|x| T::from_f64(x).unwrap());
106        let er = er.map(|x| T::from_f64(x).unwrap());
107        let bi4 = bi4.map(|row| row.map(|x| T::from_f64(x).unwrap()));
108
109        ButcherTableau {
110            c,
111            a,
112            b,
113            bh: None,
114            bi: Some(bi4),
115            er: Some(er),
116        }
117    }
118}
119
120impl<T: Real> ButcherTableau<T, 12, 16> {
121    /// Dormand-Prince 8(5,3) Tableau with dense output interpolation.
122    ///
123    /// # Overview
124    /// This provides a 12-stage Runge-Kutta method with:
125    /// - Primary order: 8
126    /// - Embedded order: 5 (for error estimation)
127    /// - Number of stages: 12 primary + 4 additional stages for interpolation
128    /// - Built-in dense output of order 7
129    ///
130    /// # Efficiency
131    /// The DOP853 method is a high-order Runge-Kutta method that provides excellent
132    /// accuracy and efficiency for non-stiff problems. It's particularly suitable when
133    /// high precision is required.
134    ///
135    /// # Interpolation
136    /// - The method provides a 7th-order interpolant using the existing stages plus
137    ///   3 additional stages for dense output
138    /// - The interpolant has continuous derivatives
139    /// - The interpolation is performed through a sophisticated continuous extension
140    ///   that maintains high accuracy throughout the integration step
141    ///
142    /// # Notes
143    /// - The method was developed by Dormand, Prince and others as an extension
144    ///   of their earlier work
145    /// - It is one of the most accurate explicit Runge-Kutta implementations
146    ///   available for solving ODEs
147    /// - The DOP853 method is widely used in scientific computing where high
148    ///   precision is required
149    ///
150    /// # References
151    /// - Hairer, E., Nørsett, S. P. & Wanner, G. (1993), "Solving Ordinary Differential Equations I:
152    ///   Nonstiff Problems", Springer Series in Computational Mathematics, Vol. 8, Springer-Verlag
153    /// - Dormand, J. R. & Prince, P. J. (1980), "A family of embedded Runge-Kutta formulae",
154    ///   Journal of Computational and Applied Mathematics, 6(1), pp. 19-26
155    pub fn dop853() -> Self {
156        let mut c = [0.0; 16];
157        let mut a = [[0.0; 16]; 16];
158        let mut b = [0.0; 12];
159        let mut bh = [0.0; 12];
160        let mut bi7 = [[0.0; 16]; 16];
161        let mut er = [0.0; 12];
162
163        c[0] = 0.0;
164        c[1] = 5.260_015_195_876_773E-2;
165        c[2] = 7.890_022_793_815_16E-2;
166        c[3] = 1.183_503_419_072_274E-1;
167        c[4] = 2.816_496_580_927_726E-1;
168        c[5] = 3.333_333_333_333_333E-1;
169        c[6] = 0.25E+00;
170        c[7] = 3.076_923_076_923_077E-1;
171        c[8] = 6.512_820_512_820_513E-1;
172        c[9] = 0.6E+00;
173        c[10] = 8.571_428_571_428_571E-1;
174        c[11] = 1.0;
175
176        a[1][0] = 5.260_015_195_876_773E-2;
177
178        a[2][0] = 1.972_505_698_453_79E-2;
179        a[2][1] = 5.917_517_095_361_37E-2;
180
181        a[3][0] = 2.958_758_547_680_685E-2;
182        a[3][1] = 0.0;
183        a[3][2] = 8.876_275_643_042_054E-2;
184
185        a[4][0] = 2.413_651_341_592_667E-1;
186        a[4][1] = 0.0;
187        a[4][2] = -8.845_494_793_282_861E-1;
188        a[4][3] = 9.248_340_032_617_92E-1;
189
190        a[5][0] = 3.703_703_703_703_703_5E-2;
191        a[5][1] = 0.0;
192        a[5][2] = 0.0;
193        a[5][3] = 1.708_286_087_294_738_6E-1;
194        a[5][4] = 1.254_676_875_668_224_2E-1;
195
196        a[6][0] = 3.7109375E-2;
197        a[6][1] = 0.0;
198        a[6][2] = 0.0;
199        a[6][3] = 1.702_522_110_195_440_5E-1;
200        a[6][4] = 6.021_653_898_045_596E-2;
201        a[6][5] = -1.7578125E-2;
202
203        a[7][0] = 3.709_200_011_850_479E-2;
204        a[7][1] = 0.0;
205        a[7][2] = 0.0;
206        a[7][3] = 1.703_839_257_122_399_8E-1;
207        a[7][4] = 1.072_620_304_463_732_8E-1;
208        a[7][5] = -1.531_943_774_862_440_2E-2;
209        a[7][6] = 8.273_789_163_814_023E-3;
210
211        a[8][0] = 6.241_109_587_160_757E-1;
212        a[8][1] = 0.0;
213        a[8][2] = 0.0;
214        a[8][3] = -3.360_892_629_446_941_4;
215        a[8][4] = -8.682_193_468_417_26E-1;
216        a[8][5] = 2.759_209_969_944_671E1;
217        a[8][6] = 2.015_406_755_047_789_4E1;
218        a[8][7] = -4.348_988_418_106_996E1;
219
220        a[9][0] = 4.776_625_364_382_643_4E-1;
221        a[9][1] = 0.0;
222        a[9][2] = 0.0;
223        a[9][3] = -2.488_114_619_971_667_7;
224        a[9][4] = -5.902_908_268_368_43E-1;
225        a[9][5] = 2.123_005_144_818_119_3E1;
226        a[9][6] = 1.527_923_363_288_242_3E1;
227        a[9][7] = -3.328_821_096_898_486E1;
228        a[9][8] = -2.033_120_170_850_862_7E-2;
229
230        a[10][0] = -9.371_424_300_859_873E-1;
231        a[10][1] = 0.0;
232        a[10][2] = 0.0;
233        a[10][3] = 5.186_372_428_844_064;
234        a[10][4] = 1.091_437_348_996_729_5;
235        a[10][5] = -8.149_787_010_746_927;
236        a[10][6] = -1.852_006_565_999_696E1;
237        a[10][7] = 2.273_948_709_935_050_5E1;
238        a[10][8] = 2.493_605_552_679_652_3;
239        a[10][9] = -3.046_764_471_898_219_6;
240
241        a[11][0] = 2.273_310_147_516_538;
242        a[11][1] = 0.0;
243        a[11][2] = 0.0;
244        a[11][3] = -1.053_449_546_673_725E1;
245        a[11][4] = -2.000_872_058_224_862_5;
246        a[11][5] = -1.795_893_186_311_88E1;
247        a[11][6] = 2.794_888_452_941_996E1;
248        a[11][7] = -2.858_998_277_135_023_5;
249        a[11][8] = -8.872_856_933_530_63;
250        a[11][9] = 1.236_056_717_579_430_3E1;
251        a[11][10] = 6.433_927_460_157_636E-1;
252
253        b[0] = 5.429_373_411_656_876_5E-2;
254        b[1] = 0.0;
255        b[2] = 0.0;
256        b[3] = 0.0;
257        b[4] = 0.0;
258        b[5] = 4.450_312_892_752_409;
259        b[6] = 1.891_517_899_314_500_3;
260        b[7] = -5.801_203_960_010_585;
261        b[8] = 3.111_643_669_578_199E-1;
262        b[9] = -1.521_609_496_625_161E-1;
263        b[10] = 2.013_654_008_040_303_4E-1;
264        b[11] = 4.471_061_572_777_259E-2;
265
266        bh[0] = 2.440_944_881_889_764E-1;
267        bh[8] = 7.338_466_882_816_118E-1;
268        bh[11] = 2.205_882_352_941_176_6E-2;
269
270        er[0] = 1.312_004_499_419_488E-2;
271        er[5] = -1.225_156_446_376_204_4;
272        er[6] = -4.957_589_496_572_502E-1;
273        er[7] = 1.664_377_182_454_986_4;
274        er[8] = -3.503_288_487_499_736_6E-1;
275        er[9] = 3.341_791_187_130_175E-1;
276        er[10] = 8.192_320_648_511_571E-2;
277        er[11] = -2.235_530_786_388_629_4E-2;
278
279        // c[12] is the derivative of the new point at c[11] (1.0)
280
281        c[13] = 0.1E+00;
282
283        a[13][0] = 5.616_750_228_304_795_4E-2;
284        a[13][6] = 2.535_002_102_166_248_3E-1;
285        a[13][7] = -2.462_390_374_708_025E-1;
286        a[13][8] = -1.241_914_232_638_163_7E-1;
287        a[13][9] = 1.532_917_982_787_656_8E-1;
288        a[13][10] = 8.201_052_295_634_69E-3;
289        a[13][11] = 7.567_897_660_545_699E-3;
290        a[13][12] = -8.298E-3;
291
292        c[14] = 0.2E+00;
293
294        a[14][0] = 3.183_464_816_350_214E-2;
295        a[14][5] = 2.830_090_967_236_677_6E-2;
296        a[14][6] = 5.354_198_830_743_856_6E-2;
297        a[14][7] = -5.492_374_857_139_099E-2;
298        a[14][10] = -1.083_473_286_972_493_2E-4;
299        a[14][11] = 3.825_710_908_356_584E-4;
300        a[14][12] = -3.404_650_086_874_045_6E-4;
301        a[14][13] = 1.413_124_436_746_325E-1;
302
303        c[15] = 7.777_777_777_777_778E-1;
304
305        a[15][0] = -4.288_963_015_837_919_4E-1;
306        a[15][5] = -4.697_621_415_361_164;
307        a[15][6] = 7.683_421_196_062_599;
308        a[15][7] = 4.068_989_818_397_11;
309        a[15][8] = 3.567_271_874_552_811E-1;
310        a[15][12] = -1.399_024_165_159_014_5E-3;
311        a[15][13] = 2.947_514_789_152_772_4;
312        a[15][14] = -9.150_958_472_179_87;
313
314        bi7[4][0] = -8.428_938_276_109_013;
315        bi7[4][5] = 5.667_149_535_193_777E-1;
316        bi7[4][6] = -3.068_949_945_949_891_7;
317        bi7[4][7] = 2.384_667_656_512_07;
318        bi7[4][8] = 2.117_034_582_445_028;
319        bi7[4][9] = -8.713_915_837_779_73E-1;
320        bi7[4][10] = 2.240_437_430_260_788_3;
321        bi7[4][11] = 6.315_787_787_694_688E-1;
322        bi7[4][12] = -8.899_033_645_133_331E-2;
323        bi7[4][13] = 1.814_850_552_085_472_7E1;
324        bi7[4][14] = -9.194_632_392_478_356;
325        bi7[4][15] = -4.436_036_387_594_894;
326
327        bi7[5][0] = 1.042_750_864_257_913_4E1;
328        bi7[5][5] = 2.422_834_917_752_581_7E2;
329        bi7[5][6] = 1.652_004_517_172_702_8E2;
330        bi7[5][7] = -3.745_467_547_226_902E2;
331        bi7[5][8] = -2.211_366_685_312_530_6E1;
332        bi7[5][9] = 7.733_432_668_472_264;
333        bi7[5][10] = -3.067_408_473_108_939_8E1;
334        bi7[5][11] = -9.332_130_526_430_229;
335        bi7[5][12] = 1.569_723_812_177_084_5E1;
336        bi7[5][13] = -3.113_940_321_956_517_8E1;
337        bi7[5][14] = -9.352_924_358_844_48;
338        bi7[5][15] = 3.581_684_148_639_408E1;
339
340        bi7[6][0] = 1.998_505_324_200_243_3E1;
341        bi7[6][5] = -3.870_373_087_493_518E2;
342        bi7[6][6] = -1.891_781_381_951_675_8E2;
343        bi7[6][7] = 5.278_081_592_054_236E2;
344        bi7[6][8] = -1.157_390_253_995_963E1;
345        bi7[6][9] = 6.881_232_694_696_3;
346        bi7[6][10] = -1.000_605_096_691_083_8;
347        bi7[6][11] = 7.777_137_798_053_443E-1;
348        bi7[6][12] = -2.778_205_752_353_508;
349        bi7[6][13] = -6.019_669_523_126_412E1;
350        bi7[6][14] = 8.432_040_550_667_716E1;
351        bi7[6][15] = 1.199_229_113_618_279E1;
352
353        bi7[7][0] = -2.569_393_346_270_375E1;
354        bi7[7][5] = -1.541_897_486_902_364_3E2;
355        bi7[7][6] = -2.315_293_791_760_455E2;
356        bi7[7][7] = 3.576_391_179_106_141E2;
357        bi7[7][8] = 9.340_532_418_362_432E1;
358        bi7[7][9] = -3.745_832_313_645_163E1;
359        bi7[7][10] = 1.040_996_495_089_623E2;
360        bi7[7][11] = 2.984_029_342_666_05E1;
361        bi7[7][12] = -4.353_345_659_001_114E1;
362        bi7[7][13] = 9.632_455_395_918_828E1;
363        bi7[7][14] = -3.917_726_167_561_544E1;
364        bi7[7][15] = -1.497_268_362_579_856_4E2;
365
366        let c = c.map(|x| T::from_f64(x).unwrap());
367        let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
368        let b = b.map(|x| T::from_f64(x).unwrap());
369        let bh = bh.map(|x| T::from_f64(x).unwrap());
370        let bi7 = bi7.map(|row| row.map(|x| T::from_f64(x).unwrap()));
371        let er = er.map(|x| T::from_f64(x).unwrap());
372
373        ButcherTableau {
374            c,
375            a,
376            b,
377            bh: Some(bh),
378            bi: Some(bi7),
379            er: Some(er),
380        }
381    }
382}
differential_equations/tableau/dorman_prince.rs

differential_equations/tableau/
dorman_prince.rs
