differential_equations/tableau/
gauss_legendre.rs

1//! Gauss-Legendre Implicit Runge-Kutta Methods
2
3use crate::{tableau::ButcherTableau, traits::Real};
4
5const SQRT_3: f64 = 1.732050808;
6
7impl<T: Real> ButcherTableau<T, 2> {
8    /// Butcher Tableau for the Gauss-Legendre method of order 4.
9    ///
10    /// # Overview
11    /// This provides a 2-stage, implicit Runge-Kutta method (Gauss-Legendre) with:
12    /// - Primary order: 4
13    /// - Number of stages: 2
14    ///
15    /// # Notes
16    /// - Gauss-Legendre methods are A-stable and symmetric.
17    /// - They are highly accurate for their number of stages.
18    /// - The `c` values are the roots of the Legendre polynomial P_s(2x-1) = 0.
19    /// - This implementation includes two sets of `b` coefficients. The primary `b` coefficients
20    ///   are used for the solution, and `bh` can represent alternative coefficients (often related to error estimation or specific properties).
21    ///
22    /// # Butcher Tableau
23    /// ```text
24    /// (1/2 - sqrt(3)/6) |  1/4                  1/4 - sqrt(3)/6
25    /// (1/2 + sqrt(3)/6) |  1/4 + sqrt(3)/6      1/4
26    /// -------------------|---------------------------------------
27    ///                    |  1/2                  1/2
28    ///                    |  1/2 + sqrt(3)/2      1/2 - sqrt(3)/2  (bh coefficients)
29    /// ```
30    ///
31    /// # References
32    /// - Hairer, E., Nørsett, S. P., & Wanner, G. (1993). *Solving Ordinary Differential Equations I: Nonstiff Problems*. Springer. (Page 200, Table 4.5)
33    pub fn gauss_legendre_4() -> Self {
34        let mut c = [0.0; 2];
35        let mut a = [[0.0; 2]; 2];
36        let mut b = [0.0; 2];
37        let mut bh = [0.0; 2];
38
39        let sqrt3_6 = SQRT_3 / 6.0;
40        let sqrt3_2 = SQRT_3 / 2.0;
41
42        c[0] = 0.5 - sqrt3_6;
43        c[1] = 0.5 + sqrt3_6;
44
45        a[0][0] = 0.25;
46        a[0][1] = 0.25 - sqrt3_6;
47        a[1][0] = 0.25 + sqrt3_6;
48        a[1][1] = 0.25;
49
50        b[0] = 0.5;
51        b[1] = 0.5;
52
53        bh[0] = 0.5 + sqrt3_2;
54        bh[1] = 0.5 - sqrt3_2;
55
56        let c = c.map(|x| T::from_f64(x).unwrap());
57        let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
58        let b = b.map(|x| T::from_f64(x).unwrap());
59        let bh = bh.map(|x| T::from_f64(x).unwrap());
60
61        Self {
62            c,
63            a,
64            b,
65            bh: Some(bh),
66            bi: None,
67            er: None,
68        }
69    }
70}
71
72impl<T: Real> ButcherTableau<T, 3> {
73    /// Butcher Tableau for the Gauss-Legendre method of order 6.
74    ///
75    /// # Overview
76    /// This provides a 3-stage, implicit Runge-Kutta method (Gauss-Legendre) with:
77    /// - Primary order: 6
78    /// - Number of stages: 3
79    ///
80    /// # Notes
81    /// - Gauss-Legendre methods are A-stable and symmetric.
82    /// - They are highly accurate for their number of stages.
83    /// - The `c` values are the roots of the Legendre polynomial P_s(2x-1) = 0.
84    /// - This implementation includes two sets of `b` coefficients. The primary `b` coefficients
85    ///   are used for the solution, and `bh` can represent alternative coefficients.
86    ///
87    /// # Butcher Tableau
88    /// ```text
89    /// (1/2 - sqrt(15)/10) |  5/36                2/9 - sqrt(15)/15    5/36 - sqrt(15)/30
90    ///  1/2                |  5/36 + sqrt(15)/24  2/9                  5/36 - sqrt(15)/24
91    /// (1/2 + sqrt(15)/10) |  5/36 + sqrt(15)/30  2/9 + sqrt(15)/15    5/36
92    /// --------------------|-------------------------------------------------------------
93    ///                     |  5/18              4/9                  5/18
94    ///                     | -5/6               8/3                 -5/6                (bh coefficients)
95    /// ```
96    ///
97    /// # References
98    /// - Hairer, E., Nørsett, S. P., & Wanner, G. (1993). *Solving Ordinary Differential Equations I: Nonstiff Problems*. Springer. (Page 200, Table 4.5)
99    pub fn gauss_legendre_6() -> Self {
100        let mut c = [0.0; 3];
101        let mut a = [[0.0; 3]; 3];
102        let mut b = [0.0; 3];
103        let mut bh = [0.0; 3];
104
105        let sqrt_15: f64 = 3.872983346207417;
106        let sqrt15_10 = sqrt_15 / 10.0;
107        let sqrt15_15 = sqrt_15 / 15.0;
108        let sqrt15_24 = sqrt_15 / 24.0;
109        let sqrt15_30 = sqrt_15 / 30.0;
110
111        c[0] = 0.5 - sqrt15_10;
112        c[1] = 0.5;
113        c[2] = 0.5 + sqrt15_10;
114
115        a[0][0] = 5.0 / 36.0;
116        a[0][1] = 2.0 / 9.0 - sqrt15_15;
117        a[0][2] = 5.0 / 36.0 - sqrt15_30;
118
119        a[1][0] = 5.0 / 36.0 + sqrt15_24;
120        a[1][1] = 2.0 / 9.0;
121        a[1][2] = 5.0 / 36.0 - sqrt15_24;
122
123        a[2][0] = 5.0 / 36.0 + sqrt15_30;
124        a[2][1] = 2.0 / 9.0 + sqrt15_15;
125        a[2][2] = 5.0 / 36.0;
126
127        b[0] = 5.0 / 18.0;
128        b[1] = 4.0 / 9.0;
129        b[2] = 5.0 / 18.0;
130
131        bh[0] = -5.0 / 6.0;
132        bh[1] = 8.0 / 3.0;
133        bh[2] = -5.0 / 6.0;
134
135        let c = c.map(|x| T::from_f64(x).unwrap());
136        let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
137        let b = b.map(|x| T::from_f64(x).unwrap());
138        let bh = bh.map(|x| T::from_f64(x).unwrap());
139
140        Self {
141            c,
142            a,
143            b,
144            bh: Some(bh),
145            bi: None,
146            er: None,
147        }
148    }
149}
differential_equations/tableau/gauss_legendre.rs

differential_equations/tableau/
gauss_legendre.rs
