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

1use crate::math::{
2    Derivative, Differentiate, Quantity, Scalar, Tensor, TensorVec,
3    integrate::{ExplicitDaeVariableStepExplicit, ode::explicit::variable_step::verner_9::*},
4};
5use std::ops::{Mul, Sub};
6
7impl<Y, Z, U, V, W, T> ExplicitDaeVariableStepExplicit<Y, Z, U, V, W, T> for Verner9
8where
9    Y: Differentiate<T> + Tensor,
10    Z: PartialEq + Tensor,
11    Derivative<Y, T>: Mul<Quantity<T>, Output = Y>,
12    U: TensorVec<Item = Y>,
13    V: TensorVec<Item = Z>,
14    W: TensorVec<Item = Derivative<Y, T>>,
15    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
16    for<'a> &'a Derivative<Y, T>:
17        Mul<Scalar, Output = Derivative<Y, T>> + Mul<Quantity<T>, Output = Y>,
18{
19    fn slopes_solve(
20        mut evolution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Derivative<Y, T>, String>,
21        mut solution: impl FnMut(Quantity<T>, &Y, &Z) -> Result<Z, String>,
22        y: &Y,
23        z: &Z,
24        t: Quantity<T>,
25        dt: Quantity<T>,
26        k: &mut [Derivative<Y, T>],
27        y_trial: &mut Y,
28        z_trial: &mut Z,
29    ) -> Result<(), String> {
30        k[0] = evolution(t, y, z)?;
31        *y_trial = &k[0] * (A_2_1 * dt) + y;
32        *z_trial = solution(t + C_2 * dt, y_trial, z)?;
33        k[1] = evolution(t + C_2 * dt, y_trial, z_trial)?;
34        *y_trial = &k[0] * (A_3_1 * dt) + &k[1] * (A_3_2 * dt) + y;
35        *z_trial = solution(t + C_3 * dt, y_trial, z_trial)?;
36        k[2] = evolution(t + C_3 * dt, y_trial, z_trial)?;
37        *y_trial = &k[0] * (A_4_1 * dt) + &k[2] * (A_4_3 * dt) + y;
38        *z_trial = solution(t + C_4 * dt, y_trial, z_trial)?;
39        k[3] = evolution(t + C_4 * dt, y_trial, z_trial)?;
40        *y_trial = &k[0] * (A_5_1 * dt) + &k[2] * (A_5_3 * dt) + &k[3] * (A_5_4 * dt) + y;
41        *z_trial = solution(t + C_5 * dt, y_trial, z_trial)?;
42        k[4] = evolution(t + C_5 * dt, y_trial, z_trial)?;
43        *y_trial = &k[0] * (A_6_1 * dt) + &k[3] * (A_6_4 * dt) + &k[4] * (A_6_5 * dt) + y;
44        *z_trial = solution(t + C_6 * dt, y_trial, z_trial)?;
45        k[5] = evolution(t + C_6 * dt, y_trial, z_trial)?;
46        *y_trial = &k[0] * (A_7_1 * dt)
47            + &k[3] * (A_7_4 * dt)
48            + &k[4] * (A_7_5 * dt)
49            + &k[5] * (A_7_6 * dt)
50            + y;
51        *z_trial = solution(t + C_7 * dt, y_trial, z_trial)?;
52        k[6] = evolution(t + C_7 * dt, y_trial, z_trial)?;
53        *y_trial = &k[0] * (A_8_1 * dt) + &k[5] * (A_8_6 * dt) + &k[6] * (A_8_7 * dt) + y;
54        *z_trial = solution(t + C_8 * dt, y_trial, z_trial)?;
55        k[7] = evolution(t + C_8 * dt, y_trial, z_trial)?;
56        *y_trial = &k[0] * (A_9_1 * dt)
57            + &k[5] * (A_9_6 * dt)
58            + &k[6] * (A_9_7 * dt)
59            + &k[7] * (A_9_8 * dt)
60            + y;
61        *z_trial = solution(t + C_9 * dt, y_trial, z_trial)?;
62        k[8] = evolution(t + C_9 * dt, y_trial, z_trial)?;
63        *y_trial = &k[0] * (A_10_1 * dt)
64            + &k[5] * (A_10_6 * dt)
65            + &k[6] * (A_10_7 * dt)
66            + &k[7] * (A_10_8 * dt)
67            + &k[8] * (A_10_9 * dt)
68            + y;
69        *z_trial = solution(t + C_10 * dt, y_trial, z_trial)?;
70        k[9] = evolution(t + C_10 * dt, y_trial, z_trial)?;
71        *y_trial = &k[0] * (A_11_1 * dt)
72            + &k[5] * (A_11_6 * dt)
73            + &k[6] * (A_11_7 * dt)
74            + &k[7] * (A_11_8 * dt)
75            + &k[8] * (A_11_9 * dt)
76            + &k[9] * (A_11_10 * dt)
77            + y;
78        *z_trial = solution(t + C_11 * dt, y_trial, z_trial)?;
79        k[10] = evolution(t + C_11 * dt, y_trial, z_trial)?;
80        *y_trial = &k[0] * (A_12_1 * dt)
81            + &k[5] * (A_12_6 * dt)
82            + &k[6] * (A_12_7 * dt)
83            + &k[7] * (A_12_8 * dt)
84            + &k[8] * (A_12_9 * dt)
85            + &k[9] * (A_12_10 * dt)
86            + &k[10] * (A_12_11 * dt)
87            + y;
88        *z_trial = solution(t + C_12 * dt, y_trial, z_trial)?;
89        k[11] = evolution(t + C_12 * dt, y_trial, z_trial)?;
90        *y_trial = &k[0] * (A_13_1 * dt)
91            + &k[5] * (A_13_6 * dt)
92            + &k[6] * (A_13_7 * dt)
93            + &k[7] * (A_13_8 * dt)
94            + &k[8] * (A_13_9 * dt)
95            + &k[9] * (A_13_10 * dt)
96            + &k[10] * (A_13_11 * dt)
97            + &k[11] * (A_13_12 * dt)
98            + y;
99        *z_trial = solution(t + C_13 * dt, y_trial, z_trial)?;
100        k[12] = evolution(t + C_13 * dt, y_trial, z_trial)?;
101        *y_trial = &k[0] * (A_14_1 * dt)
102            + &k[5] * (A_14_6 * dt)
103            + &k[6] * (A_14_7 * dt)
104            + &k[7] * (A_14_8 * dt)
105            + &k[8] * (A_14_9 * dt)
106            + &k[9] * (A_14_10 * dt)
107            + &k[10] * (A_14_11 * dt)
108            + &k[11] * (A_14_12 * dt)
109            + &k[12] * (A_14_13 * dt)
110            + y;
111        *z_trial = solution(t + C_14 * dt, y_trial, z_trial)?;
112        k[13] = evolution(t + C_14 * dt, y_trial, z_trial)?;
113        *y_trial = &k[0] * (A_15_1 * dt)
114            + &k[5] * (A_15_6 * dt)
115            + &k[6] * (A_15_7 * dt)
116            + &k[7] * (A_15_8 * dt)
117            + &k[8] * (A_15_9 * dt)
118            + &k[9] * (A_15_10 * dt)
119            + &k[10] * (A_15_11 * dt)
120            + &k[11] * (A_15_12 * dt)
121            + &k[12] * (A_15_13 * dt)
122            + &k[13] * (A_15_14 * dt)
123            + y;
124        *z_trial = solution(t + dt, y_trial, z_trial)?;
125        k[14] = evolution(t + dt, y_trial, z_trial)?;
126        *y_trial = &k[0] * (A_16_1 * dt)
127            + &k[5] * (A_16_6 * dt)
128            + &k[6] * (A_16_7 * dt)
129            + &k[7] * (A_16_8 * dt)
130            + &k[8] * (A_16_9 * dt)
131            + &k[9] * (A_16_10 * dt)
132            + &k[10] * (A_16_11 * dt)
133            + &k[11] * (A_16_12 * dt)
134            + &k[12] * (A_16_13 * dt)
135            + y;
136        *z_trial = solution(t + dt, y_trial, z_trial)?;
137        k[15] = evolution(t + dt, y_trial, z_trial)?;
138        *y_trial = (&k[0] * B_1
139            + &k[7] * B_8
140            + &k[8] * B_9
141            + &k[9] * B_10
142            + &k[10] * B_11
143            + &k[11] * B_12
144            + &k[12] * B_13
145            + &k[13] * B_14
146            + &k[14] * B_15)
147            * dt
148            + y;
149        *z_trial = solution(t + dt, y_trial, z_trial)?;
150        Ok(())
151    }
152}
conspire/math/integrate/dae/explicit/variable_step/verner_9/mod.rs

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