differential_equations/tableau/kvaerno.rs
1//! Kvaerno Diagonally Implicit Runge-Kutta (DIRK) methods
2
3use crate::{tableau::ButcherTableau, traits::Real};
4
5impl<T: Real> ButcherTableau<T, 4> {
6 /// Kvaerno(4,2,3): 4-stage, 3rd order DIRK method with embedded 2nd order
7 ///
8 /// # Overview
9 /// This provides a 4-stage, diagonally implicit Runge-Kutta method with:
10 /// - Primary order: 3
11 /// - Embedded order: 2 (for error estimation)
12 /// - Number of stages: 4
13 /// - A-stable
14 ///
15 /// # Notes
16 /// - Developed specifically for stiff differential equations
17 /// - A-stable method suitable for moderately stiff problems
18 /// - The embedded method provides efficient error estimation for adaptive stepping
19 /// - All diagonal elements are identical (γ ≈ 0.4358665215), simplifying LU factorization
20 /// - Good balance between stability and computational efficiency
21 ///
22 /// # Butcher Tableau
23 /// ```text
24 /// 0 | 0 0 0 0
25 /// .8717 | .4359 .4359 0 0
26 /// 1 | .4906 .0736 .4359 0
27 /// 1 | .3088 1.4906 -1.2352 .4359
28 /// -------|-------------------------
29 /// b³ | .3088 1.4906 -1.2352 .4359
30 /// b² | .4906 .0736 .4359 0
31 /// ```
32 /// where γ ≈ 0.4358665215
33 ///
34 /// # References
35 /// - Kvaerno, A. (2004). "Singly diagonally implicit Runge-Kutta methods with an explicit first stage"
36 ///
37 pub fn kvaerno423() -> Self {
38 // Main diagonal entry
39 let gamma = 0.4358665215;
40
41 let c = [0.0, 0.871733043, 1.0, 1.0];
42
43 let a = [
44 [0.0, 0.0, 0.0, 0.0],
45 [gamma, gamma, 0.0, 0.0],
46 [0.490563388419108, 0.073570090080892, gamma, 0.0],
47 [
48 0.308809969973036,
49 1.490563388254106,
50 -1.235239879727145,
51 gamma,
52 ],
53 ];
54
55 let b = [
56 0.308809969973036,
57 1.490563388254106,
58 -1.235239879727145,
59 gamma,
60 ];
61
62 let bh = [0.490563388419108, 0.073570090080892, gamma, 0.0];
63
64 let c = c.map(|x| T::from_f64(x).unwrap());
65 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
66 let b = b.map(|x| T::from_f64(x).unwrap());
67 let bh = Some(bh.map(|x| T::from_f64(x).unwrap()));
68
69 ButcherTableau {
70 c,
71 a,
72 b,
73 bh,
74 bi: None,
75 er: None,
76 }
77 }
78}
79
80impl<T: Real> ButcherTableau<T, 7> {
81 /// Kvaerno(7,4,5): 7-stage, 5th order DIRK method with embedded 4th order
82 ///
83 /// # Overview
84 /// This provides a 7-stage, diagonally implicit Runge-Kutta method with:
85 /// - Primary order: 5
86 /// - Embedded order: 4 (for error estimation)
87 /// - Number of stages: 7
88 /// - A-stable and B-stable
89 ///
90 /// # Notes
91 /// - High-order method designed for stiff differential equations requiring high accuracy
92 /// - A-stable and B-stable for excellent stability properties
93 /// - The embedded 4th order method provides accurate error estimation for adaptive control
94 /// - All diagonal elements are identical (γ = 0.26), enabling efficient LU factorization reuse
95 /// - Particularly effective for smooth solutions where high accuracy is required
96 /// - Higher computational cost per step but fewer steps needed due to high order
97 ///
98 /// # Butcher Tableau
99 /// ```text
100 /// 0 | 0 0 0 0 0 0 0
101 /// .52 | .26 .26 0 0 0 0 0
102 /// 1.2303 | .13 .8403 .26 0 0 0 0
103 /// .8958 | .2237 .4768 -.0647 .26 0 0 0
104 /// .4364 | .1665 .1045 .0363 -.1309 .26 0 0
105 /// 1 | .1386 .0000 -.0425 .0245 .6194 .26 0
106 /// 1 | .1366 .0000 -.0550 -.0412 .6299 .0696 .26
107 /// -------|--------------------------------------------
108 /// b⁵ | .1366 .0000 -.0550 -.0412 .6299 .0696 .26
109 /// b⁴ | .1386 .0000 -.0425 .0245 .6194 .26 0
110 /// ```
111 /// where γ = 0.26 exactly
112 ///
113 /// # References
114 /// - Kvaerno, A. (2004). "Singly diagonally implicit Runge-Kutta methods with an explicit first stage"
115 ///
116 pub fn kvaerno745() -> Self {
117 // Main diagonal entry
118 let gamma = 0.26;
119
120 let c = [
121 0.0,
122 0.52,
123 1.230333209967908,
124 0.895765984350076,
125 0.436393609858648,
126 1.0,
127 1.0,
128 ];
129
130 let a = [
131 [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
132 [gamma, gamma, 0.0, 0.0, 0.0, 0.0, 0.0],
133 [0.13, 0.840_333_209_967_908_1, gamma, 0.0, 0.0, 0.0, 0.0],
134 [
135 0.22371961478320505,
136 0.476_755_323_197_997,
137 -0.06470895363112615,
138 gamma,
139 0.0,
140 0.0,
141 0.0,
142 ],
143 [
144 0.16648564323248321,
145 0.104_500_188_415_917_2,
146 0.03631482272098715,
147 -0.13090704451073998,
148 gamma,
149 0.0,
150 0.0,
151 ],
152 [
153 0.13855640231268224,
154 0.0,
155 -0.04245337201752043,
156 0.02446657898003141,
157 0.619_430_390_724_806_8,
158 gamma,
159 0.0,
160 ],
161 [
162 0.13659751177640291,
163 0.0,
164 -0.05496908796538376,
165 -0.04118626728321046,
166 0.629_933_048_990_164,
167 0.06962479448202728,
168 gamma,
169 ],
170 ];
171
172 let b = [
173 0.13659751177640291,
174 0.0,
175 -0.05496908796538376,
176 -0.04118626728321046,
177 0.629_933_048_990_164,
178 0.06962479448202728,
179 gamma,
180 ];
181
182 let bh = [
183 0.13855640231268224,
184 0.0,
185 -0.04245337201752043,
186 0.02446657898003141,
187 0.619_430_390_724_806_8,
188 gamma,
189 0.0,
190 ];
191
192 let c = c.map(|x| T::from_f64(x).unwrap());
193 let a = a.map(|row| row.map(|x| T::from_f64(x).unwrap()));
194 let b = b.map(|x| T::from_f64(x).unwrap());
195 let bh = Some(bh.map(|x| T::from_f64(x).unwrap()));
196
197 ButcherTableau {
198 c,
199 a,
200 b,
201 bh,
202 bi: None,
203 er: None,
204 }
205 }
206}