differential_equations/tableau/runge_kutta.rs
1//! Classic or typical Runge-Kutta methods without unique properties
2
3use crate::{tableau::ButcherTableau, traits::Real};
4
5// -- Explicit Runge-Kutta methods (ERK) --
6
7impl<T: Real> ButcherTableau<T, 4> {
8 /// Classic Runge-Kutta 4th order method (RK4).
9 ///
10 /// # Overview
11 /// This provides a 4-stage, explicit Runge-Kutta method with:
12 /// - Primary order: 4
13 /// - Embedded order: None (this implementation does not include an embedded error estimate)
14 /// - Number of stages: 4
15 ///
16 /// # Interpolation
17 /// - This standard RK4 implementation does not provide coefficients for dense output (interpolation).
18 ///
19 /// # Notes
20 /// - RK4 is a widely used, general-purpose method known for its balance of accuracy and
21 /// computational efficiency for many problems.
22 /// - It is a fixed-step method as presented here, lacking an embedded formula for adaptive
23 /// step-size control. For adaptive step sizes, a method with an error estimate (e.g., RKF45)
24 /// would be required.
25 ///
26 /// # Butcher Tableau
27 /// ```text
28 /// 0 |
29 /// 1/2 | 1/2
30 /// 1/2 | 0 1/2
31 /// 1 | 0 0 1
32 /// ----|--------------------
33 /// | 1/6 1/3 1/3 1/6
34 /// ```
35 ///
36 /// # References
37 /// - Kutta, W. (1901). "Beitrag zur näherungsweisen Integration totaler Differentialgleichungen". *Zeitschrift für Mathematik und Physik*, 46, 435-453.
38 /// - Runge, C. (1895). "Über die numerische Auflösung von Differentialgleichungen". *Mathematische Annalen*, 46, 167-178.
39 pub fn rk4() -> Self {
40 let mut c = [0.0; 4];
41 let mut a = [[0.0; 4]; 4];
42 let mut b = [0.0; 4];
43
44 c[0] = 0.0;
45 c[1] = 0.5;
46 c[2] = 0.5;
47 c[3] = 1.0;
48
49 a[1][0] = 0.5;
50 a[2][1] = 0.5;
51 a[3][2] = 1.0;
52
53 b[0] = 1.0 / 6.0;
54 b[1] = 1.0 / 3.0;
55 b[2] = 1.0 / 3.0;
56 b[3] = 1.0 / 6.0;
57
58 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
59 let b = b.map(|x| T::from_f64(x).unwrap());
60 let c = c.map(|x| T::from_f64(x).unwrap());
61
62 Self {
63 c,
64 a,
65 b,
66 bh: None,
67 bi: None,
68 er: None,
69 }
70 }
71
72 /// Three-Eighths Rule 4th order method.
73 ///
74 /// # Overview
75 /// This provides a 4-stage, explicit Runge-Kutta method with:
76 /// - Primary order: 4
77 /// - Embedded order: None (no embedded error estimate)
78 /// - Number of stages: 4
79 ///
80 /// # Notes
81 /// - The primary advantage of this method is that almost all of the error coefficients
82 /// are smaller than in the classic RK4 method, but it requires slightly more FLOPs
83 /// (floating-point operations) per time step.
84 ///
85 /// # Butcher Tableau
86 /// ```text
87 /// 0 |
88 /// 1/3 | 1/3
89 /// 2/3 | -1/3 1
90 /// 1 | 1 -1 1
91 /// ----|--------------------
92 /// | 1/8 3/8 3/8 1/8
93 /// ```
94 ///
95 /// # References
96 /// - Butcher, J.C. (2008). "Numerical Methods for Ordinary Differential Equations".
97 pub fn three_eighths() -> Self {
98 let mut c = [0.0; 4];
99 let mut a = [[0.0; 4]; 4];
100 let mut b = [0.0; 4];
101
102 c[0] = 0.0;
103 c[1] = 1.0 / 3.0;
104 c[2] = 2.0 / 3.0;
105 c[3] = 1.0;
106
107 a[1][0] = 1.0 / 3.0;
108 a[2][0] = -1.0 / 3.0;
109 a[2][1] = 1.0;
110 a[3][0] = 1.0;
111 a[3][1] = -1.0;
112 a[3][2] = 1.0;
113
114 b[0] = 1.0 / 8.0;
115 b[1] = 3.0 / 8.0;
116 b[2] = 3.0 / 8.0;
117 b[3] = 1.0 / 8.0;
118
119 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
120 let b = b.map(|x| T::from_f64(x).unwrap());
121 let c = c.map(|x| T::from_f64(x).unwrap());
122
123 Self {
124 c,
125 a,
126 b,
127 bh: None,
128 bi: None,
129 er: None,
130 }
131 }
132}
133
134impl<T: Real> ButcherTableau<T, 2> {
135 /// Midpoint method (2nd order Runge-Kutta).
136 ///
137 /// # Overview
138 /// This provides a 2-stage, explicit Runge-Kutta method with:
139 /// - Primary order: 2
140 /// - Embedded order: None
141 /// - Number of stages: 2
142 ///
143 /// # Butcher Tableau
144 /// ```text
145 /// 0 |
146 /// 1/2 | 1/2
147 /// ----|--------
148 /// | 0 1
149 /// ```
150 ///
151 /// # References
152 /// - Butcher, J.C. (2008). "Numerical Methods for Ordinary Differential Equations".
153 pub fn midpoint() -> Self {
154 let mut c = [0.0; 2];
155 let mut a = [[0.0; 2]; 2];
156 let mut b = [0.0; 2];
157
158 c[0] = 0.0;
159 c[1] = 0.5;
160
161 a[1][0] = 0.5;
162
163 b[0] = 0.0;
164 b[1] = 1.0;
165
166 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
167 let b = b.map(|x| T::from_f64(x).unwrap());
168 let c = c.map(|x| T::from_f64(x).unwrap());
169
170 Self {
171 c,
172 a,
173 b,
174 bh: None,
175 bi: None,
176 er: None,
177 }
178 }
179
180 /// Heun's method (2nd order Runge-Kutta).
181 ///
182 /// # Overview
183 /// This provides a 2-stage, explicit Runge-Kutta method with:
184 /// - Primary order: 2
185 /// - Embedded order: None
186 /// - Number of stages: 2
187 ///
188 /// # Butcher Tableau
189 /// ```text
190 /// 0 |
191 /// 1 | 1
192 /// ----|--------
193 /// | 1/2 1/2
194 /// ```
195 ///
196 /// # References
197 /// - Heun, K. (1900). "Neue Methoden zur approximativen Integration der Differentialgleichungen einer unabhängigen Veränderlichen".
198 pub fn heun() -> Self {
199 let mut c = [0.0; 2];
200 let mut a = [[0.0; 2]; 2];
201 let mut b = [0.0; 2];
202
203 c[0] = 0.0;
204 c[1] = 1.0;
205
206 a[1][0] = 1.0;
207
208 b[0] = 0.5;
209 b[1] = 0.5;
210
211 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
212 let b = b.map(|x| T::from_f64(x).unwrap());
213 let c = c.map(|x| T::from_f64(x).unwrap());
214
215 Self {
216 c,
217 a,
218 b,
219 bh: None,
220 bi: None,
221 er: None,
222 }
223 }
224
225 /// Ralston's method (2nd order Runge-Kutta).
226 ///
227 /// # Overview
228 /// This provides a 2-stage, explicit Runge-Kutta method with:
229 /// - Primary order: 2
230 /// - Embedded order: None
231 /// - Number of stages: 2
232 ///
233 /// # Butcher Tableau
234 /// ```text
235 /// 0 |
236 /// 2/3 | 2/3
237 /// ----|--------
238 /// | 1/4 3/4
239 /// ```
240 ///
241 /// # References
242 /// - Ralston, A. (1962). "Runge-Kutta Methods with Minimum Error Bounds".
243 pub fn ralston() -> Self {
244 let mut c = [0.0; 2];
245 let mut a = [[0.0; 2]; 2];
246 let mut b = [0.0; 2];
247
248 c[0] = 0.0;
249 c[1] = 2.0 / 3.0;
250
251 a[1][0] = 2.0 / 3.0;
252
253 b[0] = 1.0 / 4.0;
254 b[1] = 3.0 / 4.0;
255
256 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
257 let b = b.map(|x| T::from_f64(x).unwrap());
258 let c = c.map(|x| T::from_f64(x).unwrap());
259
260 Self {
261 c,
262 a,
263 b,
264 bh: None,
265 bi: None,
266 er: None,
267 }
268 }
269}
270
271impl<T: Real> ButcherTableau<T, 1> {
272 /// Euler's method (1st order Runge-Kutta).
273 ///
274 /// # Overview
275 /// This provides a 1-stage, explicit Runge-Kutta method with:
276 /// - Primary order: 1
277 /// - Embedded order: None
278 /// - Number of stages: 1
279 ///
280 /// # Butcher Tableau
281 /// ```text
282 /// 0 |
283 /// --|--
284 /// | 1
285 /// ```
286 ///
287 /// # References
288 /// - Euler, L. (1768). "Institutionum calculi integralis".
289 pub fn euler() -> Self {
290 let c = [T::zero()];
291 let a = [[T::zero()]];
292 let b = [T::one()];
293
294 Self {
295 c,
296 a,
297 b,
298 bh: None,
299 bi: None,
300 er: None,
301 }
302 }
303}
304
305// Implementations for methods with embedded error estimators (adaptive methods)
306impl<T: Real> ButcherTableau<T, 6> {
307 /// Runge-Kutta-Fehlberg 4(5) method (RKF45).
308 ///
309 /// # Overview
310 /// This provides a 6-stage, explicit Runge-Kutta method with:
311 /// - Primary order: 5
312 /// - Embedded order: 4 (for error estimation)
313 /// - Number of stages: 6
314 ///
315 /// # Notes
316 /// - RKF45 is a widely used adaptive step size method that provides a good balance
317 /// between accuracy and computational efficiency.
318 /// - It uses the difference between 4th and 5th order approximations to estimate error.
319 ///
320 /// # Butcher Tableau
321 /// ```text
322 /// 0 |
323 /// 1/4 | 1/4
324 /// 3/8 | 3/32 9/32
325 /// 12/13 | 1932/2197 -7200/2197 7296/2197
326 /// 1 | 439/216 -8 3680/513 -845/4104
327 /// 1/2 | -8/27 2 -3544/2565 1859/4104 -11/40
328 /// -------|---------------------------------------------------------------
329 /// | 16/135 0 6656/12825 28561/56430 -9/50 2/55 (5th)
330 /// | 25/216 0 1408/2565 2197/4104 -1/5 0 (4th)
331 /// ```
332 ///
333 /// # References
334 /// - Fehlberg, E. (1969). "Low-order classical Runge-Kutta formulas with step size control and their application to some heat transfer problems".
335 pub fn rkf45() -> Self {
336 let mut c = [0.0; 6];
337 let mut a = [[0.0; 6]; 6];
338 let mut b = [0.0; 6];
339 let mut bh = [0.0; 6];
340
341 c[0] = 0.0;
342 c[1] = 1.0 / 4.0;
343 c[2] = 3.0 / 8.0;
344 c[3] = 12.0 / 13.0;
345 c[4] = 1.0;
346 c[5] = 1.0 / 2.0;
347
348 a[1][0] = 1.0 / 4.0;
349 a[2][0] = 3.0 / 32.0;
350 a[2][1] = 9.0 / 32.0;
351 a[3][0] = 1932.0 / 2197.0;
352 a[3][1] = -7200.0 / 2197.0;
353 a[3][2] = 7296.0 / 2197.0;
354 a[4][0] = 439.0 / 216.0;
355 a[4][1] = -8.0;
356 a[4][2] = 3680.0 / 513.0;
357 a[4][3] = -845.0 / 4104.0;
358 a[5][0] = -8.0 / 27.0;
359 a[5][1] = 2.0;
360 a[5][2] = -3544.0 / 2565.0;
361 a[5][3] = 1859.0 / 4104.0;
362 a[5][4] = -11.0 / 40.0;
363
364 b[0] = 16.0 / 135.0;
365 b[1] = 0.0;
366 b[2] = 6656.0 / 12825.0;
367 b[3] = 28561.0 / 56430.0;
368 b[4] = -9.0 / 50.0;
369 b[5] = 2.0 / 55.0;
370
371 bh[0] = 25.0 / 216.0;
372 bh[1] = 0.0;
373 bh[2] = 1408.0 / 2565.0;
374 bh[3] = 2197.0 / 4104.0;
375 bh[4] = -1.0 / 5.0;
376 bh[5] = 0.0;
377
378 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
379 let b = b.map(|x| T::from_f64(x).unwrap());
380 let bh = bh.map(|x| T::from_f64(x).unwrap());
381 let c = c.map(|x| T::from_f64(x).unwrap());
382
383 Self {
384 c,
385 a,
386 b,
387 bh: Some(bh),
388 bi: None,
389 er: None,
390 }
391 }
392
393 /// Cash-Karp 4(5) method.
394 ///
395 /// # Overview
396 /// This provides a 6-stage, explicit Runge-Kutta method with:
397 /// - Primary order: 5
398 /// - Embedded order: 4 (for error estimation)
399 /// - Number of stages: 6
400 ///
401 /// # Notes
402 /// - The Cash-Karp method is a variant of Runge-Kutta methods with embedded error estimation
403 /// that often provides better accuracy than RKF45 for some problem types.
404 ///
405 /// # Butcher Tableau
406 /// ```text
407 /// 0 |
408 /// 1/5 | 1/5
409 /// 3/10 | 3/40 9/40
410 /// 3/5 | 3/10 -9/10 6/5
411 /// 1 | -11/54 5/2 -70/27 35/27
412 /// 7/8 | 1631/55296 175/512 575/13824 44275/110592 253/4096
413 /// -------|---------------------------------------------------------------
414 /// | 37/378 0 250/621 125/594 0 512/1771 (5th)
415 /// | 2825/27648 0 18575/48384 13525/55296 277/14336 1/4 (4th)
416 /// ```
417 ///
418 /// # References
419 /// - Cash, J.R., Karp, A.H. (1990). "A Variable Order Runge-Kutta Method for Initial Value Problems with Rapidly Varying Right-Hand Sides".
420 pub fn cash_karp() -> Self {
421 let mut c = [0.0; 6];
422 let mut a = [[0.0; 6]; 6];
423 let mut b = [0.0; 6];
424 let mut bh = [0.0; 6];
425
426 c[0] = 0.0;
427 c[1] = 1.0 / 5.0;
428 c[2] = 3.0 / 10.0;
429 c[3] = 3.0 / 5.0;
430 c[4] = 1.0;
431 c[5] = 7.0 / 8.0;
432
433 a[1][0] = 1.0 / 5.0;
434 a[2][0] = 3.0 / 40.0;
435 a[2][1] = 9.0 / 40.0;
436 a[3][0] = 3.0 / 10.0;
437 a[3][1] = -9.0 / 10.0;
438 a[3][2] = 6.0 / 5.0;
439 a[4][0] = -11.0 / 54.0;
440 a[4][1] = 5.0 / 2.0;
441 a[4][2] = -70.0 / 27.0;
442 a[4][3] = 35.0 / 27.0;
443 a[5][0] = 1631.0 / 55296.0;
444 a[5][1] = 175.0 / 512.0;
445 a[5][2] = 575.0 / 13824.0;
446 a[5][3] = 44275.0 / 110592.0;
447 a[5][4] = 253.0 / 4096.0;
448
449 b[0] = 37.0 / 378.0;
450 b[1] = 0.0;
451 b[2] = 250.0 / 621.0;
452 b[3] = 125.0 / 594.0;
453 b[4] = 0.0;
454 b[5] = 512.0 / 1771.0;
455
456 bh[0] = 2825.0 / 27648.0;
457 bh[1] = 0.0;
458 bh[2] = 18575.0 / 48384.0;
459 bh[3] = 13525.0 / 55296.0;
460 bh[4] = 277.0 / 14336.0;
461 bh[5] = 1.0 / 4.0;
462
463 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
464 let b = b.map(|x| T::from_f64(x).unwrap());
465 let bh = bh.map(|x| T::from_f64(x).unwrap());
466 let c = c.map(|x| T::from_f64(x).unwrap());
467
468 Self {
469 c,
470 a,
471 b,
472 bh: Some(bh),
473 bi: None,
474 er: None,
475 }
476 }
477}
478
479// --- Implicit Runge-Kutta methods (IRK) ---
480
481impl<T: Real> ButcherTableau<T, 1> {
482 pub fn backward_euler() -> Self {
483 let mut c = [0.0; 1];
484 let mut a = [[0.0; 1]; 1];
485 let mut b = [0.0; 1];
486
487 c[0] = 1.0;
488 a[0][0] = 1.0;
489 b[0] = 1.0;
490
491 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
492 let b = b.map(|x| T::from_f64(x).unwrap());
493 let c = c.map(|x| T::from_f64(x).unwrap());
494
495 Self {
496 c,
497 a,
498 b,
499 bh: None,
500 bi: None,
501 er: None,
502 }
503 }
504}
505
506impl<T: Real> ButcherTableau<T, 2> {
507 pub fn trapezoidal() -> Self {
508 let mut c = [0.0; 2];
509 let mut a = [[0.0; 2]; 2];
510 let mut b = [0.0; 2];
511
512 c[0] = 0.0;
513 c[1] = 1.0;
514
515 a[1][0] = 1.0;
516
517 b[0] = 0.5;
518 b[1] = 0.5;
519
520 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
521 let b = b.map(|x| T::from_f64(x).unwrap());
522 let c = c.map(|x| T::from_f64(x).unwrap());
523
524 Self {
525 c,
526 a,
527 b,
528 bh: None,
529 bi: None,
530 er: None,
531 }
532 }
533
534 pub fn crank_nicolson() -> Self {
535 let mut c = [0.0; 2];
536 let mut a = [[0.0; 2]; 2];
537 let mut b = [0.0; 2];
538
539 c[0] = 0.0;
540 c[1] = 1.0;
541
542 a[1][0] = 1.0;
543
544 b[0] = 0.5;
545 b[1] = 0.5;
546
547 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
548 let b = b.map(|x| T::from_f64(x).unwrap());
549 let c = c.map(|x| T::from_f64(x).unwrap());
550
551 Self {
552 c,
553 a,
554 b,
555 bh: None,
556 bi: None,
557 er: None,
558 }
559 }
560}