mathru/analysis/differential_equation/ordinary/solver/explicit/
adams_bashforth.rs

1//! Solves an ODE using Adams-Bashforth method.
2use crate::{
3    algebra::{abstr::Real, linear::vector::vector::Vector},
4    analysis::differential_equation::ordinary::{ExplicitInitialValueProblem, ExplicitODE},
5};
6
7#[cfg(feature = "serde")]
8use serde::{Deserialize, Serialize};
9use std::clone::Clone;
10
11/// Adams-Bashforth method
12///
13/// # Example
14///
15/// For this example, we want to solve the following ordinary differential
16/// equation:
17/// ```math
18/// \frac{dy}{dt} = ay = f(t, y)
19/// ```
20/// The initial condition is $y(0) = 0.5$ and we solve it in the interval
21/// $\lbrack 0, 2\rbrack$ The following equation is the closed solution for
22/// this ODE:
23/// ```math
24/// y(t) = C a e^{at}
25/// ```
26/// $C$ is a parameter and depends on the initial condition $y(t_{0})$
27/// ```math
28/// C = \frac{y(t_{0})}{ae^{at_{0}}}
29/// ```
30///
31/// In this example, we set $a=2$
32/// ```
33/// use mathru::{
34///     vector,
35///     algebra::linear::vector::Vector,
36///     analysis::differential_equation::ordinary::{ExplicitODE, ExplicitInitialValueProblemBuilder, solver::explicit::AdamsBashforth},
37/// };
38///
39/// pub struct Ode
40/// {
41/// }
42///
43/// impl ExplicitODE<f64> for Ode
44/// {
45///    fn ode(&self, _x: &f64, y: &Vector<f64>) -> Vector<f64>
46///    {
47///        y * &2.0f64
48///    }
49/// }
50///
51/// let ode = Ode{};
52///
53/// let problem = ExplicitInitialValueProblemBuilder::new(
54///     &ode,
55///     0.0,
56///     vector![0.5]
57/// )
58/// .t_end(2.0)
59/// .build();
60///
61/// let step_size: f64 = 0.001;
62/// let solver: AdamsBashforth<f64> = AdamsBashforth::new(1, step_size);
63///
64/// // Solve the ODE
65/// let (t, y): (Vec<f64>, Vec<Vector<f64>>) = solver.solve(&problem).unwrap();
66///
67/// ```
68#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
69#[derive(Clone, Copy, Debug)]
70pub struct AdamsBashforth<T> {
71    k: u8,
72    step_size: T,
73}
74
75impl<T> AdamsBashforth<T>
76where
77    T: Real,
78{
79    pub fn new(k: u8, step_size: T) -> AdamsBashforth<T> {
80        if k == 0 || k > 5 {
81            panic!();
82        }
83        if step_size <= T::zero() {
84            panic!();
85        }
86
87        AdamsBashforth { k, step_size }
88    }
89}
90
91impl<T> AdamsBashforth<T>
92where
93    T: Real,
94{
95    /// Solves `prob` using Adams' method.
96    ///
97    /// # Arguments
98    ///
99    /// # Return
100    ///
101    /// The solver returns a vector and a matrix, containing the times used in
102    /// each step of the algorithm and the respectful values for that time.
103    ///
104    /// # Panic
105    ///
106    pub fn solve<O>(
107        &self,
108        prob: &ExplicitInitialValueProblem<T, O>,
109    ) -> Result<(Vec<T>, Vec<Vector<T>>), String>
110    where
111        O: ExplicitODE<T>,
112    {
113        let t_start: T = prob.t_start();
114        let t_stop: T = prob.t_end().unwrap();
115
116        if t_start > t_stop {
117            panic!();
118        }
119
120        let mut x_n: Vector<T> = prob.init_cond();
121        let mut t_n: T = t_start;
122
123        let limit = ((t_stop - t_start) / self.step_size).ceil() + T::one();
124        let steps: usize = limit.to_u64() as usize;
125        let mut t_vec: Vec<T> = Vec::with_capacity(steps);
126        let mut res_vec: Vec<Vector<T>> = Vec::with_capacity(steps);
127
128        t_vec.push(t_n);
129        res_vec.push(x_n.clone());
130
131        //calculate initial steps
132        if self.k >= 2 {
133            let h: T = self.step_size.min(t_stop - t_n);
134            let x_n = AdamsBashforth::step_s1(prob, &t_vec, &res_vec, h);
135            t_n += h;
136
137            t_vec.push(t_n);
138            res_vec.push(x_n);
139        }
140
141        if self.k >= 3 {
142            let h: T = self.step_size.min(t_stop - t_n);
143            let x_n = AdamsBashforth::step_s2(prob, &t_vec, &res_vec, h);
144            t_n += h;
145
146            t_vec.push(t_n);
147            res_vec.push(x_n);
148        }
149
150        if self.k >= 4 {
151            let h: T = self.step_size.min(t_stop - t_n);
152            let x_n = AdamsBashforth::step_s3(prob, &t_vec, &res_vec, h);
153            t_n += h;
154
155            t_vec.push(t_n);
156            res_vec.push(x_n);
157        }
158
159        if self.k >= 5 {
160            let h: T = self.step_size.min(t_stop - t_n);
161            let x_n = AdamsBashforth::step_s4(prob, &t_vec, &res_vec, h);
162            t_n += h;
163
164            t_vec.push(t_n);
165            res_vec.push(x_n);
166        }
167
168        let step = match self.k {
169            1 => AdamsBashforth::step_s1,
170            2 => AdamsBashforth::step_s2,
171            3 => AdamsBashforth::step_s3,
172            4 => AdamsBashforth::step_s4,
173            5 => AdamsBashforth::step_s5,
174            _ => panic!(),
175        };
176
177        while (t_n - t_stop).abs() > T::from_f64(0.0000000001) {
178            //Step size
179            let h: T = self.step_size.min(t_stop - t_n);
180
181            x_n = step(prob, &t_vec, &res_vec, h);
182            t_n += h;
183
184            t_vec.push(t_n);
185            res_vec.push(x_n.clone());
186        }
187
188        Ok((t_vec, res_vec))
189    }
190}
191
192impl<T> AdamsBashforth<T>
193where
194    T: Real,
195{
196    fn step_s1<O>(
197        prob: &ExplicitInitialValueProblem<T, O>,
198        t: &[T],
199        x: &[Vector<T>],
200        h: T,
201    ) -> Vector<T>
202    where
203        O: ExplicitODE<T>,
204    {
205        let n: usize = x.len() - 1;
206        let x_n: &Vector<T> = &x[n];
207        let t_n: &T = &t[n];
208        let ode = prob.ode();
209        x_n + &(&ode.ode(t_n, x_n) * &h)
210    }
211
212    fn step_s2<O>(
213        prob: &ExplicitInitialValueProblem<T, O>,
214        t: &[T],
215        x: &[Vector<T>],
216        h: T,
217    ) -> Vector<T>
218    where
219        O: ExplicitODE<T>,
220    {
221        let n: usize = x.len() - 1;
222        let x_n: &Vector<T> = &x[n];
223        let t_n: &T = &t[n];
224        let x_n1: &Vector<T> = &x[n - 1];
225        let t_n1: &T = &t[n - 1];
226        let ode = prob.ode();
227        x_n + &((ode.ode(t_n, x_n) * T::from_f64(3.0 / 2.0)
228            + ode.ode(t_n1, x_n1) * T::from_f64(-0.5))
229            * h)
230    }
231
232    fn step_s3<O>(
233        prob: &ExplicitInitialValueProblem<T, O>,
234        t: &[T],
235        x: &[Vector<T>],
236        h: T,
237    ) -> Vector<T>
238    where
239        O: ExplicitODE<T>,
240    {
241        let n: usize = x.len() - 1;
242        let x_n: &Vector<T> = &x[n];
243        let t_n: &T = &t[n];
244        let x_n1: &Vector<T> = &x[n - 1];
245        let t_n1: &T = &t[n - 1];
246        let x_n2: &Vector<T> = &x[n - 2];
247        let t_n2: &T = &t[n - 2];
248        let ode = prob.ode();
249
250        x_n + &((ode.ode(t_n, x_n) * T::from_f64(23.0 / 12.0)
251            + ode.ode(t_n1, x_n1) * T::from_f64(-16.0 / 12.0)
252            + ode.ode(t_n2, x_n2) * T::from_f64(5.0 / 12.0))
253            * h)
254    }
255
256    fn step_s4<O>(
257        prob: &ExplicitInitialValueProblem<T, O>,
258        t: &[T],
259        x: &[Vector<T>],
260        h: T,
261    ) -> Vector<T>
262    where
263        O: ExplicitODE<T>,
264    {
265        let n: usize = x.len() - 1;
266        let x_n: &Vector<T> = &x[n];
267        let t_n: &T = &t[n];
268        let x_n1: &Vector<T> = &x[n - 1];
269        let t_n1: &T = &t[n - 1];
270        let x_n2: &Vector<T> = &x[n - 2];
271        let t_n2: &T = &t[n - 2];
272        let x_n3: &Vector<T> = &x[n - 3];
273        let t_n3: &T = &t[n - 3];
274        let ode = prob.ode();
275
276        x_n + &((ode.ode(t_n, x_n) * T::from_f64(55.0 / 24.0)
277            + ode.ode(t_n1, x_n1) * T::from_f64(-59.0 / 24.0)
278            + ode.ode(t_n2, x_n2) * T::from_f64(37.0 / 24.0)
279            + ode.ode(t_n3, x_n3) * T::from_f64(-9.0 / 24.0))
280            * h)
281    }
282
283    fn step_s5<O>(
284        prob: &ExplicitInitialValueProblem<T, O>,
285        t: &[T],
286        x: &[Vector<T>],
287        h: T,
288    ) -> Vector<T>
289    where
290        O: ExplicitODE<T>,
291    {
292        let n: usize = x.len() - 1;
293        let x_n: &Vector<T> = &x[n];
294        let t_n: &T = &t[n];
295        let x_n1: &Vector<T> = &x[n - 1];
296        let t_n1: &T = &t[n - 1];
297        let x_n2: &Vector<T> = &x[n - 2];
298        let t_n2: &T = &t[n - 2];
299        let x_n3: &Vector<T> = &x[n - 3];
300        let t_n3: &T = &t[n - 3];
301        let x_n4: &Vector<T> = &x[n - 4];
302        let t_n4: &T = &t[n - 4];
303        let ode = prob.ode();
304
305        x_n + &((ode.ode(t_n, x_n) * T::from_f64(1901.0 / 720.0)
306            + ode.ode(t_n1, x_n1) * T::from_f64(-2774.0 / 720.0)
307            + ode.ode(t_n2, x_n2) * T::from_f64(2616.0 / 720.0)
308            + ode.ode(t_n3, x_n3) * T::from_f64(-1274.0 / 720.0)
309            + ode.ode(t_n4, x_n4) * T::from_f64(251.0 / 720.0))
310            * h)
311    }
312}
mathru/analysis/differential_equation/ordinary/solver/explicit/adams_bashforth.rs

mathru/analysis/differential_equation/ordinary/solver/explicit/
adams_bashforth.rs
