differential_equations/tableau/
lobatto.rs

1//! Lobatto Implicit Runge-Kutta Methods
2
3use crate::{tableau::ButcherTableau, traits::Real};
4
5impl<T: Real> ButcherTableau<T, 2> {
6    /// Butcher Tableau for the Lobatto IIIC method of order 2.
7    ///
8    /// # Overview
9    /// This provides a 2-stage, implicit Runge-Kutta method (Lobatto IIIC) with:
10    /// - Primary order: 2
11    /// - Number of stages: 2
12    ///
13    /// # Notes
14    /// - Lobatto IIIC methods are L-stable.
15    /// - They are algebraically stable and thus B-stable, making them suitable for stiff problems.
16    /// - This implementation includes two sets of `b` coefficients. The primary `b` coefficients
17    ///   are used for the solution, and `bh` can represent alternative coefficients if needed
18    ///   (though their specific use here as a second row is characteristic of Lobatto IIIC).
19    ///
20    /// # Butcher Tableau
21    /// ```text
22    /// 0   |  1/2  -1/2
23    /// 1   |  1/2   1/2
24    /// ----|------------
25    ///     |  1/2   1/2   (b coefficients)
26    ///     |  1     0     (bh coefficients)
27    /// ```
28    ///
29    /// # References
30    /// - Hairer, E., & Wanner, G. (1996). *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems*. Springer. (Page 80, Table 5.7)
31    pub fn lobatto_iiic_2() -> Self {
32        let mut c = [0.0; 2];
33        let mut a = [[0.0; 2]; 2];
34        let mut b = [0.0; 2];
35        let mut bh = [0.0; 2];
36
37        c[0] = 0.0;
38        c[1] = 1.0;
39
40        a[0][0] = 1.0 / 2.0;
41        a[0][1] = -1.0 / 2.0;
42        a[1][0] = 1.0 / 2.0;
43        a[1][1] = 1.0 / 2.0;
44
45        b[0] = 1.0 / 2.0;
46        b[1] = 1.0 / 2.0;
47
48        bh[0] = 1.0;
49        bh[1] = 0.0;
50
51        let c = c.map(|x| T::from_f64(x).unwrap());
52        let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
53        let b = b.map(|x| T::from_f64(x).unwrap());
54        let bh = bh.map(|x| T::from_f64(x).unwrap());
55
56        Self {
57            c,
58            a,
59            b,
60            bh: Some(bh),
61            bi: None,
62            er: None,
63        }
64    }
65}
66
67impl<T: Real> ButcherTableau<T, 3> {
68    /// Butcher Tableau for the Lobatto IIIC method of order 4.
69    ///
70    /// # Overview
71    /// This provides a 3-stage, implicit Runge-Kutta method (Lobatto IIIC) with:
72    /// - Primary order: 4
73    /// - Number of stages: 3
74    ///
75    /// # Notes
76    /// - Lobatto IIIC methods are L-stable.
77    /// - They are algebraically stable and thus B-stable, making them suitable for stiff problems.
78    /// - This implementation includes two sets of `b` coefficients. The primary `b` coefficients
79    ///   are used for the solution, and `bh` represents the second row of coefficients from the tableau.
80    ///
81    /// # Butcher Tableau
82    /// ```text
83    /// 0   |  1/6  -1/3   1/6
84    /// 1/2 |  1/6   5/12 -1/12
85    /// 1   |  1/6   2/3   1/6
86    /// ----|------------------
87    ///     |  1/6   2/3   1/6   (b coefficients)
88    ///     | -1/2   2    -1/2   (bh coefficients)
89    /// ```
90    ///
91    /// # References
92    /// - Hairer, E., & Wanner, G. (1996). *Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems*. Springer. (Page 80, Table 5.7)
93    pub fn lobatto_iiic_4() -> Self {
94        let mut c = [0.0; 3];
95        let mut a = [[0.0; 3]; 3];
96        let mut b = [0.0; 3];
97        let mut bh = [0.0; 3];
98
99        c[0] = 0.0;
100        c[1] = 1.0 / 2.0;
101        c[2] = 1.0;
102
103        a[0][0] = 1.0 / 6.0;
104        a[0][1] = -1.0 / 3.0;
105        a[0][2] = 1.0 / 6.0;
106
107        a[1][0] = 1.0 / 6.0;
108        a[1][1] = 5.0 / 12.0;
109        a[1][2] = -1.0 / 12.0;
110
111        a[2][0] = 1.0 / 6.0;
112        a[2][1] = 2.0 / 3.0;
113        a[2][2] = 1.0 / 6.0;
114
115        b[0] = 1.0 / 6.0;
116        b[1] = 2.0 / 3.0;
117        b[2] = 1.0 / 6.0;
118
119        bh[0] = -1.0 / 2.0;
120        bh[1] = 2.0;
121        bh[2] = -1.0 / 2.0;
122
123        let c = c.map(|x| T::from_f64(x).unwrap());
124        let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
125        let b = b.map(|x| T::from_f64(x).unwrap());
126        let bh = bh.map(|x| T::from_f64(x).unwrap());
127
128        Self {
129            c,
130            a,
131            b,
132            bh: Some(bh),
133            bi: None,
134            er: None,
135        }
136    }
137}