conspire/math/integrate/dae/explicit/variable_step/verner_8/
mod.rs

1use crate::math::{
2    Banded, Scalar, Tensor, TensorVec, Vector,
3    integrate::{
4        ExplicitDaeVariableStep, IntegrationError, ode::explicit::variable_step::verner_8::*,
5    },
6    optimize::{
7        EqualityConstraint, FirstOrderOptimization, FirstOrderRootFinding, SecondOrderOptimization,
8        ZerothOrderRootFinding,
9    },
10};
11use std::ops::{Mul, Sub};
12
13impl<Y, Z, U, V> ExplicitDaeVariableStep<Y, Z, U, V> for Verner8
14where
15    Y: Tensor,
16    Z: Tensor,
17    U: TensorVec<Item = Y>,
18    V: TensorVec<Item = Z>,
19    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
20{
21    fn slopes_solve(
22        mut evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
23        mut solution: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
24        y: &Y,
25        z: &Z,
26        t: Scalar,
27        dt: Scalar,
28        k: &mut [Y],
29        y_trial: &mut Y,
30        z_trial: &mut Z,
31    ) -> Result<(), String> {
32        k[0] = evolution(t, y, z)?;
33        *y_trial = &k[0] * (A_2_1 * dt) + y;
34        *z_trial = solution(t + C_2 * dt, y_trial, z)?;
35        k[1] = evolution(t + C_2 * dt, y_trial, z_trial)?;
36        *y_trial = &k[0] * (A_3_1 * dt) + &k[1] * (A_3_2 * dt) + y;
37        *z_trial = solution(t + C_3 * dt, y_trial, z_trial)?;
38        k[2] = evolution(t + C_3 * dt, y_trial, z_trial)?;
39        *y_trial = &k[0] * (A_4_1 * dt) + &k[2] * (A_4_3 * dt) + y;
40        *z_trial = solution(t + C_4 * dt, y_trial, z_trial)?;
41        k[3] = evolution(t + C_4 * dt, y_trial, z_trial)?;
42        *y_trial = &k[0] * (A_5_1 * dt) + &k[2] * (A_5_3 * dt) + &k[3] * (A_5_4 * dt) + y;
43        *z_trial = solution(t + C_5 * dt, y_trial, z_trial)?;
44        k[4] = evolution(t + C_5 * dt, y_trial, z_trial)?;
45        *y_trial = &k[0] * (A_6_1 * dt) + &k[3] * (A_6_4 * dt) + &k[4] * (A_6_5 * dt) + y;
46        *z_trial = solution(t + C_6 * dt, y_trial, z_trial)?;
47        k[5] = evolution(t + C_6 * dt, y_trial, z_trial)?;
48        *y_trial = &k[0] * (A_7_1 * dt)
49            + &k[3] * (A_7_4 * dt)
50            + &k[4] * (A_7_5 * dt)
51            + &k[5] * (A_7_6 * dt)
52            + y;
53        *z_trial = solution(t + C_7 * dt, y_trial, z_trial)?;
54        k[6] = evolution(t + C_7 * dt, y_trial, z_trial)?;
55        *y_trial = &k[0] * (A_8_1 * dt)
56            + &k[3] * (A_8_4 * dt)
57            + &k[4] * (A_8_5 * dt)
58            + &k[5] * (A_8_6 * dt)
59            + &k[6] * (A_8_7 * dt)
60            + y;
61        *z_trial = solution(t + C_8 * dt, y_trial, z_trial)?;
62        k[7] = evolution(t + C_8 * dt, y_trial, z_trial)?;
63        *y_trial = &k[0] * (A_9_1 * dt)
64            + &k[3] * (A_9_4 * dt)
65            + &k[4] * (A_9_5 * dt)
66            + &k[5] * (A_9_6 * dt)
67            + &k[6] * (A_9_7 * dt)
68            + &k[7] * (A_9_8 * dt)
69            + y;
70        *z_trial = solution(t + C_9 * dt, y_trial, z_trial)?;
71        k[8] = evolution(t + C_9 * dt, y_trial, z_trial)?;
72        *y_trial = &k[0] * (A_10_1 * dt)
73            + &k[3] * (A_10_4 * dt)
74            + &k[4] * (A_10_5 * dt)
75            + &k[5] * (A_10_6 * dt)
76            + &k[6] * (A_10_7 * dt)
77            + &k[7] * (A_10_8 * dt)
78            + &k[8] * (A_10_9 * dt)
79            + y;
80        *z_trial = solution(t + C_10 * dt, y_trial, z_trial)?;
81        k[9] = evolution(t + C_10 * dt, y_trial, z_trial)?;
82        *y_trial = &k[0] * (A_11_1 * dt)
83            + &k[3] * (A_11_4 * dt)
84            + &k[4] * (A_11_5 * dt)
85            + &k[5] * (A_11_6 * dt)
86            + &k[6] * (A_11_7 * dt)
87            + &k[7] * (A_11_8 * dt)
88            + &k[8] * (A_11_9 * dt)
89            + &k[9] * (A_11_10 * dt)
90            + y;
91        *z_trial = solution(t + C_11 * dt, y_trial, z_trial)?;
92        k[10] = evolution(t + C_11 * dt, y_trial, z_trial)?;
93        *y_trial = &k[0] * (A_12_1 * dt)
94            + &k[3] * (A_12_4 * dt)
95            + &k[4] * (A_12_5 * dt)
96            + &k[5] * (A_12_6 * dt)
97            + &k[6] * (A_12_7 * dt)
98            + &k[7] * (A_12_8 * dt)
99            + &k[8] * (A_12_9 * dt)
100            + &k[9] * (A_12_10 * dt)
101            + &k[10] * (A_12_11 * dt)
102            + y;
103        *z_trial = solution(t + dt, y_trial, z_trial)?;
104        k[11] = evolution(t + dt, y_trial, z_trial)?;
105        *y_trial = &k[0] * (A_13_1 * dt)
106            + &k[3] * (A_13_4 * dt)
107            + &k[4] * (A_13_5 * dt)
108            + &k[5] * (A_13_6 * dt)
109            + &k[6] * (A_13_7 * dt)
110            + &k[7] * (A_13_8 * dt)
111            + &k[8] * (A_13_9 * dt)
112            + &k[9] * (A_13_10 * dt)
113            + y;
114        *z_trial = solution(t + dt, y_trial, z_trial)?;
115        k[12] = evolution(t + dt, y_trial, z_trial)?;
116        *y_trial = (&k[0] * B_1
117            + &k[5] * B_6
118            + &k[6] * B_7
119            + &k[7] * B_8
120            + &k[8] * B_9
121            + &k[9] * B_10
122            + &k[10] * B_11
123            + &k[11] * B_12)
124            * dt
125            + y;
126        *z_trial = solution(t + dt, y_trial, z_trial)?;
127        Ok(())
128    }
129}
130
131super::implement_solvers!(Verner8);
conspire/math/integrate/dae/explicit/variable_step/verner_8/mod.rs

conspire/math/integrate/dae/explicit/variable_step/verner_8/
mod.rs
