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
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
//! Diagonally Implicit Runge-Kutta (DIRK) tableau
use crate::{tableau::ButcherTableau, traits::Real};
impl<T: Real> ButcherTableau<T, 2> {
/// SDIRK-2-1: 2-stage, 2nd order SDIRK method with 1st order embedding
///
/// # Overview
/// This provides a 2-stage, singly diagonally implicit Runge-Kutta method with:
/// - Primary order: 2
/// - Embedded order: 1 (for error estimation)
/// - Number of stages: 2
/// - A-stable and B-stable
///
/// # Notes
/// - This is a simple SDIRK method where all diagonal entries are equal (γ = 1)
/// - Good for basic adaptive stepping with stiff problems
/// - The embedded method provides basic error estimation for step size control
/// - Particularly useful as a starting method for more complex stiff systems
///
/// # Butcher Tableau
/// ```text
/// 1 | 1 0
/// 0 | -1 1
/// -----|--------
/// | 1/2 1/2 (2nd order)
/// | 1 0 (1st order embedding)
/// ```
///
/// # References
/// - Hairer, E., Wanner, G. (1996). "Solving Ordinary Differential Equations II: Stiff and Differential-Algebraic Problems"
///
pub fn sdirk21() -> Self {
let c = [1.0, 0.0];
let a = [[1.0, 0.0], [-1.0, 1.0]];
let b = [0.5, 0.5];
let bh = [1.0, 0.0]; // 1st order embedding
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 = Some(bh.map(|x| T::from_f64(x).unwrap()));
ButcherTableau {
c,
a,
b,
bh,
bi: None,
er: None,
}
}
}
impl<T: Real> ButcherTableau<T, 3> {
/// ESDIRK-3-3: 3-stage, 3rd order ESDIRK method
///
/// # Overview
/// This provides a 3-stage, explicit singly diagonally implicit Runge-Kutta method with:
/// - Primary order: 3
/// - Number of stages: 3
/// - A-stable
///
/// # Notes
/// - This method pairs with SSPRK(3,3)-Shu-Osher-ERK to make a 3rd order IMEX method
/// - Has an explicit first stage (ESDIRK property) making it computationally efficient
/// - The first stage being explicit reduces the computational cost per step
/// - Suitable for problems with both stiff and non-stiff components
///
/// # Butcher Tableau
/// ```text
/// 0 | 0 0 0
/// 1 | 4γ+2β 1-4γ-2β 0
/// 1/2 | α₃₁ γ β
/// -------|------------------
/// | 1/6 1/6 2/3
/// ```
/// where:
/// - β = √3/6 + 1/2 ≈ 0.7886751346
/// - γ = (-1/8)(√3 + 1) ≈ -0.3416407865
/// - α₃₁ = 1/2 - β - γ ≈ 0.0529656519
///
/// # References
/// - Conde, S., et al. (2017). "Implicit and implicit-explicit strong stability preserving Runge-Kutta methods"
///
pub fn esdirk33() -> Self {
// Parameters (computed in f64 for precision)
let sqrt3 = 3.0_f64.sqrt();
let beta = sqrt3 / 6.0 + 0.5;
let gamma = (-1.0 / 8.0) * (sqrt3 + 1.0);
let c = [0.0, 1.0, 0.5];
let a = [
[0.0, 0.0, 0.0],
[
4.0 * gamma + 2.0 * beta,
1.0 - 4.0 * gamma - 2.0 * beta,
0.0,
],
[0.5 - beta - gamma, gamma, beta],
];
let b = [1.0 / 6.0, 1.0 / 6.0, 2.0 / 3.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());
ButcherTableau {
c,
a,
b,
bh: None,
bi: None,
er: None,
}
}
}
impl<T: Real> ButcherTableau<T, 4> {
/// ESDIRK3(2)4L[2]SA: 4-stage, 3rd order ESDIRK method with embedded 2nd order
///
/// # Overview
/// This provides a 4-stage, explicit singly diagonally implicit Runge-Kutta method with:
/// - Primary order: 3
/// - Embedded order: 2 (for error estimation)
/// - Number of stages: 4
/// - A-stable and B-stable
///
/// # References
/// - Kennedy, C.A. and Carpenter, M.H. (2003). "Additive Runge-Kutta schemes for convection-diffusion-reaction equations"
///
pub fn esdirk324l2sa() -> Self {
// Gamma parameter and derived values (computed in f64 for precision)
let g = 0.435_866_521_508_459;
let g2 = g * g;
let g3 = g2 * g;
let g4 = g3 * g;
let g5 = g4 * g;
let c3 = 3.0 / 5.0;
// Compute coefficients in f64
let a32 = c3 * (c3 - 2.0 * g) / (4.0 * g);
let a31 = c3 - g - a32;
let b2 = (-2.0 + 3.0 * c3 + 6.0 * g * (1.0 - c3)) / (12.0 * g * (c3 - 2.0 * g));
let b3 = (1.0 - 6.0 * g + 6.0 * g2) / (3.0 * c3 * (c3 - 2.0 * g));
let b1 = 1.0 - g - b2 - b3;
// Embedding coefficients
let d2_term1 = c3 * (-1.0 + 6.0 * g - 24.0 * g3 + 12.0 * g4 - 6.0 * g5)
/ (4.0 * g * (2.0 * g - c3) * (1.0 - 6.0 * g + 6.0 * g2));
let d2_term2 = (3.0 - 27.0 * g + 68.0 * g2 - 55.0 * g3 + 21.0 * g4 - 6.0 * g5)
/ (2.0 * (2.0 * g - c3) * (1.0 - 6.0 * g + 6.0 * g2));
let d2 = d2_term1 + d2_term2;
let d3 = -g * (-2.0 + 21.0 * g - 68.0 * g2 + 79.0 * g3 - 33.0 * g4 + 12.0 * g5)
/ (c3 * (c3 - 2.0 * g) * (1.0 - 6.0 * g + 6.0 * g2));
let d4 = -3.0 * g2 * (-1.0 + 4.0 * g - 2.0 * g2 + g3) / (1.0 - 6.0 * g + 6.0 * g2);
let d1 = 1.0 - d2 - d3 - d4;
let c = [0.0, 2.0 * g, 3.0 / 5.0, 1.0];
let a = [
[0.0, 0.0, 0.0, 0.0],
[g, g, 0.0, 0.0],
[a31, a32, g, 0.0],
[b1, b2, b3, g],
];
let b = [b1, b2, b3, g];
let bh = [d1, d2, d3, d4]; // 2nd order embedding
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 = Some(bh.map(|x| T::from_f64(x).unwrap()));
ButcherTableau {
c,
a,
b,
bh,
bi: None,
er: None,
}
}
}