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