differential_equations/tableau/
radau.rs

1//! Radau Implicit Runge-Kutta methods.
2
3use crate::{tableau::ButcherTableau, traits::Real};
4
5const SQRT_6: f64 = 2.449489743;
6
7impl<T: Real> ButcherTableau<T, 2> {
8    /// Butcher Tableau for the Radau IIA method of order 3.
9    ///
10    /// # Overview
11    /// This provides a 2-stage, implicit Runge-Kutta method (Radau IIA) with:
12    /// - Primary order: 3
13    /// - Embedded order: None (this implementation does not include an embedded error estimate)
14    /// - Number of stages: 2
15    ///
16    /// # Interpolation
17    /// - This standard Radau IIA order 3 implementation does not provide coefficients for dense output (interpolation).
18    ///
19    /// # Notes
20    /// - Radau IIA methods are A-stable and L-stable, making them suitable for stiff differential equations.
21    /// - Being implicit, they generally require solving a system of algebraic equations at each step.
22    ///
23    /// # Butcher Tableau
24    /// ```text
25    /// 1/3 |  5/12  -1/12
26    /// 1   |  3/4    1/4
27    /// ----|---------------
28    ///     |  3/4    1/4
29    /// ```
30    ///
31    /// # References
32    /// - Hairer, E., Nørsett, S. P., & Wanner, G. (1993). *Solving Ordinary Differential Equations I: Nonstiff Problems*. Springer.
33    /// - Hairer, E., & Wanner, G. (1996). *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems*. Springer.
34    pub fn radau_iia_3() -> Self {
35        let mut c = [0.0; 2];
36        let mut a = [[0.0; 2]; 2];
37        let mut b = [0.0; 2];
38
39        c[0] = 1.0 / 3.0;
40        c[1] = 1.0;
41
42        a[0][0] = 5.0 / 12.0;
43        a[0][1] = -1.0 / 12.0;
44        a[1][0] = 3.0 / 4.0;
45        a[1][1] = 1.0 / 4.0;
46
47        b[0] = 3.0 / 4.0;
48        b[1] = 1.0 / 4.0;
49
50        let c = c.map(|x| T::from_f64(x).unwrap());
51        let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
52        let b = b.map(|x| T::from_f64(x).unwrap());
53
54        Self {
55            c,
56            a,
57            b,
58            bh: None,
59            bi: None,
60            er: None,
61        }
62    }
63}
64
65impl<T: Real> ButcherTableau<T, 3> {
66    /// Butcher Tableau for the Radau IIA method of order 5.
67    ///
68    /// # Overview
69    /// This provides a 3-stage, implicit Runge-Kutta method (Radau IIA) with:
70    /// - Primary order: 5
71    /// - Embedded order: None (this implementation does not include an embedded error estimate)
72    /// - Number of stages: 3
73    ///
74    /// # Interpolation
75    /// - This standard Radau IIA order 5 implementation does not provide coefficients for dense output (interpolation).
76    ///
77    /// # Notes
78    /// - Radau IIA methods are A-stable and L-stable, making them highly suitable for stiff differential equations.
79    /// - Being implicit, they generally require solving a system of algebraic equations at each step.
80    /// - The `c` values are the roots of `P_s(2x-1) - P_{s-1}(2x-1) = 0` where `P_s` is the s-th Legendre polynomial.
81    ///   For s=3, the c values are `(2/5 - sqrt(6)/10)`, `(2/5 + sqrt(6)/10)`, and `1`.
82    ///
83    /// # Butcher Tableau
84    /// ```text
85    /// c1 | a11 a12 a13
86    /// c2 | a21 a22 a23
87    /// c3 | a31 a32 a33
88    /// ---|------------
89    ///    | b1  b2  b3
90    /// ```
91    /// Where:
92    /// c1 = 2/5 - sqrt(6)/10
93    /// c2 = 2/5 + sqrt(6)/10
94    /// c3 = 1
95    ///
96    /// a11 = 11/45 - 7*sqrt(6)/360
97    /// a12 = 37/225 - 169*sqrt(6)/1800
98    /// a13 = -2/225 + sqrt(6)/75
99    ///
100    /// a21 = 37/225 + 169*sqrt(6)/1800
101    /// a22 = 11/45 + 7*sqrt(6)/360
102    /// a23 = -2/225 - sqrt(6)/75
103    ///
104    /// a31 = 4/9 - sqrt(6)/36
105    /// a32 = 4/9 + sqrt(6)/36
106    /// a33 = 1/9
107    ///
108    /// b1 = 4/9 - sqrt(6)/36
109    /// b2 = 4/9 + sqrt(6)/36
110    /// b3 = 1/9
111    ///
112    /// # References
113    /// - Hairer, E., Nørsett, S. P., & Wanner, G. (1993). *Solving Ordinary Differential Equations I: Nonstiff Problems*. Springer.
114    /// - Hairer, E., & Wanner, G. (1996). *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems*. Springer.
115    pub fn radau_iia_5() -> Self {
116        let mut c = [0.0; 3];
117        let mut a = [[0.0; 3]; 3];
118        let mut b = [0.0; 3];
119
120        c[0] = 2.0 / 5.0 - SQRT_6 / 10.0;
121        c[1] = 2.0 / 5.0 + SQRT_6 / 10.0;
122        c[2] = 1.0;
123
124        a[0][0] = 11.0 / 45.0 - 7.0 * SQRT_6 / 360.0;
125        a[0][1] = 37.0 / 225.0 - 169.0 * SQRT_6 / 1800.0;
126        a[0][2] = -2.0 / 225.0 + SQRT_6 / 75.0;
127
128        a[1][0] = 37.0 / 225.0 + 169.0 * SQRT_6 / 1800.0;
129        a[1][1] = 11.0 / 45.0 + 7.0 * SQRT_6 / 360.0;
130        a[1][2] = -2.0 / 225.0 - SQRT_6 / 75.0;
131
132        a[2][0] = 4.0 / 9.0 - SQRT_6 / 36.0;
133        a[2][1] = 4.0 / 9.0 + SQRT_6 / 36.0;
134        a[2][2] = 1.0 / 9.0;
135
136        b[0] = 4.0 / 9.0 - SQRT_6 / 36.0;
137        b[1] = 4.0 / 9.0 + SQRT_6 / 36.0;
138        b[2] = 1.0 / 9.0;
139
140        let c = c.map(|x| T::from_f64(x).unwrap());
141        let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
142        let b = b.map(|x| T::from_f64(x).unwrap());
143
144        Self {
145            c,
146            a,
147            b,
148            bh: None,
149            bi: None,
150            er: None,
151        }
152    }
153}
differential_equations/tableau/radau.rs

differential_equations/tableau/
radau.rs
