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