bacon_sci_1/ivp/
adams.rs

1/* This file is part of bacon.
2 * Copyright (c) Wyatt Campbell.
3 *
4 * See repository LICENSE for information.
5 */
6
7use super::{IVPSolver, IVPStatus};
8use nalgebra::{
9    allocator::Allocator, ComplexField, DefaultAllocator, DimName, RealField, VectorN, U3, U5,
10};
11use num_traits::{FromPrimitive, Zero};
12use std::collections::VecDeque;
13
14/// This trait allows a struct to be used in the Adams Predictor-Corrector solver
15///
16/// Things implementing AdamsSolver should have an AdamsInfo to handle the actual
17/// IVP solving.
18///
19/// # Examples
20/// See `struct Adams` for an example of implementing this trait
21pub trait AdamsSolver<N: ComplexField, S: DimName, O: DimName>: Sized
22where
23    DefaultAllocator: Allocator<N, S>,
24    DefaultAllocator: Allocator<N::RealField, O>,
25{
26    /// The polynomial interpolation coefficients for the predictor. Should start
27    /// with the coefficient for n - 1
28    fn predictor_coefficients() -> VectorN<N::RealField, O>;
29    /// The polynomial interpolation coefficients for the corrector. Must be
30    /// the same length as predictor_coefficients. Should start with the
31    /// implicit coefficient.
32    fn corrector_coefficients() -> VectorN<N::RealField, O>;
33    /// Coefficient for multiplying error by.
34    fn error_coefficient() -> N::RealField;
35
36    /// Use AdamsInfo to solve an initial value problem
37    fn solve_ivp<T: Clone, F: FnMut(N::RealField, &[N], &mut T) -> Result<VectorN<N, S>, String>>(
38        self,
39        f: F,
40        params: &mut T,
41    ) -> super::Path<N, N::RealField, S>;
42
43    /// Set the error tolerance for this solver.
44    fn with_tolerance(self, tol: N::RealField) -> Result<Self, String>;
45    /// Set the maximum time step for this solver.
46    fn with_dt_max(self, max: N::RealField) -> Result<Self, String>;
47    /// Set the minimum time step for this solver.
48    fn with_dt_min(self, min: N::RealField) -> Result<Self, String>;
49    /// Set the initial time for this solver.
50    fn with_start(self, t_initial: N::RealField) -> Result<Self, String>;
51    /// Set the end time for this solver.
52    fn with_end(self, t_final: N::RealField) -> Result<Self, String>;
53    /// Set the initial conditions for this solver.
54    fn with_initial_conditions(self, start: &[N]) -> Result<Self, String>;
55    /// Build this solver.
56    fn build(self) -> Self;
57}
58
59/// Provides an IVPSolver implementation for AdamsSolver, based on
60/// the predictor and correctorr coefficients. It is up to the AdamsSolver
61/// to set up AdamsInfo with the correct coefficients.
62#[derive(Debug, Clone)]
63#[cfg_attr(feature = "serialize", derive(Serialize, Deserialize))]
64pub struct AdamsInfo<N: ComplexField, S: DimName, O: DimName>
65where
66    DefaultAllocator: Allocator<N, S>,
67    DefaultAllocator: Allocator<N, O>,
68{
69    dt: Option<N::RealField>,
70    time: Option<N::RealField>,
71    end: Option<N::RealField>,
72    state: Option<VectorN<N, S>>,
73    dt_max: Option<N::RealField>,
74    dt_min: Option<N::RealField>,
75    tolerance: Option<N::RealField>,
76    predictor_coefficients: VectorN<N, O>,
77    corrector_coefficients: VectorN<N, O>,
78    error_coefficient: N::RealField,
79    memory: VecDeque<VectorN<N, S>>,
80    states: VecDeque<(N::RealField, VectorN<N, S>)>,
81    nflag: bool,
82    last: bool,
83}
84
85impl<N: ComplexField, S: DimName, O: DimName> AdamsInfo<N, S, O>
86where
87    DefaultAllocator: Allocator<N, S>,
88    DefaultAllocator: Allocator<N, O>,
89{
90    pub fn new() -> Self {
91        AdamsInfo {
92            dt: None,
93            time: None,
94            end: None,
95            state: None,
96            dt_max: None,
97            dt_min: None,
98            tolerance: None,
99            predictor_coefficients: VectorN::<N, O>::zero(),
100            corrector_coefficients: VectorN::<N, O>::zero(),
101            error_coefficient: N::RealField::zero(),
102            memory: VecDeque::new(),
103            states: VecDeque::new(),
104            nflag: false,
105            last: false,
106        }
107    }
108}
109
110#[allow(clippy::too_many_arguments)]
111fn rk4<
112    N: ComplexField,
113    S: DimName,
114    T: Clone,
115    F: FnMut(N::RealField, &[N], &mut T) -> Result<VectorN<N, S>, String>,
116>(
117    time: N::RealField,
118    dt: N::RealField,
119    initial: &[N],
120    states: &mut VecDeque<(N::RealField, VectorN<N, S>)>,
121    derivs: &mut VecDeque<VectorN<N, S>>,
122    mut f: F,
123    params: &mut T,
124    num: usize,
125) -> Result<(), String>
126where
127    DefaultAllocator: Allocator<N, S>,
128{
129    let mut state = VectorN::from_column_slice(initial);
130    let mut time = time;
131    for i in 0..num {
132        let k1 = f(time, state.as_slice(), &mut params.clone())? * N::from_real(dt);
133        let intermediate = &state + &k1 * N::from_f64(0.5).unwrap();
134        let k2 = f(
135            time + N::RealField::from_f64(0.5).unwrap() * dt,
136            intermediate.as_slice(),
137            &mut params.clone(),
138        )? * N::from_real(dt);
139        let intermediate = &state + &k2 * N::from_f64(0.5).unwrap();
140        let k3 = f(
141            time + N::RealField::from_f64(0.5).unwrap() * dt,
142            intermediate.as_slice(),
143            &mut params.clone(),
144        )? * N::from_real(dt);
145        let intermediate = &state + &k3;
146        let k4 = f(time + dt, intermediate.as_slice(), &mut params.clone())? * N::from_real(dt);
147        if i != 0 {
148            derivs.push_back(f(time, state.as_slice(), params)?);
149            states.push_back((time, state.clone()));
150        }
151        state += (k1 + k2 * N::from_f64(2.0).unwrap() + k3 * N::from_f64(2.0).unwrap() + k4)
152            * N::from_f64(1.0 / 6.0).unwrap();
153        time += dt;
154    }
155    derivs.push_back(f(time, state.as_slice(), params)?);
156    states.push_back((time, state));
157
158    Ok(())
159}
160
161impl<N: ComplexField, S: DimName, O: DimName> Default for AdamsInfo<N, S, O>
162where
163    DefaultAllocator: Allocator<N, S>,
164    DefaultAllocator: Allocator<N, O>,
165{
166    fn default() -> Self {
167        Self::new()
168    }
169}
170
171impl<N: ComplexField, S: DimName, O: DimName> IVPSolver<N, S> for AdamsInfo<N, S, O>
172where
173    DefaultAllocator: Allocator<N, S>,
174    DefaultAllocator: Allocator<N, O>,
175{
176    fn step<T: Clone, F: FnMut(N::RealField, &[N], &mut T) -> Result<VectorN<N, S>, String>>(
177        &mut self,
178        mut f: F,
179        params: &mut T,
180    ) -> Result<IVPStatus<N, S>, String> {
181        if self.time.unwrap() >= self.end.unwrap() {
182            return Ok(IVPStatus::Done);
183        }
184
185        let mut output = vec![];
186
187        if self.time.unwrap() + self.dt.unwrap() >= self.end.unwrap() {
188            self.dt = Some(self.end.unwrap() - self.time.unwrap());
189            rk4(
190                self.time.unwrap(),
191                self.dt.unwrap(),
192                self.state.as_ref().unwrap().as_slice(),
193                &mut self.states,
194                &mut self.memory,
195                &mut f,
196                params,
197                1,
198            )?;
199            *self.time.get_or_insert(N::RealField::zero()) += self.dt.unwrap();
200            return Ok(IVPStatus::Ok(vec![(
201                self.time.unwrap(),
202                self.states.back().unwrap().1.clone(),
203            )]));
204        }
205
206        if self.memory.is_empty() {
207            rk4(
208                self.time.unwrap(),
209                self.dt.unwrap(),
210                self.state.as_ref().unwrap().as_slice(),
211                &mut self.states,
212                &mut self.memory,
213                &mut f,
214                params,
215                self.predictor_coefficients.len() - 1,
216            )?;
217            self.time = Some(
218                self.time.unwrap()
219                    + N::RealField::from_usize(self.predictor_coefficients.len() - 1).unwrap()
220                        * self.dt.unwrap(),
221            );
222            self.state = Some(self.states.back().unwrap().1.clone());
223        }
224
225        let tenth_real = N::RealField::from_f64(0.1).unwrap();
226        let two_real = N::RealField::from_i32(2).unwrap();
227        let four_real = N::RealField::from_i32(4).unwrap();
228
229        let wp = &self.state.as_ref().unwrap();
230        let wp =
231            VectorN::<N, S>::from_iterator(wp.as_slice().iter().enumerate().map(|(ind, y)| {
232                let mut acc = *y;
233                let dt = N::from_real(self.dt.unwrap());
234                for (j, mem) in self.memory.iter().enumerate() {
235                    acc += mem[ind]
236                        * self.predictor_coefficients[self.predictor_coefficients.len() - j - 2]
237                        * dt;
238                }
239                acc
240            }));
241
242        let implicit = f(
243            self.time.unwrap() + self.dt.unwrap(),
244            wp.as_slice(),
245            &mut params.clone(),
246        )?;
247        let wc = &self.state.as_ref().unwrap();
248        let wc =
249            VectorN::<N, S>::from_iterator(wc.as_slice().iter().enumerate().map(|(ind, y)| {
250                let dt = N::from_real(self.dt.unwrap());
251                let mut acc = implicit[ind] * self.corrector_coefficients[0] * dt;
252                for (j, mem) in self.memory.iter().enumerate() {
253                    acc += mem[ind]
254                        * self.corrector_coefficients[self.corrector_coefficients.len() - j - 1]
255                        * dt;
256                }
257                acc + *y
258            }));
259
260        let diff = &wc - &wp;
261        let error = self.error_coefficient / self.dt.unwrap() * diff.dot(&diff).sqrt().abs();
262
263        if error <= self.tolerance.unwrap() {
264            self.state = Some(wc);
265            self.time = Some(self.time.unwrap() + self.dt.unwrap());
266            if self.nflag {
267                for state in self.states.iter() {
268                    output.push((state.0, state.1.clone()));
269                }
270                self.nflag = false;
271            }
272            output.push((self.time.unwrap(), self.state.as_ref().unwrap().clone()));
273
274            self.memory.push_back(implicit);
275            self.states
276                .push_back((self.time.unwrap(), self.state.as_ref().unwrap().clone()));
277            self.memory.pop_front();
278            self.states.pop_front();
279
280            if self.last {
281                return Ok(IVPStatus::Ok(output));
282            }
283
284            if error < tenth_real * self.tolerance.unwrap()
285                || self.time.unwrap() > self.end.unwrap()
286            {
287                let q = (self.tolerance.unwrap() / (two_real * error)).powf(
288                    N::RealField::from_f64(1.0 / self.predictor_coefficients.len() as f64).unwrap(),
289                );
290                if q > four_real {
291                    self.dt = Some(self.dt.unwrap() * four_real);
292                } else {
293                    self.dt = Some(self.dt.unwrap() * q);
294                }
295
296                if self.dt.unwrap() > self.dt_max.unwrap() {
297                    self.dt = Some(self.dt_max.unwrap());
298                }
299
300                if self.time.unwrap()
301                    + N::RealField::from_usize(self.predictor_coefficients.len()).unwrap()
302                        * self.dt.unwrap()
303                    > self.end.unwrap()
304                {
305                    self.dt = Some(
306                        (self.end.unwrap() - self.time.unwrap())
307                            / N::RealField::from_usize(self.predictor_coefficients.len()).unwrap(),
308                    );
309                    self.last = true;
310                }
311
312                self.memory.clear();
313                self.states.clear();
314            }
315
316            return Ok(IVPStatus::Ok(output));
317        }
318
319        let q = (self.tolerance.unwrap() / (N::RealField::from_f64(2.0).unwrap() * error)).powf(
320            N::RealField::from_f64(1.0 / (self.predictor_coefficients.len() as f64)).unwrap(),
321        );
322
323        if q < tenth_real {
324            self.dt = Some(self.dt.unwrap() * tenth_real);
325        } else {
326            self.dt = Some(self.dt.unwrap() * q);
327        }
328
329        if self.dt.unwrap() < self.dt_min.unwrap() {
330            return Err("AdamsInfo step: minimum dt exceeded".to_owned());
331        }
332
333        self.memory.clear();
334        self.states.clear();
335        Ok(IVPStatus::Redo)
336    }
337
338    fn with_tolerance(mut self, tol: N::RealField) -> Result<Self, String> {
339        if !tol.is_sign_positive() {
340            return Err("AdamsInfo with_tolerance: tolerance must be postive".to_owned());
341        }
342        self.tolerance = Some(tol);
343        Ok(self)
344    }
345
346    fn with_dt_max(mut self, max: N::RealField) -> Result<Self, String> {
347        if !max.is_sign_positive() {
348            return Err("AdamsInfo with_dt_max: dt_max must be positive".to_owned());
349        }
350        if let Some(min) = self.dt_min {
351            if max <= min {
352                return Err("AdamsInfo with_dt_max: dt_max must be greater than dt_min".to_owned());
353            }
354        }
355        self.dt_max = Some(max);
356        self.dt = Some(max);
357        Ok(self)
358    }
359
360    fn with_dt_min(mut self, min: N::RealField) -> Result<Self, String> {
361        if !min.is_sign_positive() {
362            return Err("AdamsInfo with_dt_min: dt_min must be positive".to_owned());
363        }
364        if let Some(max) = self.dt_max {
365            if min >= max {
366                return Err("AdamsInfo with_dt_min: dt_min must be less than dt_max".to_owned());
367            }
368        }
369        self.dt_min = Some(min);
370        Ok(self)
371    }
372
373    fn with_start(mut self, t_initial: N::RealField) -> Result<Self, String> {
374        if let Some(end) = self.end {
375            if end <= t_initial {
376                return Err("AdamsInfo with_start: Start must be before end".to_owned());
377            }
378        }
379        self.time = Some(t_initial);
380        Ok(self)
381    }
382
383    fn with_end(mut self, t_final: N::RealField) -> Result<Self, String> {
384        if let Some(start) = self.time {
385            if t_final <= start {
386                return Err("AdamsInfo with_end: Start must be before end".to_owned());
387            }
388        }
389        self.end = Some(t_final);
390        Ok(self)
391    }
392
393    fn with_initial_conditions(mut self, start: &[N]) -> Result<Self, String> {
394        self.state = Some(VectorN::<N, S>::from_column_slice(start));
395        Ok(self)
396    }
397
398    fn build(self) -> Self {
399        self
400    }
401
402    fn get_initial_conditions(&self) -> Option<VectorN<N, S>> {
403        if let Some(state) = &self.state {
404            Some(state.clone())
405        } else {
406            None
407        }
408    }
409
410    fn get_time(&self) -> Option<N::RealField> {
411        if let Some(time) = &self.time {
412            Some(*time)
413        } else {
414            None
415        }
416    }
417
418    fn check_start(&self) -> Result<(), String> {
419        if self.time == None {
420            Err("AdamsInfo check_start: No initial time".to_owned())
421        } else if self.end == None {
422            Err("AdamsInfo check_start: No end time".to_owned())
423        } else if self.tolerance == None {
424            Err("AdamsInfo check_start: No tolerance".to_owned())
425        } else if self.state == None {
426            Err("AdamsInfo check_start: No initial conditions".to_owned())
427        } else if self.dt_max == None {
428            Err("AdamsInfo check_start: No dt_max".to_owned())
429        } else if self.dt_min == None {
430            Err("AdamsInfo check_start: No dt_min".to_owned())
431        } else {
432            Ok(())
433        }
434    }
435}
436
437/// 5th order Adams predictor-corrector method for solving an IVP.
438///
439/// Defines the predictor and corrector coefficients, as well as
440/// the error coefficient. Uses AdamsInfo for the actual solving.
441///
442/// # Examples
443/// ```
444/// use nalgebra::{VectorN, U1};
445/// use bacon_sci::ivp::{Adams, AdamsSolver};
446/// fn derivatives(_t: f64, state: &[f64], _p: &mut ()) -> Result<VectorN<f64, U1>, String> {
447///     Ok(VectorN::<f64, U1>::from_column_slice(state))
448/// }
449///
450///
451/// fn example() -> Result<(), String> {
452///     let adams = Adams::new()
453///         .with_dt_max(0.1)?
454///         .with_dt_min(0.00001)?
455///         .with_tolerance(0.00001)?
456///         .with_start(0.0)?
457///         .with_end(1.0)?
458///         .with_initial_conditions(&[1.0])?
459///         .build();
460///     let path = adams.solve_ivp(derivatives, &mut ())?;
461///     for (time, state) in &path {
462///         assert!((time.exp() - state.column(0)[0]).abs() < 0.001);
463///     }
464///     Ok(())
465/// }
466/// ```
467#[derive(Debug, Clone)]
468#[cfg_attr(feature = "serialize", derive(Serialize, Deserialize))]
469pub struct Adams<N: ComplexField, S: DimName>
470where
471    DefaultAllocator: Allocator<N, S>,
472    DefaultAllocator: Allocator<N::RealField, U5>,
473    DefaultAllocator: Allocator<N, U5>,
474{
475    info: AdamsInfo<N, S, U5>,
476}
477
478impl<N: ComplexField, S: DimName> Adams<N, S>
479where
480    DefaultAllocator: Allocator<N, S>,
481    DefaultAllocator: Allocator<N, U5>,
482    DefaultAllocator: Allocator<N::RealField, U5>,
483{
484    pub fn new() -> Self {
485        let mut info = AdamsInfo::new();
486        info.corrector_coefficients = VectorN::<N, U5>::from_iterator(
487            Self::corrector_coefficients()
488                .iter()
489                .map(|&x| N::from_real(x)),
490        );
491        info.predictor_coefficients = VectorN::<N, U5>::from_iterator(
492            Self::predictor_coefficients()
493                .iter()
494                .map(|&x| N::from_real(x)),
495        );
496        info.error_coefficient = Self::error_coefficient();
497
498        Adams { info }
499    }
500}
501
502impl<N: ComplexField, S: DimName> Default for Adams<N, S>
503where
504    DefaultAllocator: Allocator<N, S>,
505    DefaultAllocator: Allocator<N::RealField, U5>,
506    DefaultAllocator: Allocator<N, U5>,
507{
508    fn default() -> Self {
509        Self::new()
510    }
511}
512
513impl<N: ComplexField, S: DimName> AdamsSolver<N, S, U5> for Adams<N, S>
514where
515    DefaultAllocator: Allocator<N, S>,
516    DefaultAllocator: Allocator<N::RealField, U5>,
517    DefaultAllocator: Allocator<N, U5>,
518{
519    fn predictor_coefficients() -> VectorN<N::RealField, U5> {
520        VectorN::<N::RealField, U5>::from_column_slice(&[
521            N::RealField::from_f64(55.0 / 24.0).unwrap(),
522            N::RealField::from_f64(-59.0 / 24.0).unwrap(),
523            N::RealField::from_f64(37.0 / 24.0).unwrap(),
524            N::RealField::from_f64(-9.0 / 24.0).unwrap(),
525            N::RealField::zero(),
526        ])
527    }
528
529    fn corrector_coefficients() -> VectorN<N::RealField, U5> {
530        VectorN::<N::RealField, U5>::from_column_slice(&[
531            N::RealField::from_f64(251.0 / 720.0).unwrap(),
532            N::RealField::from_f64(646.0 / 720.0).unwrap(),
533            N::RealField::from_f64(-264.0 / 720.0).unwrap(),
534            N::RealField::from_f64(106.0 / 720.0).unwrap(),
535            N::RealField::from_f64(-19.0 / 720.0).unwrap(),
536        ])
537    }
538
539    fn error_coefficient() -> N::RealField {
540        N::RealField::from_f64(19.0 / 270.0).unwrap()
541    }
542
543    fn solve_ivp<
544        T: Clone,
545        F: FnMut(N::RealField, &[N], &mut T) -> Result<VectorN<N, S>, String>,
546    >(
547        self,
548        f: F,
549        params: &mut T,
550    ) -> super::Path<N, N::RealField, S> {
551        self.info.solve_ivp(f, params)
552    }
553
554    fn with_tolerance(mut self, tol: N::RealField) -> Result<Self, String> {
555        self.info = self.info.with_tolerance(tol)?;
556        Ok(self)
557    }
558
559    fn with_dt_max(mut self, max: N::RealField) -> Result<Self, String> {
560        self.info = self.info.with_dt_max(max)?;
561        Ok(self)
562    }
563
564    fn with_dt_min(mut self, min: N::RealField) -> Result<Self, String> {
565        self.info = self.info.with_dt_min(min)?;
566        Ok(self)
567    }
568
569    fn with_start(mut self, t_initial: N::RealField) -> Result<Self, String> {
570        self.info = self.info.with_start(t_initial)?;
571        Ok(self)
572    }
573
574    fn with_end(mut self, t_final: N::RealField) -> Result<Self, String> {
575        self.info = self.info.with_end(t_final)?;
576        Ok(self)
577    }
578
579    fn with_initial_conditions(mut self, start: &[N]) -> Result<Self, String> {
580        self.info = self.info.with_initial_conditions(start)?;
581        Ok(self)
582    }
583
584    fn build(mut self) -> Self {
585        self.info = self.info.build();
586        self
587    }
588}
589
590impl<N: ComplexField, S: DimName> From<Adams<N, S>> for AdamsInfo<N, S, U5>
591where
592    DefaultAllocator: Allocator<N, S>,
593    DefaultAllocator: Allocator<N, U5>,
594    DefaultAllocator: Allocator<N::RealField, U5>,
595{
596    fn from(adams: Adams<N, S>) -> AdamsInfo<N, S, U5> {
597        adams.info
598    }
599}
600
601/// 3rd order Adams predictor-corrector method for solving an IVP.
602///
603/// Defines the predictor and corrector coefficients, as well as
604/// the error coefficient. Uses AdamsInfo for the actual solving.
605///
606/// # Examples
607/// ```
608/// use nalgebra::{VectorN, U1};
609/// use bacon_sci::ivp::{Adams2, AdamsSolver};
610/// fn derivatives(_t: f64, state: &[f64], _p: &mut ()) -> Result<VectorN<f64, U1>, String> {
611///     Ok(VectorN::<f64, U1>::from_column_slice(state))
612/// }
613///
614///
615/// fn example() -> Result<(), String> {
616///     let adams = Adams2::new()
617///         .with_dt_max(0.1)?
618///         .with_dt_min(0.00001)?
619///         .with_tolerance(0.00001)?
620///         .with_start(0.0)?
621///         .with_end(1.0)?
622///         .with_initial_conditions(&[1.0])?
623///         .build();
624///     let path = adams.solve_ivp(derivatives, &mut ())?;
625///     for (time, state) in &path {
626///         assert!((time.exp() - state.column(0)[0]).abs() < 0.001);
627///     }
628///     Ok(())
629/// }
630/// ```
631#[derive(Debug, Clone)]
632#[cfg_attr(feature = "serialize", derive(Serialize, Deserialize))]
633pub struct Adams2<N: ComplexField, S: DimName>
634where
635    DefaultAllocator: Allocator<N, S>,
636    DefaultAllocator: Allocator<N::RealField, U3>,
637    DefaultAllocator: Allocator<N, U3>,
638{
639    info: AdamsInfo<N, S, U3>,
640}
641
642impl<N: ComplexField, S: DimName> Adams2<N, S>
643where
644    DefaultAllocator: Allocator<N, S>,
645    DefaultAllocator: Allocator<N, U3>,
646    DefaultAllocator: Allocator<N::RealField, U3>,
647{
648    pub fn new() -> Self {
649        let mut info = AdamsInfo::new();
650        info.corrector_coefficients = VectorN::<N, U3>::from_iterator(
651            Self::corrector_coefficients()
652                .iter()
653                .map(|&x| N::from_real(x)),
654        );
655        info.predictor_coefficients = VectorN::<N, U3>::from_iterator(
656            Self::predictor_coefficients()
657                .iter()
658                .map(|&x| N::from_real(x)),
659        );
660        info.error_coefficient = Self::error_coefficient();
661
662        Adams2 { info }
663    }
664}
665
666impl<N: ComplexField, S: DimName> Default for Adams2<N, S>
667where
668    DefaultAllocator: Allocator<N, S>,
669    DefaultAllocator: Allocator<N::RealField, U3>,
670    DefaultAllocator: Allocator<N, U3>,
671{
672    fn default() -> Self {
673        Self::new()
674    }
675}
676
677impl<N: ComplexField, S: DimName> AdamsSolver<N, S, U3> for Adams2<N, S>
678where
679    DefaultAllocator: Allocator<N, S>,
680    DefaultAllocator: Allocator<N::RealField, U3>,
681    DefaultAllocator: Allocator<N, U3>,
682{
683    fn predictor_coefficients() -> VectorN<N::RealField, U3> {
684        VectorN::<N::RealField, U3>::from_column_slice(&[
685            N::RealField::from_f64(1.5).unwrap(),
686            N::RealField::from_f64(-0.5).unwrap(),
687            N::RealField::zero(),
688        ])
689    }
690
691    fn corrector_coefficients() -> VectorN<N::RealField, U3> {
692        VectorN::<N::RealField, U3>::from_column_slice(&[
693            N::RealField::from_f64(5.0 / 12.0).unwrap(),
694            N::RealField::from_f64(2.0 / 3.0).unwrap(),
695            N::RealField::from_f64(-1.0 / 12.0).unwrap(),
696        ])
697    }
698
699    fn error_coefficient() -> N::RealField {
700        N::RealField::from_f64(19.0 / 270.0).unwrap()
701    }
702
703    fn solve_ivp<
704        T: Clone,
705        F: FnMut(N::RealField, &[N], &mut T) -> Result<VectorN<N, S>, String>,
706    >(
707        self,
708        f: F,
709        params: &mut T,
710    ) -> super::Path<N, N::RealField, S> {
711        self.info.solve_ivp(f, params)
712    }
713
714    fn with_tolerance(mut self, tol: N::RealField) -> Result<Self, String> {
715        self.info = self.info.with_tolerance(tol)?;
716        Ok(self)
717    }
718
719    fn with_dt_max(mut self, max: N::RealField) -> Result<Self, String> {
720        self.info = self.info.with_dt_max(max)?;
721        Ok(self)
722    }
723
724    fn with_dt_min(mut self, min: N::RealField) -> Result<Self, String> {
725        self.info = self.info.with_dt_min(min)?;
726        Ok(self)
727    }
728
729    fn with_start(mut self, t_initial: N::RealField) -> Result<Self, String> {
730        self.info = self.info.with_start(t_initial)?;
731        Ok(self)
732    }
733
734    fn with_end(mut self, t_final: N::RealField) -> Result<Self, String> {
735        self.info = self.info.with_end(t_final)?;
736        Ok(self)
737    }
738
739    fn with_initial_conditions(mut self, start: &[N]) -> Result<Self, String> {
740        self.info = self.info.with_initial_conditions(start)?;
741        Ok(self)
742    }
743
744    fn build(mut self) -> Self {
745        self.info = self.info.build();
746        self
747    }
748}
749
750impl<N: ComplexField, S: DimName> From<Adams2<N, S>> for AdamsInfo<N, S, U3>
751where
752    DefaultAllocator: Allocator<N, S>,
753    DefaultAllocator: Allocator<N, U3>,
754    DefaultAllocator: Allocator<N::RealField, U3>,
755{
756    fn from(adams: Adams2<N, S>) -> AdamsInfo<N, S, U3> {
757        adams.info
758    }
759}
bacon_sci_1/ivp/adams.rs

bacon_sci_1/ivp/
adams.rs
