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

bacon_sci_1/ivp/
rk.rs
