conspire/math/integrate/verner_9/
mod.rs1#[cfg(test)]
2mod test;
3
4use super::{
5 super::{Scalar, Tensor, TensorVec, Vector, interpolate::InterpolateSolution},
6 Explicit, IntegrationError, OdeSolver,
7};
8use crate::{ABS_TOL, REL_TOL};
9use std::ops::{Mul, Sub};
10
11const C_2: Scalar = 0.03462;
12const C_3: Scalar = 0.097_024_350_638_780_44;
13const C_4: Scalar = 0.145_536_525_958_170_67;
14const C_5: Scalar = 0.561;
15const C_6: Scalar = 0.229_007_911_590_485;
16const C_7: Scalar = 0.544_992_088_409_515;
17const C_8: Scalar = 0.645;
18const C_9: Scalar = 0.48375;
19const C_10: Scalar = 0.06757;
20const C_11: Scalar = 0.2500;
21const C_12: Scalar = 0.659_065_061_873_099_9;
22const C_13: Scalar = 0.8206;
23const C_14: Scalar = 0.9012;
24
25const A_2_1: Scalar = 0.03462;
26const A_3_1: Scalar = -0.03893354388572875;
27const A_3_2: Scalar = 0.13595789452450918;
28const A_4_1: Scalar = 0.03638413148954267;
29const A_4_3: Scalar = 0.10915239446862801;
30const A_5_1: Scalar = 2.0257639143939694;
31const A_5_3: Scalar = -7.638023836496291;
32const A_5_4: Scalar = 6.173259922102322;
33const A_6_1: Scalar = 0.05112275589406061;
34const A_6_4: Scalar = 0.17708237945550218;
35const A_6_5: Scalar = 0.0008027762409222536;
36const A_7_1: Scalar = 0.13160063579752163;
37const A_7_4: Scalar = -0.2957276252669636;
38const A_7_5: Scalar = 0.08781378035642955;
39const A_7_6: Scalar = 0.6213052975225274;
40const A_8_1: Scalar = 0.07166666666666667;
41const A_8_6: Scalar = 0.33055335789153195;
42const A_8_7: Scalar = 0.2427799754418014;
43const A_9_1: Scalar = 0.071806640625;
44const A_9_6: Scalar = 0.3294380283228177;
45const A_9_7: Scalar = 0.1165190029271823;
46const A_9_8: Scalar = -0.034013671875;
47const A_10_1: Scalar = 0.04836757646340646;
48const A_10_6: Scalar = 0.03928989925676164;
49const A_10_7: Scalar = 0.10547409458903446;
50const A_10_8: Scalar = -0.021438652846483126;
51const A_10_9: Scalar = -0.10412291746271944;
52const A_11_1: Scalar = -0.026645614872014785;
53const A_11_6: Scalar = 0.03333333333333333;
54const A_11_7: Scalar = -0.1631072244872467;
55const A_11_8: Scalar = 0.03396081684127761;
56const A_11_9: Scalar = 0.1572319413814626;
57const A_11_10: Scalar = 0.21522674780318796;
58const A_12_1: Scalar = 0.03689009248708622;
59const A_12_6: Scalar = -0.1465181576725543;
60const A_12_7: Scalar = 0.2242577768172024;
61const A_12_8: Scalar = 0.02294405717066073;
62const A_12_9: Scalar = -0.0035850052905728597;
63const A_12_10: Scalar = 0.08669223316444385;
64const A_12_11: Scalar = 0.43838406519683376;
65const A_13_1: Scalar = -0.4866012215113341;
66const A_13_6: Scalar = -6.304602650282853;
67const A_13_7: Scalar = -0.2812456182894729;
68const A_13_8: Scalar = -2.679019236219849;
69const A_13_9: Scalar = 0.5188156639241577;
70const A_13_10: Scalar = 1.3653531876033418;
71const A_13_11: Scalar = 5.8850910885039465;
72const A_13_12: Scalar = 2.8028087862720628;
73const A_14_1: Scalar = 0.4185367457753472;
74const A_14_6: Scalar = 6.724547581906459;
75const A_14_7: Scalar = -0.42544428016461133;
76const A_14_8: Scalar = 3.3432791530012653;
77const A_14_9: Scalar = 0.6170816631175374;
78const A_14_10: Scalar = -0.9299661239399329;
79const A_14_11: Scalar = -6.099948804751011;
80const A_14_12: Scalar = -3.002206187889399;
81const A_14_13: Scalar = 0.2553202529443446;
82const A_15_1: Scalar = -0.7793740861228848;
83const A_15_6: Scalar = -13.937342538107776;
84const A_15_7: Scalar = 1.2520488533793563;
85const A_15_8: Scalar = -14.691500408016868;
86const A_15_9: Scalar = -0.494705058533141;
87const A_15_10: Scalar = 2.2429749091462368;
88const A_15_11: Scalar = 13.367893803828643;
89const A_15_12: Scalar = 14.396650486650687;
90const A_15_13: Scalar = -0.79758133317768;
91const A_15_14: Scalar = 0.4409353709534278;
92const A_16_1: Scalar = 2.0580513374668867;
93const A_16_6: Scalar = 22.357937727968032;
94const A_16_7: Scalar = 0.9094981099755646;
95const A_16_8: Scalar = 35.89110098240264;
96const A_16_9: Scalar = -3.442515027624454;
97const A_16_10: Scalar = -4.865481358036369;
98const A_16_11: Scalar = -18.909803813543427;
99const A_16_12: Scalar = -34.26354448030452;
100const A_16_13: Scalar = 1.2647565216956427;
101
102const B_1: Scalar = 0.014611976858423152;
103const B_8: Scalar = -0.3915211862331339;
104const B_9: Scalar = 0.23109325002895065;
105const B_10: Scalar = 0.12747667699928525;
106const B_11: Scalar = 0.2246434176204158;
107const B_12: Scalar = 0.5684352689748513;
108const B_13: Scalar = 0.058258715572158275;
109const B_14: Scalar = 0.13643174034822156;
110const B_15: Scalar = 0.030570139830827976;
111
112const D_1: Scalar = -0.005357988290444578;
113const D_8: Scalar = -2.583020491182464;
114const D_9: Scalar = 0.14252253154686625;
115const D_10: Scalar = 0.013420653512688676;
116const D_11: Scalar = -0.02867296291409493;
117const D_12: Scalar = 2.624999655215792;
118const D_13: Scalar = -0.2825509643291537;
119const D_14: Scalar = 0.13643174034822156;
120const D_15: Scalar = 0.030570139830827976;
121const D_16: Scalar = -0.04834231373823958;
122
123#[doc = include_str!("doc.md")]
124#[derive(Debug)]
125pub struct Verner9 {
126 pub abs_tol: Scalar,
128 pub rel_tol: Scalar,
130 pub dt_beta: Scalar,
132 pub dt_expn: Scalar,
134 pub dt_cut: Scalar,
136 pub dt_min: Scalar,
138}
139
140impl Default for Verner9 {
141 fn default() -> Self {
142 Self {
143 abs_tol: ABS_TOL,
144 rel_tol: REL_TOL,
145 dt_beta: 0.9,
146 dt_expn: 9.0,
147 dt_cut: 0.5,
148 dt_min: ABS_TOL,
149 }
150 }
151}
152
153impl<Y, U> OdeSolver<Y, U> for Verner9
154where
155 Y: Tensor,
156 U: TensorVec<Item = Y>,
157{
158 fn abs_tol(&self) -> Scalar {
159 self.abs_tol
160 }
161 fn dt_cut(&self) -> Scalar {
162 self.dt_cut
163 }
164 fn dt_min(&self) -> Scalar {
165 self.dt_min
166 }
167}
168
169impl<Y, U> Explicit<Y, U> for Verner9
170where
171 Self: OdeSolver<Y, U>,
172 Y: Tensor,
173 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
174 U: TensorVec<Item = Y>,
175{
176 const SLOPES: usize = 16;
177 fn dt_beta(&self) -> Scalar {
178 self.dt_beta
179 }
180 fn dt_expn(&self) -> Scalar {
181 self.dt_expn
182 }
183 fn slopes(
184 &self,
185 mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
186 y: &Y,
187 t: Scalar,
188 dt: Scalar,
189 k: &mut [Y],
190 y_trial: &mut Y,
191 ) -> Result<Scalar, String> {
192 k[0] = function(t, y)?;
193 *y_trial = &k[0] * (A_2_1 * dt) + y;
194 k[1] = function(t + C_2 * dt, y_trial)?;
195 *y_trial = &k[0] * (A_3_1 * dt) + &k[1] * (A_3_2 * dt) + y;
196 k[2] = function(t + C_3 * dt, y_trial)?;
197 *y_trial = &k[0] * (A_4_1 * dt) + &k[2] * (A_4_3 * dt) + y;
198 k[3] = function(t + C_4 * dt, y_trial)?;
199 *y_trial = &k[0] * (A_5_1 * dt) + &k[2] * (A_5_3 * dt) + &k[3] * (A_5_4 * dt) + y;
200 k[4] = function(t + C_5 * dt, y_trial)?;
201 *y_trial = &k[0] * (A_6_1 * dt) + &k[3] * (A_6_4 * dt) + &k[4] * (A_6_5 * dt) + y;
202 k[5] = function(t + C_6 * dt, y_trial)?;
203 *y_trial = &k[0] * (A_7_1 * dt)
204 + &k[3] * (A_7_4 * dt)
205 + &k[4] * (A_7_5 * dt)
206 + &k[5] * (A_7_6 * dt)
207 + y;
208 k[6] = function(t + C_7 * dt, y_trial)?;
209 *y_trial = &k[0] * (A_8_1 * dt) + &k[5] * (A_8_6 * dt) + &k[6] * (A_8_7 * dt) + y;
210 k[7] = function(t + C_8 * dt, y_trial)?;
211 *y_trial = &k[0] * (A_9_1 * dt)
212 + &k[5] * (A_9_6 * dt)
213 + &k[6] * (A_9_7 * dt)
214 + &k[7] * (A_9_8 * dt)
215 + y;
216 k[8] = function(t + C_9 * dt, y_trial)?;
217 *y_trial = &k[0] * (A_10_1 * dt)
218 + &k[5] * (A_10_6 * dt)
219 + &k[6] * (A_10_7 * dt)
220 + &k[7] * (A_10_8 * dt)
221 + &k[8] * (A_10_9 * dt)
222 + y;
223 k[9] = function(t + C_10 * dt, y_trial)?;
224 *y_trial = &k[0] * (A_11_1 * dt)
225 + &k[5] * (A_11_6 * dt)
226 + &k[6] * (A_11_7 * dt)
227 + &k[7] * (A_11_8 * dt)
228 + &k[8] * (A_11_9 * dt)
229 + &k[9] * (A_11_10 * dt)
230 + y;
231 k[10] = function(t + C_11 * dt, y_trial)?;
232 *y_trial = &k[0] * (A_12_1 * dt)
233 + &k[5] * (A_12_6 * dt)
234 + &k[6] * (A_12_7 * dt)
235 + &k[7] * (A_12_8 * dt)
236 + &k[8] * (A_12_9 * dt)
237 + &k[9] * (A_12_10 * dt)
238 + &k[10] * (A_12_11 * dt)
239 + y;
240 k[11] = function(t + C_12 * dt, y_trial)?;
241 *y_trial = &k[0] * (A_13_1 * dt)
242 + &k[5] * (A_13_6 * dt)
243 + &k[6] * (A_13_7 * dt)
244 + &k[7] * (A_13_8 * dt)
245 + &k[8] * (A_13_9 * dt)
246 + &k[9] * (A_13_10 * dt)
247 + &k[10] * (A_13_11 * dt)
248 + &k[11] * (A_13_12 * dt)
249 + y;
250 k[12] = function(t + C_13 * dt, y_trial)?;
251 *y_trial = &k[0] * (A_14_1 * dt)
252 + &k[5] * (A_14_6 * dt)
253 + &k[6] * (A_14_7 * dt)
254 + &k[7] * (A_14_8 * dt)
255 + &k[8] * (A_14_9 * dt)
256 + &k[9] * (A_14_10 * dt)
257 + &k[10] * (A_14_11 * dt)
258 + &k[11] * (A_14_12 * dt)
259 + &k[12] * (A_14_13 * dt)
260 + y;
261 k[13] = function(t + C_14 * dt, y_trial)?;
262 *y_trial = &k[0] * (A_15_1 * dt)
263 + &k[5] * (A_15_6 * dt)
264 + &k[6] * (A_15_7 * dt)
265 + &k[7] * (A_15_8 * dt)
266 + &k[8] * (A_15_9 * dt)
267 + &k[9] * (A_15_10 * dt)
268 + &k[10] * (A_15_11 * dt)
269 + &k[11] * (A_15_12 * dt)
270 + &k[12] * (A_15_13 * dt)
271 + &k[13] * (A_15_14 * dt)
272 + y;
273 k[14] = function(t + dt, y_trial)?;
274 *y_trial = &k[0] * (A_16_1 * dt)
275 + &k[5] * (A_16_6 * dt)
276 + &k[6] * (A_16_7 * dt)
277 + &k[7] * (A_16_8 * dt)
278 + &k[8] * (A_16_9 * dt)
279 + &k[9] * (A_16_10 * dt)
280 + &k[10] * (A_16_11 * dt)
281 + &k[11] * (A_16_12 * dt)
282 + &k[12] * (A_16_13 * dt)
283 + y;
284 k[15] = function(t + dt, y_trial)?;
285 *y_trial = (&k[0] * B_1
286 + &k[7] * B_8
287 + &k[8] * B_9
288 + &k[9] * B_10
289 + &k[10] * B_11
290 + &k[11] * B_12
291 + &k[12] * B_13
292 + &k[13] * B_14
293 + &k[14] * B_15)
294 * dt
295 + y;
296 Ok(((&k[0] * D_1
297 + &k[7] * D_8
298 + &k[8] * D_9
299 + &k[9] * D_10
300 + &k[10] * D_11
301 + &k[11] * D_12
302 + &k[12] * D_13
303 + &k[13] * D_14
304 + &k[14] * D_15
305 + &k[15] * D_16)
306 * dt)
307 .norm_inf())
308 }
309 fn step(
310 &self,
311 mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
312 y: &mut Y,
313 t: &mut Scalar,
314 y_sol: &mut U,
315 t_sol: &mut Vector,
316 dydt_sol: &mut U,
317 dt: &mut Scalar,
318 _k: &mut [Y],
319 y_trial: &Y,
320 e: Scalar,
321 ) -> Result<(), String> {
322 if e < self.abs_tol || e / y_trial.norm_inf() < self.rel_tol {
323 *t += *dt;
324 *y = y_trial.clone();
325 t_sol.push(*t);
326 y_sol.push(y.clone());
327 dydt_sol.push(function(*t, y)?);
328 }
329 self.time_step(e, dt);
330 Ok(())
331 }
332}
333
334impl<Y, U> InterpolateSolution<Y, U> for Verner9
335where
336 Y: Tensor,
337 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
338 U: TensorVec<Item = Y>,
339{
340 fn interpolate(
341 &self,
342 time: &Vector,
343 tp: &Vector,
344 yp: &U,
345 mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
346 ) -> Result<(U, U), IntegrationError> {
347 let mut dt;
348 let mut i;
349 let mut k_1;
350 let mut k_2;
351 let mut k_3;
352 let mut k_4;
353 let mut k_5;
354 let mut k_6;
355 let mut k_7;
356 let mut k_8;
357 let mut k_9;
358 let mut k_10;
359 let mut k_11;
360 let mut k_12;
361 let mut k_13;
362 let mut k_14;
363 let mut k_15;
364 let mut t;
365 let mut y;
366 let mut y_int = U::new();
367 let mut dydt_int = U::new();
368 let mut y_trial;
369 for time_k in time.iter() {
370 i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
371 if time_k == &tp[i] {
372 t = tp[i];
373 y_trial = yp[i].clone();
374 dt = 0.0;
375 } else {
376 t = tp[i - 1];
377 y = yp[i - 1].clone();
378 dt = time_k - t;
379 k_1 = function(t, &y)?;
380 y_trial = &k_1 * (A_2_1 * dt) + &y;
381 k_2 = function(t + C_2 * dt, &y_trial)?;
382 y_trial = &k_1 * (A_3_1 * dt) + &k_2 * (A_3_2 * dt) + &y;
383 k_3 = function(t + C_3 * dt, &y_trial)?;
384 y_trial = &k_1 * (A_4_1 * dt) + &k_3 * (A_4_3 * dt) + &y;
385 k_4 = function(t + C_4 * dt, &y_trial)?;
386 y_trial = &k_1 * (A_5_1 * dt) + &k_3 * (A_5_3 * dt) + &k_4 * (A_5_4 * dt) + &y;
387 k_5 = function(t + C_5 * dt, &y_trial)?;
388 y_trial = &k_1 * (A_6_1 * dt) + &k_4 * (A_6_4 * dt) + &k_5 * (A_6_5 * dt) + &y;
389 k_6 = function(t + C_6 * dt, &y_trial)?;
390 y_trial = &k_1 * (A_7_1 * dt)
391 + &k_4 * (A_7_4 * dt)
392 + &k_5 * (A_7_5 * dt)
393 + &k_6 * (A_7_6 * dt)
394 + &y;
395 k_7 = function(t + C_7 * dt, &y_trial)?;
396 y_trial = &k_1 * (A_8_1 * dt) + &k_6 * (A_8_6 * dt) + &k_7 * (A_8_7 * dt) + &y;
397 k_8 = function(t + C_8 * dt, &y_trial)?;
398 y_trial = &k_1 * (A_9_1 * dt)
399 + &k_6 * (A_9_6 * dt)
400 + &k_7 * (A_9_7 * dt)
401 + &k_8 * (A_9_8 * dt)
402 + &y;
403 k_9 = function(t + C_9 * dt, &y_trial)?;
404 y_trial = &k_1 * (A_10_1 * dt)
405 + &k_6 * (A_10_6 * dt)
406 + &k_7 * (A_10_7 * dt)
407 + &k_8 * (A_10_8 * dt)
408 + &k_9 * (A_10_9 * dt)
409 + &y;
410 k_10 = function(t + C_10 * dt, &y_trial)?;
411 y_trial = &k_1 * (A_11_1 * dt)
412 + &k_6 * (A_11_6 * dt)
413 + &k_7 * (A_11_7 * dt)
414 + &k_8 * (A_11_8 * dt)
415 + &k_9 * (A_11_9 * dt)
416 + &k_10 * (A_11_10 * dt)
417 + &y;
418 k_11 = function(t + C_11 * dt, &y_trial)?;
419 y_trial = &k_1 * (A_12_1 * dt)
420 + &k_6 * (A_12_6 * dt)
421 + &k_7 * (A_12_7 * dt)
422 + &k_8 * (A_12_8 * dt)
423 + &k_9 * (A_12_9 * dt)
424 + &k_10 * (A_12_10 * dt)
425 + &k_11 * (A_12_11 * dt)
426 + &y;
427 k_12 = function(t + C_12 * dt, &y_trial)?;
428 y_trial = &k_1 * (A_13_1 * dt)
429 + &k_6 * (A_13_6 * dt)
430 + &k_7 * (A_13_7 * dt)
431 + &k_8 * (A_13_8 * dt)
432 + &k_9 * (A_13_9 * dt)
433 + &k_10 * (A_13_10 * dt)
434 + &k_11 * (A_13_11 * dt)
435 + &k_12 * (A_13_12 * dt)
436 + &y;
437 k_13 = function(t + C_13 * dt, &y_trial)?;
438 y_trial = &k_1 * (A_14_1 * dt)
439 + &k_6 * (A_14_6 * dt)
440 + &k_7 * (A_14_7 * dt)
441 + &k_8 * (A_14_8 * dt)
442 + &k_9 * (A_14_9 * dt)
443 + &k_10 * (A_14_10 * dt)
444 + &k_11 * (A_14_11 * dt)
445 + &k_12 * (A_14_12 * dt)
446 + &k_13 * (A_14_13 * dt)
447 + &y;
448 k_14 = function(t + C_14 * dt, &y_trial)?;
449 y_trial = &k_1 * (A_15_1 * dt)
450 + &k_6 * (A_15_6 * dt)
451 + &k_7 * (A_15_7 * dt)
452 + &k_8 * (A_15_8 * dt)
453 + &k_9 * (A_15_9 * dt)
454 + &k_10 * (A_15_10 * dt)
455 + &k_11 * (A_15_11 * dt)
456 + &k_12 * (A_15_12 * dt)
457 + &k_13 * (A_15_13 * dt)
458 + &k_14 * (A_15_14 * dt)
459 + &y;
460 k_15 = function(t + dt, &y_trial)?;
461 y_trial = (&k_1 * B_1
462 + &k_8 * B_8
463 + &k_9 * B_9
464 + &k_10 * B_10
465 + &k_11 * B_11
466 + &k_12 * B_12
467 + &k_13 * B_13
468 + &k_14 * B_14
469 + &k_15 * B_15)
470 * dt
471 + &y;
472 }
473 dydt_int.push(function(t + dt, &y_trial)?);
474 y_int.push(y_trial);
475 }
476 Ok((y_int, dydt_int))
477 }
478}