conspire/math/integrate/ode/explicit/fixed_step/
mod.rs1#[cfg(test)]
2mod test;
3
4use crate::math::{
5 Scalar, Tensor, TensorVec, Vector,
6 integrate::{Explicit, FixedStep, IntegrationError},
7};
8
9pub mod bogacki_shampine;
10pub mod dormand_prince;
11pub mod euler;
12pub mod heun;
13pub mod midpoint;
14pub mod ralston;
15pub mod verner_8;
16pub mod verner_9;
17
18pub trait FixedStepExplicit<Y, U>
20where
21 Self: Explicit<Y, U> + FixedStep,
22 Y: Tensor,
23 U: TensorVec<Item = Y>,
24{
25 fn integrate_fixed_step(
26 &self,
27 mut function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
28 time: &[Scalar],
29 initial_condition: Y,
30 ) -> Result<(Vector, U, U), IntegrationError> {
31 let t_0 = time[0];
32 let t_f = time[time.len() - 1];
33 let mut t_sol: Vector;
34 if time.len() < 2 {
35 return Err(IntegrationError::LengthTimeLessThanTwo);
36 } else if t_0 >= t_f {
37 return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
38 } else if time.len() == 2 {
39 if self.dt() <= 0.0 || self.dt().is_nan() {
40 return Err(IntegrationError::TimeStepNotSet(
41 time[0],
42 time[1],
43 format!("{self:?}"),
44 ));
45 } else {
46 let num_steps = ((t_f - t_0) / self.dt()).ceil() as usize;
47 t_sol = (0..num_steps)
48 .map(|step| t_0 + (step as Scalar) * self.dt())
49 .collect();
50 t_sol.push(t_f);
51 }
52 } else {
53 t_sol = time.iter().copied().collect();
54 }
55 let mut index = 0;
56 let mut t = t_0;
57 let mut dt;
58 let mut t_trial;
59 let mut k = vec![Y::default(); Self::SLOPES];
60 k[0] = function(t, &initial_condition)?;
61 let mut y = initial_condition.clone();
62 let mut y_sol = U::new();
63 y_sol.push(initial_condition.clone());
64 let mut dydt_sol = U::new();
65 dydt_sol.push(function(t, &y.clone())?);
66 let mut y_trial = Y::default();
67 while t < t_f {
68 t_trial = t_sol[index + 1];
69 dt = t_trial - t;
70 if let Err(error) = self.step(&mut function, &y, t, dt, &mut k, &mut y_trial) {
71 return Err(IntegrationError::Upstream(error, format!("{self:?}")));
72 } else {
73 t += dt;
74 y = y_trial.clone();
75 y_sol.push(y.clone());
76 dydt_sol.push(k[0].clone());
77 index += 1;
78 }
79 }
80 Ok((t_sol, y_sol, dydt_sol))
81 }
82 #[allow(clippy::too_many_arguments)]
83 fn step(
84 &self,
85 function: impl FnMut(Scalar, &Y) -> Result<Y, String>,
86 y: &Y,
87 t: Scalar,
88 dt: Scalar,
89 k: &mut [Y],
90 y_trial: &mut Y,
91 ) -> Result<(), String>;
92}