conspire/math/integrate/bogacki_shampine/
mod.rs1#[cfg(test)]
2mod test;
3
4use super::{
5 super::{Tensor, TensorRank0, TensorVec, Vector, interpolate::InterpolateSolution},
6 Explicit, IntegrationError,
7};
8use crate::{ABS_TOL, REL_TOL};
9use std::ops::{Mul, Sub};
10
11#[derive(Debug)]
43pub struct BogackiShampine {
44 pub abs_tol: TensorRank0,
46 pub rel_tol: TensorRank0,
48 pub dt_beta: TensorRank0,
50 pub dt_expn: TensorRank0,
52}
53
54impl Default for BogackiShampine {
55 fn default() -> Self {
56 Self {
57 abs_tol: ABS_TOL,
58 rel_tol: REL_TOL,
59 dt_beta: 0.9,
60 dt_expn: 3.0,
61 }
62 }
63}
64
65impl<Y, U> Explicit<Y, U> for BogackiShampine
66where
67 Self: InterpolateSolution<Y, U>,
68 Y: Tensor,
69 for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
70 U: TensorVec<Item = Y>,
71{
72 fn integrate(
73 &self,
74 mut function: impl FnMut(TensorRank0, &Y) -> Result<Y, IntegrationError>,
75 time: &[TensorRank0],
76 initial_condition: Y,
77 ) -> Result<(Vector, U, U), IntegrationError> {
78 let t_0 = time[0];
79 let t_f = time[time.len() - 1];
80 if time.len() < 2 {
81 return Err(IntegrationError::LengthTimeLessThanTwo);
82 } else if t_0 >= t_f {
83 return Err(IntegrationError::InitialTimeNotLessThanFinalTime);
84 }
85 let mut t = t_0;
86 let mut dt = t_f;
87 let mut e;
88 let mut k_1 = function(t, &initial_condition)?;
89 let mut k_2;
90 let mut k_3;
91 let mut k_4;
92 let mut t_sol = Vector::zero(0);
93 t_sol.push(t_0);
94 let mut y = initial_condition.clone();
95 let mut y_sol = U::zero(0);
96 y_sol.push(initial_condition.clone());
97 let mut dydt_sol = U::zero(0);
98 dydt_sol.push(k_1.clone());
99 let mut y_trial;
100 while t < t_f {
101 k_2 = function(t + 0.5 * dt, &(&k_1 * (0.5 * dt) + &y))?;
102 k_3 = function(t + 0.75 * dt, &(&k_2 * (0.75 * dt) + &y))?;
103 y_trial = (&k_1 * 2.0 + &k_2 * 3.0 + &k_3 * 4.0) * (dt / 9.0) + &y;
104 k_4 = function(t + dt, &y_trial)?;
105 e = ((&k_1 * -5.0 + k_2 * 6.0 + k_3 * 8.0 + &k_4 * -9.0) * (dt / 72.0)).norm_inf();
106 if e < self.abs_tol || e / y_trial.norm_inf() < self.rel_tol {
107 k_1 = k_4;
108 t += dt;
109 y = y_trial;
110 t_sol.push(t);
111 y_sol.push(y.clone());
112 dydt_sol.push(k_1.clone());
113 }
114 if e > 0.0 {
115 dt *= self.dt_beta * (self.abs_tol / e).powf(1.0 / self.dt_expn)
116 }
117 dt = dt.min(t_f - t)
118 }
119 if time.len() > 2 {
120 let t_int = Vector::new(time);
121 let (y_int, dydt_int) = self.interpolate(&t_int, &t_sol, &y_sol, function)?;
122 Ok((t_int, y_int, dydt_int))
123 } else {
124 Ok((t_sol, y_sol, dydt_sol))
125 }
126 }
127}
128
129impl<Y, U> InterpolateSolution<Y, U> for BogackiShampine
130where
131 Y: Tensor,
132 for<'a> &'a Y: Mul<TensorRank0, Output = Y> + Sub<&'a Y, Output = Y>,
133 U: TensorVec<Item = Y>,
134{
135 fn interpolate(
136 &self,
137 time: &Vector,
138 tp: &Vector,
139 yp: &U,
140 mut function: impl FnMut(TensorRank0, &Y) -> Result<Y, IntegrationError>,
141 ) -> Result<(U, U), IntegrationError> {
142 let mut dt;
143 let mut i;
144 let mut k_1;
145 let mut k_2;
146 let mut k_3;
147 let mut t;
148 let mut y;
149 let mut y_int = U::zero(0);
150 let mut dydt_int = U::zero(0);
151 let mut y_trial;
152 for time_k in time.iter() {
153 i = tp.iter().position(|tp_i| tp_i >= time_k).unwrap();
154 if time_k == &tp[i] {
155 t = tp[i];
156 y_trial = yp[i].clone();
157 dt = 0.0;
158 } else {
159 t = tp[i - 1];
160 y = yp[i - 1].clone();
161 dt = time_k - t;
162 k_1 = function(t, &y)?;
163 k_2 = function(t + 0.5 * dt, &(&k_1 * (0.5 * dt) + &y))?;
164 k_3 = function(t + 0.75 * dt, &(&k_2 * (0.75 * dt) + &y))?;
165 y_trial = (&k_1 * 2.0 + &k_2 * 3.0 + &k_3 * 4.0) * (dt / 9.0) + &y;
166 }
167 dydt_int.push(function(t + dt, &y_trial)?);
168 y_int.push(y_trial);
169 }
170 Ok((y_int, dydt_int))
171 }
172}