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

1use crate::math::{
2    Banded, Scalar, Tensor, TensorVec, Vector,
3    integrate::{
4        ImplicitDaeFirstOrderMinimize, ImplicitDaeFirstOrderRoot, ImplicitDaeSecondOrderMinimize,
5        ImplicitDaeZerothOrderRoot, IntegrationError, VariableStepExplicit,
6    },
7    optimize::{
8        EqualityConstraint, FirstOrderOptimization, FirstOrderRootFinding, SecondOrderOptimization,
9        ZerothOrderRootFinding,
10    },
11};
12use std::ops::{Mul, Sub};
13
14/// Variable-step explicit integrators for implicit differential-algebraic equations.
15pub trait ImplicitDaeVariableStepExplicit<Y, U>
16where
17    Self: VariableStepExplicit<Y, U>,
18    Y: Tensor,
19    U: TensorVec<Item = Y>,
20    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
21{
22    fn integrate_implicit_dae_variable_step(
23        &self,
24        mut evolution: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
25        time: &[Scalar],
26        initial_condition: Y,
27    ) -> Result<(Vector, U, U), IntegrationError> {
28        let t_0 = time[0];
29        let t_f = time[time.len() - 1];
30        if time.len() < 2 {
31            return Err(IntegrationError::LengthTimeLessThanTwo);
32        } else if t_0 >= t_f {
33            return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
34        }
35        let mut t = t_0;
36        let mut dt = t_f - t_0;
37        let mut t_sol = Vector::new();
38        t_sol.push(t_0);
39        let mut dydt = &initial_condition * 0.0;
40        let mut y = initial_condition;
41        let mut k = vec![Y::default(); Self::SLOPES];
42        k[0] = evolution(t, &y, &dydt)?;
43        let mut y_sol = U::new();
44        y_sol.push(y.clone());
45        let mut dydt_sol = U::new();
46        dydt_sol.push(k[0].clone());
47        let mut y_trial = Y::default();
48        while t < t_f {
49            match self.slopes_and_error(
50                |t: Scalar, y: &Y| evolution(t, y, &dydt),
51                &y,
52                t,
53                dt,
54                &mut k,
55                &mut y_trial,
56            ) {
57                Ok(e) => {
58                    if let Some(error) = self
59                        .step(
60                            |t: Scalar, y: &Y| evolution(t, y, &dydt),
61                            &mut y,
62                            &mut t,
63                            &mut y_sol,
64                            &mut t_sol,
65                            &mut dydt_sol,
66                            &mut dt,
67                            &mut k,
68                            &y_trial,
69                            e,
70                        )
71                        .err()
72                    {
73                        dt *= self.dt_cut();
74                        if dt < self.dt_min() {
75                            return Err(IntegrationError::MinimumStepSizeUpstream(
76                                self.dt_min(),
77                                error,
78                                format!("{:?}", self),
79                            ));
80                        }
81                    } else {
82                        dydt = k[0].clone();
83                        dt = dt.min(t_f - t);
84                        if dt < self.dt_min() && t < t_f {
85                            return Err(IntegrationError::MinimumStepSizeReached(
86                                self.dt_min(),
87                                format!("{:?}", self),
88                            ));
89                        }
90                    }
91                }
92                Err(error) => {
93                    dt *= self.dt_cut();
94                    if dt < self.dt_min() {
95                        return Err(IntegrationError::MinimumStepSizeUpstream(
96                            self.dt_min(),
97                            error,
98                            format!("{:?}", self),
99                        ));
100                    }
101                }
102            }
103        }
104        if time.len() > 2 {
105            let t_int = Vector::from(time);
106            let (y_int, dydt_int) = self.interpolate_implicit_dae_variable_step(
107                evolution, &t_int, &t_sol, &y_sol, &dydt_sol,
108            )?;
109            Ok((t_int, y_int, dydt_int))
110        } else {
111            Ok((t_sol, y_sol, dydt_sol))
112        }
113    }
114    fn interpolate_implicit_dae_variable_step(
115        &self,
116        mut evolution: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
117        time: &Vector,
118        tp: &Vector,
119        yp: &U,
120        dydtp: &U,
121    ) -> Result<(U, U), IntegrationError> {
122        let mut dt;
123        let mut i;
124        let mut k = vec![Y::default(); Self::SLOPES];
125        let mut t;
126        let mut y;
127        let mut dydt;
128        let mut y_int = U::new();
129        let mut dydt_int = U::new();
130        let mut y_trial = Y::default();
131        for time_k in time.iter() {
132            i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
133            if time_k == &tp[i] {
134                t = tp[i];
135                y_trial = yp[i].clone();
136                dt = 0.0;
137            } else {
138                t = tp[i - 1];
139                y = &yp[i - 1];
140                dydt = &dydtp[i - 1];
141                dt = time_k - t;
142                k[0] = evolution(t, y, dydt)?;
143                Self::slopes(
144                    |t: Scalar, y: &Y| evolution(t, y, dydt),
145                    y,
146                    t,
147                    dt,
148                    &mut k,
149                    &mut y_trial,
150                )?;
151            }
152            dydt_int.push(evolution(t + dt, &y_trial, &k[0])?);
153            y_int.push(y_trial.clone());
154        }
155        Ok((y_int, dydt_int))
156    }
157}
158
159impl<I, Y, U> ImplicitDaeVariableStepExplicit<Y, U> for I
160where
161    Self: VariableStepExplicit<Y, U>,
162    Y: Tensor,
163    U: TensorVec<Item = Y>,
164    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
165{
166}
167
168/// Variable-step explicit integrators for implicit differential-algebraic equations using zeroth-order root-finding.
169pub trait ImplicitDaeVariableStepExplicitZerothOrderRoot<Y, U>
170where
171    Self: ImplicitDaeVariableStepExplicit<Y, U>,
172    Y: Tensor,
173    U: TensorVec<Item = Y>,
174    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
175{
176    fn integrate_implicit_dae_variable_step_explicit_root_0(
177        &self,
178        mut function: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
179        solver: impl ZerothOrderRootFinding<Y>,
180        time: &[Scalar],
181        initial_condition: Y,
182        mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
183    ) -> Result<(Vector, U, U), IntegrationError> {
184        let evolution = |t: Scalar, y: &Y, dydt_0: &Y| -> Result<Y, String> {
185            Ok(solver.root(
186                |dydt| function(t, y, dydt),
187                dydt_0.clone(),
188                equality_constraint(t),
189            )?)
190        };
191        self.integrate_implicit_dae_variable_step(evolution, time, initial_condition)
192    }
193}
194
195impl<I, Y, U> ImplicitDaeVariableStepExplicitZerothOrderRoot<Y, U> for I
196where
197    I: ImplicitDaeVariableStepExplicit<Y, U>,
198    Y: Tensor,
199    U: TensorVec<Item = Y>,
200    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
201{
202}
203
204impl<I, Y, U> ImplicitDaeZerothOrderRoot<Y, U> for I
205where
206    Self: ImplicitDaeVariableStepExplicitZerothOrderRoot<Y, U>,
207    Y: Tensor,
208    U: TensorVec<Item = Y>,
209    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
210{
211    fn integrate(
212        &self,
213        function: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
214        solver: impl ZerothOrderRootFinding<Y>,
215        time: &[Scalar],
216        initial_condition: Y,
217        equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
218    ) -> Result<(Vector, U, U), IntegrationError> {
219        self.integrate_implicit_dae_variable_step_explicit_root_0(
220            function,
221            solver,
222            time,
223            initial_condition,
224            equality_constraint,
225        )
226    }
227}
228
229/// Variable-step explicit integrators for implicit differential-algebraic equations using first-order root-finding.
230pub trait ImplicitDaeVariableStepExplicitFirstOrderRoot<F, J, Y, U>
231where
232    Self: ImplicitDaeVariableStepExplicit<Y, U>,
233    Y: Tensor,
234    U: TensorVec<Item = Y>,
235    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
236{
237    fn integrate_implicit_dae_variable_step_explicit_root_1(
238        &self,
239        mut function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
240        mut jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<J, String>,
241        solver: impl FirstOrderRootFinding<F, J, Y>,
242        time: &[Scalar],
243        initial_condition: Y,
244        mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
245    ) -> Result<(Vector, U, U), IntegrationError> {
246        let evolution = |t: Scalar, y: &Y, dydt_0: &Y| -> Result<Y, String> {
247            Ok(solver.root(
248                |dydt| function(t, y, dydt),
249                |dydt| jacobian(t, y, dydt),
250                dydt_0.clone(),
251                equality_constraint(t),
252            )?)
253        };
254        self.integrate_implicit_dae_variable_step(evolution, time, initial_condition)
255    }
256}
257
258impl<I, F, J, Y, U> ImplicitDaeVariableStepExplicitFirstOrderRoot<F, J, Y, U> for I
259where
260    I: ImplicitDaeVariableStepExplicit<Y, U>,
261    Y: Tensor,
262    U: TensorVec<Item = Y>,
263    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
264{
265}
266
267impl<I, F, J, Y, U> ImplicitDaeFirstOrderRoot<F, J, Y, U> for I
268where
269    Self: ImplicitDaeVariableStepExplicitFirstOrderRoot<F, J, Y, U>,
270    Y: Tensor,
271    U: TensorVec<Item = Y>,
272    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
273{
274    fn integrate(
275        &self,
276        function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
277        jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<J, String>,
278        solver: impl FirstOrderRootFinding<F, J, Y>,
279        time: &[Scalar],
280        initial_condition: Y,
281        equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
282    ) -> Result<(Vector, U, U), IntegrationError> {
283        self.integrate_implicit_dae_variable_step_explicit_root_1(
284            function,
285            jacobian,
286            solver,
287            time,
288            initial_condition,
289            equality_constraint,
290        )
291    }
292}
293
294/// Variable-step explicit integrators for implicit differential-algebraic equations using first-order minimization.
295pub trait ImplicitDaeVariableStepExplicitFirstOrderMinimize<F, Y, U>
296where
297    Self: ImplicitDaeVariableStepExplicit<Y, U>,
298    Y: Tensor,
299    U: TensorVec<Item = Y>,
300    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
301{
302    #[allow(clippy::too_many_arguments)]
303    fn integrate_implicit_dae_variable_step_explicit_minimize_1(
304        &self,
305        mut function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
306        mut jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
307        solver: impl FirstOrderOptimization<F, Y>,
308        time: &[Scalar],
309        initial_condition: Y,
310        mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
311    ) -> Result<(Vector, U, U), IntegrationError> {
312        let evolution = |t: Scalar, y: &Y, dydt_0: &Y| -> Result<Y, String> {
313            Ok(solver.minimize(
314                |dydt| function(t, y, dydt),
315                |dydt| jacobian(t, y, dydt),
316                dydt_0.clone(),
317                equality_constraint(t),
318            )?)
319        };
320        self.integrate_implicit_dae_variable_step(evolution, time, initial_condition)
321    }
322}
323
324impl<I, F, Y, U> ImplicitDaeVariableStepExplicitFirstOrderMinimize<F, Y, U> for I
325where
326    I: ImplicitDaeVariableStepExplicit<Y, U>,
327    Y: Tensor,
328    U: TensorVec<Item = Y>,
329    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
330{
331}
332
333impl<I, F, Y, U> ImplicitDaeFirstOrderMinimize<F, Y, U> for I
334where
335    Self: ImplicitDaeVariableStepExplicitFirstOrderMinimize<F, Y, U>,
336    Y: Tensor,
337    U: TensorVec<Item = Y>,
338    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
339{
340    fn integrate(
341        &self,
342        function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
343        jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
344        solver: impl FirstOrderOptimization<F, Y>,
345        time: &[Scalar],
346        initial_condition: Y,
347        equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
348    ) -> Result<(Vector, U, U), IntegrationError> {
349        self.integrate_implicit_dae_variable_step_explicit_minimize_1(
350            function,
351            jacobian,
352            solver,
353            time,
354            initial_condition,
355            equality_constraint,
356        )
357    }
358}
359
360/// Variable-step explicit integrators for implicit differential-algebraic equations using second-order minimization.
361pub trait ImplicitDaeVariableStepExplicitSecondOrderMinimize<F, J, H, Y, U>
362where
363    Self: ImplicitDaeVariableStepExplicit<Y, U>,
364    Y: Tensor,
365    U: TensorVec<Item = Y>,
366    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
367{
368    #[allow(clippy::too_many_arguments)]
369    fn integrate_implicit_dae_variable_step_explicit_minimize_2(
370        &self,
371        mut function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
372        mut jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<J, String>,
373        mut hessian: impl FnMut(Scalar, &Y, &Y) -> Result<H, String>,
374        solver: impl SecondOrderOptimization<F, J, H, Y>,
375        time: &[Scalar],
376        initial_condition: Y,
377        mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
378        banded: Option<Banded>,
379    ) -> Result<(Vector, U, U), IntegrationError> {
380        let evolution = |t: Scalar, y: &Y, dydt_0: &Y| -> Result<Y, String> {
381            Ok(solver.minimize(
382                |dydt| function(t, y, dydt),
383                |dydt| jacobian(t, y, dydt),
384                |dydt| hessian(t, y, dydt),
385                dydt_0.clone(),
386                equality_constraint(t),
387                banded.clone(),
388            )?)
389        };
390        self.integrate_implicit_dae_variable_step(evolution, time, initial_condition)
391    }
392}
393
394impl<I, F, J, H, Y, U> ImplicitDaeVariableStepExplicitSecondOrderMinimize<F, J, H, Y, U> for I
395where
396    I: ImplicitDaeVariableStepExplicit<Y, U>,
397    Y: Tensor,
398    U: TensorVec<Item = Y>,
399    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
400{
401}
402
403impl<I, F, J, H, Y, U> ImplicitDaeSecondOrderMinimize<F, J, H, Y, U> for I
404where
405    Self: ImplicitDaeVariableStepExplicitSecondOrderMinimize<F, J, H, Y, U>,
406    Y: Tensor,
407    U: TensorVec<Item = Y>,
408    for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
409{
410    fn integrate(
411        &self,
412        function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
413        jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<J, String>,
414        hessian: impl FnMut(Scalar, &Y, &Y) -> Result<H, String>,
415        solver: impl SecondOrderOptimization<F, J, H, Y>,
416        time: &[Scalar],
417        initial_condition: Y,
418        equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
419        banded: Option<Banded>,
420    ) -> Result<(Vector, U, U), IntegrationError> {
421        self.integrate_implicit_dae_variable_step_explicit_minimize_2(
422            function,
423            jacobian,
424            hessian,
425            solver,
426            time,
427            initial_condition,
428            equality_constraint,
429            banded,
430        )
431    }
432}
conspire/math/integrate/dae/explicit/variable_step/implicit/mod.rs

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