1use crate::math::{
2 Banded, Scalar, Tensor, TensorVec, Vector, assert_eq_within_tols,
3 integrate::{IntegrationError, VariableStepExplicit},
4 optimize::{
5 EqualityConstraint, FirstOrderOptimization, FirstOrderRootFinding, SecondOrderOptimization,
6 ZerothOrderRootFinding,
7 },
8};
9use std::ops::{Mul, Sub};
10
11pub mod bogacki_shampine;
12pub mod dormand_prince;
13pub mod verner_8;
14pub mod verner_9;
15
16pub trait ExplicitDaeVariableStep<Y, Z, U, V>
17where
18 Self: VariableStepExplicit<Y, U>,
19 Y: Tensor,
20 Z: Tensor,
21 U: TensorVec<Item = Y>,
22 V: TensorVec<Item = Z>,
23 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
24{
25 fn integrate_explicit_dae_variable_step(
26 &self,
27 mut evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
28 mut solution: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
29 time: &[Scalar],
30 initial_condition: (Y, Z),
31 ) -> Result<(Vector, U, U, V), IntegrationError> {
32 let t_0 = time[0];
33 let t_f = time[time.len() - 1];
34 if time.len() < 2 {
35 return Err(IntegrationError::LengthTimeLessThanTwo);
36 } else if t_0 >= t_f {
37 return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
38 }
39 let mut t = t_0;
40 let mut dt = t_f - t_0;
41 let mut t_sol = Vector::new();
42 t_sol.push(t_0);
43 let (mut y, mut z) = initial_condition;
44 if assert_eq_within_tols(&solution(t, &y, &z)?, &z).is_err() {
45 return Err(IntegrationError::InconsistentInitialConditions);
46 }
47 let mut k = vec![Y::default(); Self::SLOPES];
48 k[0] = evolution(t, &y, &z)?;
49 let mut y_sol = U::new();
50 y_sol.push(y.clone());
51 let mut z_sol = V::new();
52 z_sol.push(z.clone());
53 let mut dydt_sol = U::new();
54 dydt_sol.push(k[0].clone());
55 let mut y_trial = Y::default();
56 let mut z_trial = Z::default();
57 while t < t_f {
58 match self.slopes_solve_and_error(
59 &mut evolution,
60 &mut solution,
61 &y,
62 &z,
63 t,
64 dt,
65 &mut k,
66 &mut y_trial,
67 &mut z_trial,
68 ) {
69 Ok(e) => {
70 if let Some(error) = self
71 .step_solve(
72 &mut evolution,
73 &mut y,
74 &mut z,
75 &mut t,
76 &mut y_sol,
77 &mut z_sol,
78 &mut t_sol,
79 &mut dydt_sol,
80 &mut dt,
81 &mut k,
82 &y_trial,
83 &z_trial,
84 e,
85 )
86 .err()
87 {
88 dt *= self.dt_cut();
89 if dt < self.dt_min() {
90 return Err(IntegrationError::MinimumStepSizeUpstream(
91 self.dt_min(),
92 error,
93 format!("{:?}", self),
94 ));
95 }
96 } else {
97 dt = dt.min(t_f - t);
98 if dt < self.dt_min() && t < t_f {
99 return Err(IntegrationError::MinimumStepSizeReached(
100 self.dt_min(),
101 format!("{:?}", self),
102 ));
103 }
104 }
105 }
106 Err(error) => {
107 dt *= self.dt_cut();
108 if dt < self.dt_min() {
109 return Err(IntegrationError::MinimumStepSizeUpstream(
110 self.dt_min(),
111 error,
112 format!("{:?}", self),
113 ));
114 }
115 }
116 }
117 }
118 if time.len() > 2 {
119 let t_int = Vector::from(time);
120 let (y_int, dydt_int, z_int) = self.interpolate_explicit_dae_variable_step(
121 evolution, solution, &t_int, &t_sol, &y_sol, &z_sol,
122 )?;
123 Ok((t_int, y_int, dydt_int, z_int))
124 } else {
125 Ok((t_sol, y_sol, dydt_sol, z_sol))
126 }
127 }
128 fn interpolate_explicit_dae_variable_step(
129 &self,
130 mut evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
131 mut solution: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
132 time: &Vector,
133 tp: &Vector,
134 yp: &U,
135 zp: &V,
136 ) -> Result<(U, U, V), IntegrationError> {
137 let mut dt;
138 let mut i;
139 let mut k = vec![Y::default(); Self::SLOPES];
140 let mut t;
141 let mut y;
142 let mut z;
143 let mut y_int = U::new();
144 let mut z_int = V::new();
145 let mut dydt_int = U::new();
146 let mut y_trial = Y::default();
147 let mut z_trial = Z::default();
148 for time_k in time.iter() {
149 i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
150 if time_k == &tp[i] {
151 t = tp[i];
152 y_trial = yp[i].clone();
153 z_trial = zp[i].clone();
154 dt = 0.0;
155 } else {
156 t = tp[i - 1];
157 y = &yp[i - 1];
158 z = &zp[i - 1];
159 dt = time_k - t;
160 k[0] = evolution(t, y, z)?;
161 Self::slopes_solve(
162 &mut evolution,
163 &mut solution,
164 y,
165 z,
166 t,
167 dt,
168 &mut k,
169 &mut y_trial,
170 &mut z_trial,
171 )?;
172 }
173 dydt_int.push(evolution(t + dt, &y_trial, &z_trial)?);
174 y_int.push(y_trial.clone());
175 z_int.push(z_trial.clone());
176 }
177 Ok((y_int, dydt_int, z_int))
178 }
179 #[allow(clippy::too_many_arguments)]
180 fn slopes_solve(
181 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
182 solution: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
183 y: &Y,
184 z: &Z,
185 t: Scalar,
186 dt: Scalar,
187 k: &mut [Y],
188 y_trial: &mut Y,
189 z_trial: &mut Z,
190 ) -> Result<(), String>;
191 #[allow(clippy::too_many_arguments)]
192 fn slopes_solve_and_error(
193 &self,
194 mut evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
195 mut solution: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
196 y: &Y,
197 z: &Z,
198 t: Scalar,
199 dt: Scalar,
200 k: &mut [Y],
201 y_trial: &mut Y,
202 z_trial: &mut Z,
203 ) -> Result<Scalar, String> {
204 Self::slopes_solve(
205 &mut evolution,
206 &mut solution,
207 y,
208 z,
209 t,
210 dt,
211 k,
212 y_trial,
213 z_trial,
214 )?;
215 Self::error(dt, k)
216 }
217 #[allow(clippy::too_many_arguments)]
218 fn step_solve(
219 &self,
220 mut evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
221 y: &mut Y,
222 z: &mut Z,
223 t: &mut Scalar,
224 y_sol: &mut U,
225 z_sol: &mut V,
226 t_sol: &mut Vector,
227 dydt_sol: &mut U,
228 dt: &mut Scalar,
229 _k: &mut [Y],
230 y_trial: &Y,
231 z_trial: &Z,
232 e: Scalar,
233 ) -> Result<(), String> {
234 if e < self.abs_tol() || e / y_trial.norm_inf() < self.rel_tol() {
235 *t += *dt;
236 *y = y_trial.clone();
237 *z = z_trial.clone();
238 t_sol.push(*t);
239 y_sol.push(y.clone());
240 z_sol.push(z.clone());
241 dydt_sol.push(evolution(*t, y, z)?);
242 }
243 self.time_step(e, dt);
244 Ok(())
245 }
246}
247
248pub trait ImplicitDaeVariableStep<Y, U>
249where
250 Self: VariableStepExplicit<Y, U>,
251 Y: Tensor,
252 U: TensorVec<Item = Y>,
253 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
254{
255 fn integrate_implicit_dae_variable_step(
256 &self,
257 mut evolution: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
258 time: &[Scalar],
259 initial_condition: Y,
260 ) -> Result<(Vector, U, U), IntegrationError> {
261 let t_0 = time[0];
262 let t_f = time[time.len() - 1];
263 if time.len() < 2 {
264 return Err(IntegrationError::LengthTimeLessThanTwo);
265 } else if t_0 >= t_f {
266 return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
267 }
268 let mut t = t_0;
269 let mut dt = t_f - t_0;
270 let mut t_sol = Vector::new();
271 t_sol.push(t_0);
272 let mut dydt = &initial_condition * 0.0;
273 let mut y = initial_condition;
274 let mut k = vec![Y::default(); Self::SLOPES];
275 k[0] = evolution(t, &y, &dydt)?;
276 let mut y_sol = U::new();
277 y_sol.push(y.clone());
278 let mut dydt_sol = U::new();
279 dydt_sol.push(k[0].clone());
280 let mut y_trial = Y::default();
281 while t < t_f {
282 match self.slopes_and_error(
283 |t: Scalar, y: &Y| evolution(t, y, &dydt),
284 &y,
285 t,
286 dt,
287 &mut k,
288 &mut y_trial,
289 ) {
290 Ok(e) => {
291 if let Some(error) = self
292 .step(
293 |t: Scalar, y: &Y| evolution(t, y, &dydt),
294 &mut y,
295 &mut t,
296 &mut y_sol,
297 &mut t_sol,
298 &mut dydt_sol,
299 &mut dt,
300 &mut k,
301 &y_trial,
302 e,
303 )
304 .err()
305 {
306 dt *= self.dt_cut();
307 if dt < self.dt_min() {
308 return Err(IntegrationError::MinimumStepSizeUpstream(
309 self.dt_min(),
310 error,
311 format!("{:?}", self),
312 ));
313 }
314 } else {
315 dydt = k[0].clone();
316 dt = dt.min(t_f - t);
317 if dt < self.dt_min() && t < t_f {
318 return Err(IntegrationError::MinimumStepSizeReached(
319 self.dt_min(),
320 format!("{:?}", self),
321 ));
322 }
323 }
324 }
325 Err(error) => {
326 dt *= self.dt_cut();
327 if dt < self.dt_min() {
328 return Err(IntegrationError::MinimumStepSizeUpstream(
329 self.dt_min(),
330 error,
331 format!("{:?}", self),
332 ));
333 }
334 }
335 }
336 }
337 if time.len() > 2 {
338 let t_int = Vector::from(time);
339 let (y_int, dydt_int) = self.interpolate_implicit_dae_variable_step(
340 evolution, &t_int, &t_sol, &y_sol, &dydt_sol,
341 )?;
342 Ok((t_int, y_int, dydt_int))
343 } else {
344 Ok((t_sol, y_sol, dydt_sol))
345 }
346 }
347 fn interpolate_implicit_dae_variable_step(
348 &self,
349 mut evolution: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
350 time: &Vector,
351 tp: &Vector,
352 yp: &U,
353 dydtp: &U,
354 ) -> Result<(U, U), IntegrationError> {
355 let mut dt;
356 let mut i;
357 let mut k = vec![Y::default(); Self::SLOPES];
358 let mut t;
359 let mut y;
360 let mut dydt;
361 let mut y_int = U::new();
362 let mut dydt_int = U::new();
363 let mut y_trial = Y::default();
364 for time_k in time.iter() {
365 i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
366 if time_k == &tp[i] {
367 t = tp[i];
368 y_trial = yp[i].clone();
369 dt = 0.0;
370 } else {
371 t = tp[i - 1];
372 y = &yp[i - 1];
373 dydt = &dydtp[i - 1];
374 dt = time_k - t;
375 k[0] = evolution(t, y, dydt)?;
376 Self::slopes(
377 |t: Scalar, y: &Y| evolution(t, y, dydt),
378 y,
379 t,
380 dt,
381 &mut k,
382 &mut y_trial,
383 )?;
384 }
385 dydt_int.push(evolution(t + dt, &y_trial, &k[0])?);
386 y_int.push(y_trial.clone());
387 }
388 Ok((y_int, dydt_int))
389 }
390}
391
392pub trait ExplicitDaeVariableStepFirstSameAsLast<Y, Z, U, V>
393where
394 Self: ExplicitDaeVariableStep<Y, Z, U, V>,
395 Y: Tensor,
396 Z: Tensor,
397 U: TensorVec<Item = Y>,
398 V: TensorVec<Item = Z>,
399 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
400{
401 #[allow(clippy::too_many_arguments)]
402 fn slopes_solve_and_error_fsal(
403 mut evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
404 mut solution: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
405 y: &Y,
406 z: &Z,
407 t: Scalar,
408 dt: Scalar,
409 k: &mut [Y],
410 y_trial: &mut Y,
411 z_trial: &mut Z,
412 ) -> Result<Scalar, String> {
413 Self::slopes_solve(
414 &mut evolution,
415 &mut solution,
416 y,
417 z,
418 t,
419 dt,
420 k,
421 y_trial,
422 z_trial,
423 )?;
424 k[Self::SLOPES - 1] = evolution(t + dt, y_trial, z_trial)?;
425 Self::error(dt, k)
426 }
427 #[allow(clippy::too_many_arguments)]
428 fn step_solve_fsal(
429 &self,
430 y: &mut Y,
431 z: &mut Z,
432 t: &mut Scalar,
433 y_sol: &mut U,
434 z_sol: &mut V,
435 t_sol: &mut Vector,
436 dydt_sol: &mut U,
437 dt: &mut Scalar,
438 k: &mut [Y],
439 y_trial: &Y,
440 z_trial: &Z,
441 e: Scalar,
442 ) -> Result<(), String> {
443 if e < self.abs_tol() || e / y_trial.norm_inf() < self.rel_tol() {
444 k[0] = k[Self::SLOPES - 1].clone();
445 *t += *dt;
446 *y = y_trial.clone();
447 *z = z_trial.clone();
448 t_sol.push(*t);
449 y_sol.push(y.clone());
450 z_sol.push(z.clone());
451 dydt_sol.push(k[0].clone());
452 }
453 self.time_step(e, dt);
454 Ok(())
455 }
456}
457
458pub trait ExplicitDaeVariableStepExplicitZerothOrderRoot<Y, Z, U, V>
459where
460 Self: ExplicitDaeVariableStep<Y, Z, U, V>,
461 Y: Tensor,
462 Z: Tensor,
463 U: TensorVec<Item = Y>,
464 V: TensorVec<Item = Z>,
465 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
466{
467 fn integrate_explicit_dae_variable_step_explicit_root_0(
468 &self,
469 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
470 mut function: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
471 solver: impl ZerothOrderRootFinding<Z>,
472 time: &[Scalar],
473 initial_condition: (Y, Z),
474 mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
475 ) -> Result<(Vector, U, U, V), IntegrationError> {
476 let solution = |t: Scalar, y: &Y, z_0: &Z| -> Result<Z, String> {
477 Ok(solver.root(|z| function(t, y, z), z_0.clone(), equality_constraint(t))?)
478 };
479 self.integrate_explicit_dae_variable_step(evolution, solution, time, initial_condition)
480 }
481}
482
483pub trait ExplicitDaeVariableStepExplicitFirstOrderRoot<F, J, Y, Z, U, V>
484where
485 Self: ExplicitDaeVariableStep<Y, Z, U, V>,
486 Y: Tensor,
487 Z: Tensor,
488 U: TensorVec<Item = Y>,
489 V: TensorVec<Item = Z>,
490 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
491{
492 #[allow(clippy::too_many_arguments)]
493 fn integrate_explicit_dae_variable_step_explicit_root_1(
494 &self,
495 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
496 mut function: impl FnMut(Scalar, &Y, &Z) -> Result<F, String>,
497 mut jacobian: impl FnMut(Scalar, &Y, &Z) -> Result<J, String>,
498 solver: impl FirstOrderRootFinding<F, J, Z>,
499 time: &[Scalar],
500 initial_condition: (Y, Z),
501 mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
502 ) -> Result<(Vector, U, U, V), IntegrationError> {
503 let solution = |t: Scalar, y: &Y, z_0: &Z| -> Result<Z, String> {
504 Ok(solver.root(
505 |z| function(t, y, z),
506 |z| jacobian(t, y, z),
507 z_0.clone(),
508 equality_constraint(t),
509 )?)
510 };
511 self.integrate_explicit_dae_variable_step(evolution, solution, time, initial_condition)
512 }
513}
514
515pub trait ExplicitDaeVariableStepExplicitFirstOrderMinimize<F, Y, Z, U, V>
516where
517 Self: ExplicitDaeVariableStep<Y, Z, U, V>,
518 Y: Tensor,
519 Z: Tensor,
520 U: TensorVec<Item = Y>,
521 V: TensorVec<Item = Z>,
522 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
523{
524 #[allow(clippy::too_many_arguments)]
525 fn integrate_explicit_dae_variable_step_explicit_minimize_1(
526 &self,
527 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
528 mut function: impl FnMut(Scalar, &Y, &Z) -> Result<F, String>,
529 mut jacobian: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
530 solver: impl FirstOrderOptimization<F, Z>,
531 time: &[Scalar],
532 initial_condition: (Y, Z),
533 mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
534 ) -> Result<(Vector, U, U, V), IntegrationError> {
535 let solution = |t: Scalar, y: &Y, z_0: &Z| -> Result<Z, String> {
536 Ok(solver.minimize(
537 |z| function(t, y, z),
538 |z| jacobian(t, y, z),
539 z_0.clone(),
540 equality_constraint(t),
541 )?)
542 };
543 self.integrate_explicit_dae_variable_step(evolution, solution, time, initial_condition)
544 }
545}
546
547pub trait ExplicitDaeVariableStepExplicitSecondOrderMinimize<F, J, H, Y, Z, U, V>
548where
549 Self: ExplicitDaeVariableStep<Y, Z, U, V>,
550 Y: Tensor,
551 Z: Tensor,
552 U: TensorVec<Item = Y>,
553 V: TensorVec<Item = Z>,
554 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
555{
556 #[allow(clippy::too_many_arguments)]
557 fn integrate_explicit_dae_variable_step_explicit_minimize_2(
558 &self,
559 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
560 mut function: impl FnMut(Scalar, &Y, &Z) -> Result<F, String>,
561 mut jacobian: impl FnMut(Scalar, &Y, &Z) -> Result<J, String>,
562 mut hessian: impl FnMut(Scalar, &Y, &Z) -> Result<H, String>,
563 solver: impl SecondOrderOptimization<F, J, H, Z>,
564 time: &[Scalar],
565 initial_condition: (Y, Z),
566 mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
567 banded: Option<Banded>,
568 ) -> Result<(Vector, U, U, V), IntegrationError> {
569 let solution = |t: Scalar, y: &Y, z_0: &Z| -> Result<Z, String> {
570 Ok(solver.minimize(
571 |z| function(t, y, z),
572 |z| jacobian(t, y, z),
573 |z| hessian(t, y, z),
574 z_0.clone(),
575 equality_constraint(t),
576 banded.clone(),
577 )?)
578 };
579 self.integrate_explicit_dae_variable_step(evolution, solution, time, initial_condition)
580 }
581}
582
583pub trait ImplicitDaeVariableStepExplicitZerothOrderRoot<Y, U>
584where
585 Self: ImplicitDaeVariableStep<Y, U>,
586 Y: Tensor,
587 U: TensorVec<Item = Y>,
588 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
589{
590 fn integrate_implicit_dae_variable_step_explicit_root_0(
591 &self,
592 mut function: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
593 solver: impl ZerothOrderRootFinding<Y>,
594 time: &[Scalar],
595 initial_condition: Y,
596 mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
597 ) -> Result<(Vector, U, U), IntegrationError> {
598 let evolution = |t: Scalar, y: &Y, dydt_0: &Y| -> Result<Y, String> {
599 Ok(solver.root(
600 |dydt| function(t, y, dydt),
601 dydt_0.clone(),
602 equality_constraint(t),
603 )?)
604 };
605 self.integrate_implicit_dae_variable_step(evolution, time, initial_condition)
606 }
607}
608
609pub trait ImplicitDaeVariableStepExplicitFirstOrderRoot<F, J, Y, U>
610where
611 Self: ImplicitDaeVariableStep<Y, U>,
612 Y: Tensor,
613 U: TensorVec<Item = Y>,
614 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
615{
616 fn integrate_implicit_dae_variable_step_explicit_root_1(
617 &self,
618 mut function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
619 mut jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<J, String>,
620 solver: impl FirstOrderRootFinding<F, J, Y>,
621 time: &[Scalar],
622 initial_condition: Y,
623 mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
624 ) -> Result<(Vector, U, U), IntegrationError> {
625 let evolution = |t: Scalar, y: &Y, dydt_0: &Y| -> Result<Y, String> {
626 Ok(solver.root(
627 |dydt| function(t, y, dydt),
628 |dydt| jacobian(t, y, dydt),
629 dydt_0.clone(),
630 equality_constraint(t),
631 )?)
632 };
633 self.integrate_implicit_dae_variable_step(evolution, time, initial_condition)
634 }
635}
636
637pub trait ImplicitDaeVariableStepExplicitFirstOrderMinimize<F, Y, U>
638where
639 Self: ImplicitDaeVariableStep<Y, U>,
640 Y: Tensor,
641 U: TensorVec<Item = Y>,
642 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
643{
644 #[allow(clippy::too_many_arguments)]
645 fn integrate_implicit_dae_variable_step_explicit_minimize_1(
646 &self,
647 mut function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
648 mut jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
649 solver: impl FirstOrderOptimization<F, Y>,
650 time: &[Scalar],
651 initial_condition: Y,
652 mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
653 ) -> Result<(Vector, U, U), IntegrationError> {
654 let evolution = |t: Scalar, y: &Y, dydt_0: &Y| -> Result<Y, String> {
655 Ok(solver.minimize(
656 |dydt| function(t, y, dydt),
657 |dydt| jacobian(t, y, dydt),
658 dydt_0.clone(),
659 equality_constraint(t),
660 )?)
661 };
662 self.integrate_implicit_dae_variable_step(evolution, time, initial_condition)
663 }
664}
665
666pub trait ImplicitDaeVariableStepExplicitSecondOrderMinimize<F, J, H, Y, U>
667where
668 Self: ImplicitDaeVariableStep<Y, U>,
669 Y: Tensor,
670 U: TensorVec<Item = Y>,
671 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
672{
673 #[allow(clippy::too_many_arguments)]
674 fn integrate_implicit_dae_variable_step_explicit_minimize_2(
675 &self,
676 mut function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
677 mut jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<J, String>,
678 mut hessian: impl FnMut(Scalar, &Y, &Y) -> Result<H, String>,
679 solver: impl SecondOrderOptimization<F, J, H, Y>,
680 time: &[Scalar],
681 initial_condition: Y,
682 mut equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
683 banded: Option<Banded>,
684 ) -> Result<(Vector, U, U), IntegrationError> {
685 let evolution = |t: Scalar, y: &Y, dydt_0: &Y| -> Result<Y, String> {
686 Ok(solver.minimize(
687 |dydt| function(t, y, dydt),
688 |dydt| jacobian(t, y, dydt),
689 |dydt| hessian(t, y, dydt),
690 dydt_0.clone(),
691 equality_constraint(t),
692 banded.clone(),
693 )?)
694 };
695 self.integrate_implicit_dae_variable_step(evolution, time, initial_condition)
696 }
697}
698
699macro_rules! implement_solvers {
700 ($integrator:ident) => {
701 use crate::math::integrate::{
702 ExplicitDaeFirstOrderMinimize, ExplicitDaeFirstOrderRoot,
703 ExplicitDaeSecondOrderMinimize, ExplicitDaeVariableStepExplicitFirstOrderMinimize,
704 ExplicitDaeVariableStepExplicitFirstOrderRoot,
705 ExplicitDaeVariableStepExplicitSecondOrderMinimize,
706 ExplicitDaeVariableStepExplicitZerothOrderRoot, ExplicitDaeZerothOrderRoot,
707 ImplicitDaeFirstOrderMinimize, ImplicitDaeFirstOrderRoot,
708 ImplicitDaeSecondOrderMinimize, ImplicitDaeVariableStep,
709 ImplicitDaeVariableStepExplicitFirstOrderMinimize,
710 ImplicitDaeVariableStepExplicitFirstOrderRoot,
711 ImplicitDaeVariableStepExplicitSecondOrderMinimize,
712 ImplicitDaeVariableStepExplicitZerothOrderRoot, ImplicitDaeZerothOrderRoot,
713 };
714 impl<Y, Z, U, V> ExplicitDaeZerothOrderRoot<Y, Z, U, V> for $integrator
715 where
716 Y: Tensor,
717 Z: Tensor,
718 U: TensorVec<Item = Y>,
719 V: TensorVec<Item = Z>,
720 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
721 {
722 fn integrate(
723 &self,
724 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
725 function: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
726 solver: impl ZerothOrderRootFinding<Z>,
727 time: &[Scalar],
728 initial_condition: (Y, Z),
729 equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
730 ) -> Result<(Vector, U, U, V), IntegrationError> {
731 self.integrate_explicit_dae_variable_step_explicit_root_0(
732 evolution,
733 function,
734 solver,
735 time,
736 initial_condition,
737 equality_constraint,
738 )
739 }
740 }
741 impl<Y, Z, U, V> ExplicitDaeVariableStepExplicitZerothOrderRoot<Y, Z, U, V> for $integrator
742 where
743 Y: Tensor,
744 Z: Tensor,
745 U: TensorVec<Item = Y>,
746 V: TensorVec<Item = Z>,
747 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
748 {
749 }
750 impl<F, J, Y, Z, U, V> ExplicitDaeFirstOrderRoot<F, J, Y, Z, U, V> for $integrator
751 where
752 Y: Tensor,
753 Z: Tensor,
754 U: TensorVec<Item = Y>,
755 V: TensorVec<Item = Z>,
756 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
757 {
758 fn integrate(
759 &self,
760 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
761 function: impl FnMut(Scalar, &Y, &Z) -> Result<F, String>,
762 jacobian: impl FnMut(Scalar, &Y, &Z) -> Result<J, String>,
763 solver: impl FirstOrderRootFinding<F, J, Z>,
764 time: &[Scalar],
765 initial_condition: (Y, Z),
766 equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
767 ) -> Result<(Vector, U, U, V), IntegrationError> {
768 self.integrate_explicit_dae_variable_step_explicit_root_1(
769 evolution,
770 function,
771 jacobian,
772 solver,
773 time,
774 initial_condition,
775 equality_constraint,
776 )
777 }
778 }
779 impl<F, J, Y, Z, U, V> ExplicitDaeVariableStepExplicitFirstOrderRoot<F, J, Y, Z, U, V>
780 for $integrator
781 where
782 Y: Tensor,
783 Z: Tensor,
784 U: TensorVec<Item = Y>,
785 V: TensorVec<Item = Z>,
786 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
787 {
788 }
789 impl<F, Y, Z, U, V> ExplicitDaeFirstOrderMinimize<F, Y, Z, U, V> for $integrator
790 where
791 Y: Tensor,
792 Z: Tensor,
793 U: TensorVec<Item = Y>,
794 V: TensorVec<Item = Z>,
795 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
796 {
797 fn integrate(
798 &self,
799 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
800 function: impl FnMut(Scalar, &Y, &Z) -> Result<F, String>,
801 jacobian: impl FnMut(Scalar, &Y, &Z) -> Result<Z, String>,
802 solver: impl FirstOrderOptimization<F, Z>,
803 time: &[Scalar],
804 initial_condition: (Y, Z),
805 equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
806 ) -> Result<(Vector, U, U, V), IntegrationError> {
807 self.integrate_explicit_dae_variable_step_explicit_minimize_1(
808 evolution,
809 function,
810 jacobian,
811 solver,
812 time,
813 initial_condition,
814 equality_constraint,
815 )
816 }
817 }
818 impl<F, Y, Z, U, V> ExplicitDaeVariableStepExplicitFirstOrderMinimize<F, Y, Z, U, V>
819 for $integrator
820 where
821 Y: Tensor,
822 Z: Tensor,
823 U: TensorVec<Item = Y>,
824 V: TensorVec<Item = Z>,
825 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
826 {
827 }
828 impl<F, J, H, Y, Z, U, V> ExplicitDaeSecondOrderMinimize<F, J, H, Y, Z, U, V>
829 for $integrator
830 where
831 Y: Tensor,
832 Z: Tensor,
833 U: TensorVec<Item = Y>,
834 V: TensorVec<Item = Z>,
835 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
836 {
837 fn integrate(
838 &self,
839 evolution: impl FnMut(Scalar, &Y, &Z) -> Result<Y, String>,
840 function: impl FnMut(Scalar, &Y, &Z) -> Result<F, String>,
841 jacobian: impl FnMut(Scalar, &Y, &Z) -> Result<J, String>,
842 hessian: impl FnMut(Scalar, &Y, &Z) -> Result<H, String>,
843 solver: impl SecondOrderOptimization<F, J, H, Z>,
844 time: &[Scalar],
845 initial_condition: (Y, Z),
846 equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
847 banded: Option<Banded>,
848 ) -> Result<(Vector, U, U, V), IntegrationError> {
849 self.integrate_explicit_dae_variable_step_explicit_minimize_2(
850 evolution,
851 function,
852 jacobian,
853 hessian,
854 solver,
855 time,
856 initial_condition,
857 equality_constraint,
858 banded,
859 )
860 }
861 }
862 impl<F, J, H, Y, Z, U, V>
863 ExplicitDaeVariableStepExplicitSecondOrderMinimize<F, J, H, Y, Z, U, V> for $integrator
864 where
865 Y: Tensor,
866 Z: Tensor,
867 U: TensorVec<Item = Y>,
868 V: TensorVec<Item = Z>,
869 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
870 {
871 }
872 impl<Y, U> ImplicitDaeVariableStep<Y, U> for $integrator
873 where
874 Y: Tensor,
875 U: TensorVec<Item = Y>,
876 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
877 {
878 }
879 impl<Y, U> ImplicitDaeZerothOrderRoot<Y, U> for $integrator
880 where
881 Y: Tensor,
882 U: TensorVec<Item = Y>,
883 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
884 {
885 fn integrate(
886 &self,
887 function: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
888 solver: impl ZerothOrderRootFinding<Y>,
889 time: &[Scalar],
890 initial_condition: Y,
891 equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
892 ) -> Result<(Vector, U, U), IntegrationError> {
893 self.integrate_implicit_dae_variable_step_explicit_root_0(
894 function,
895 solver,
896 time,
897 initial_condition,
898 equality_constraint,
899 )
900 }
901 }
902 impl<Y, U> ImplicitDaeVariableStepExplicitZerothOrderRoot<Y, U> for $integrator
903 where
904 Y: Tensor,
905 U: TensorVec<Item = Y>,
906 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
907 {
908 }
909 impl<F, J, Y, U> ImplicitDaeFirstOrderRoot<F, J, Y, U> for $integrator
910 where
911 Y: Tensor,
912 U: TensorVec<Item = Y>,
913 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
914 {
915 fn integrate(
916 &self,
917 function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
918 jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<J, String>,
919 solver: impl FirstOrderRootFinding<F, J, Y>,
920 time: &[Scalar],
921 initial_condition: Y,
922 equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
923 ) -> Result<(Vector, U, U), IntegrationError> {
924 self.integrate_implicit_dae_variable_step_explicit_root_1(
925 function,
926 jacobian,
927 solver,
928 time,
929 initial_condition,
930 equality_constraint,
931 )
932 }
933 }
934 impl<F, J, Y, U> ImplicitDaeVariableStepExplicitFirstOrderRoot<F, J, Y, U> for $integrator
935 where
936 Y: Tensor,
937 U: TensorVec<Item = Y>,
938 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
939 {
940 }
941 impl<F, Y, U> ImplicitDaeFirstOrderMinimize<F, Y, U> for $integrator
942 where
943 Y: Tensor,
944 U: TensorVec<Item = Y>,
945 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
946 {
947 fn integrate(
948 &self,
949 function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
950 jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<Y, String>,
951 solver: impl FirstOrderOptimization<F, Y>,
952 time: &[Scalar],
953 initial_condition: Y,
954 equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
955 ) -> Result<(Vector, U, U), IntegrationError> {
956 self.integrate_implicit_dae_variable_step_explicit_minimize_1(
957 function,
958 jacobian,
959 solver,
960 time,
961 initial_condition,
962 equality_constraint,
963 )
964 }
965 }
966 impl<F, Y, U> ImplicitDaeVariableStepExplicitFirstOrderMinimize<F, Y, U> for $integrator
967 where
968 Y: Tensor,
969 U: TensorVec<Item = Y>,
970 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
971 {
972 }
973 impl<F, J, H, Y, U> ImplicitDaeSecondOrderMinimize<F, J, H, Y, U> for $integrator
974 where
975 Y: Tensor,
976 U: TensorVec<Item = Y>,
977 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
978 {
979 fn integrate(
980 &self,
981 function: impl FnMut(Scalar, &Y, &Y) -> Result<F, String>,
982 jacobian: impl FnMut(Scalar, &Y, &Y) -> Result<J, String>,
983 hessian: impl FnMut(Scalar, &Y, &Y) -> Result<H, String>,
984 solver: impl SecondOrderOptimization<F, J, H, Y>,
985 time: &[Scalar],
986 initial_condition: Y,
987 equality_constraint: impl FnMut(Scalar) -> EqualityConstraint,
988 banded: Option<Banded>,
989 ) -> Result<(Vector, U, U), IntegrationError> {
990 self.integrate_implicit_dae_variable_step_explicit_minimize_2(
991 function,
992 jacobian,
993 hessian,
994 solver,
995 time,
996 initial_condition,
997 equality_constraint,
998 banded,
999 )
1000 }
1001 }
1002 impl<F, J, H, Y, U> ImplicitDaeVariableStepExplicitSecondOrderMinimize<F, J, H, Y, U>
1003 for $integrator
1004 where
1005 Y: Tensor,
1006 U: TensorVec<Item = Y>,
1007 for<'a> &'a Y: Mul<Scalar, Output = Y> + Sub<&'a Y, Output = Y>,
1008 {
1009 }
1010 };
1011}
1012pub(crate) use implement_solvers;