differential_equations/ode/methods/runge_kutta/explicit/
dop853.rs

1//! DOP853 ODENumericalMethod for Ordinary Differential Equations.
2
3use crate::{
4    Error, Status,
5    alias::Evals,
6    interpolate::Interpolation,
7    ode::{ODENumericalMethod, ODE, methods::h_init},
8    traits::{CallBackData, Real, State},
9    utils::{constrain_step_size, validate_step_size_parameters},
10};
11
12/// Dormand Prince 8(5, 3) Method for solving ordinary differential equations.
13/// 8th order Dormand Prince method with embedded 5th order error estimation and 3rd order interpolation.
14/// The resulting interpolant is of order 7.
15///
16/// Builds should begin with weight, normal, dense, or even methods.
17/// and then chain the other methods to set the parameters.
18/// The defaults should be great for most cases.
19///
20/// # Example
21/// ```
22/// use differential_equations::prelude::*;
23/// use nalgebra::{SVector, vector};
24///
25/// let mut dop853 = DOP853::new()
26///    .rtol(1e-12)
27///    .atol(1e-12);
28///
29/// let t0 = 0.0;
30/// let tf = 10.0;
31/// let y0 = vector![1.0, 0.0];
32/// struct Example;
33/// impl ODE<f64, SVector<f64, 2>> for Example {
34///    fn diff(&self, _t: f64, y: &SVector<f64, 2>, dydt: &mut SVector<f64, 2>) {
35///       dydt[0] = y[1];
36///       dydt[1] = -y[0];
37///   }
38/// }
39/// let solution = ODEProblem::new(Example, t0, tf, y0).solve(&mut dop853).unwrap();
40///
41/// let (t, y) = solution.last().unwrap();
42/// println!("Solution: ({}, {})", t, y);
43/// ```
44///
45/// # Settings
46/// * `rtol`   - Relative tolerance for the solver.
47/// * `atol`   - Absolute tolerance for the solver.
48/// * `h0`     - Initial step size.
49/// * `h_max`   - Maximum step size for the solver.
50/// * `max_steps` - Maximum number of steps for the solver.
51/// * `n_stiff` - Number of steps to check for stiffness.
52/// * `safe`   - Safety factor for step size prediction.
53/// * `fac1`   - Parameter for step size selection.
54/// * `fac2`   - Parameter for step size selection.
55/// * `beta`   - Beta for stabilized step size control.
56///
57/// # Default Settings
58/// * `rtol`   - 1e-3
59/// * `atol`   - 1e-6
60/// * `h0`     - None (Calculated by solver if None)
61/// * `h_max`   - None (Calculated by tf - t0 if None)
62/// * `h_min`   - 0.0
63/// * `max_steps` - 1_000_000
64/// * `n_stiff` - 100
65/// * `safe`   - 0.9
66/// * `fac1`   - 0.33
67/// * `fac2`   - 6.0
68/// * `beta`   - 0.0
69///
70pub struct DOP853<T: Real, V: State<T>, D: CallBackData> {
71    // Initial Conditions
72    pub h0: T, // Initial Step Size
73
74    // Current iteration
75    t: T,
76    y: V,
77    h: T,
78
79    // Tolerances
80    pub rtol: T,
81    pub atol: T,
82
83    // Settings
84    pub h_max: T,
85    pub h_min: T,
86    pub max_steps: usize,
87    pub n_stiff: usize,
88
89    // DOP853 Specific Settings
90    pub safe: T,
91    pub fac1: T,
92    pub fac2: T,
93    pub beta: T,
94
95    // Derived Settings
96    expo1: T,
97    facc1: T,
98    facc2: T,
99    facold: T,
100    fac11: T,
101    fac: T,
102
103    // Iteration Tracking
104    status: Status<T, V, D>,
105    steps: usize,      // Number of Steps
106    n_accepted: usize, // Number of Accepted Steps
107
108    // Stiffness Detection
109    h_lamb: T,
110    non_stiff_counter: usize,
111    stiffness_counter: usize,
112
113    // Butcher tableau coefficients (converted to type T)
114    a: [[T; 12]; 12],
115    b: [T; 12],
116    c: [T; 12],
117    er: [T; 12],
118    bhh: [T; 3],
119
120    // Dense output coefficients
121    a_dense: [[T; 16]; 3],
122    c_dense: [T; 3],
123    b_dense: [[T; 16]; 4],
124
125    // Derivatives - using array instead of individually numbered variables
126    k: [V; 12], // k[0] is derivative at t, others are stage derivatives
127
128    // For Interpolation - using array instead of individually numbered variables
129    y_old: V,     // State at Previous Step
130    t_old: T,     // Time of Previous Step
131    h_old: T,     // Step Size of Previous Step
132    cont: [V; 8], // Interpolation coefficients
133}
134
135impl<T: Real, V: State<T>, D: CallBackData> ODENumericalMethod<T, V, D> for DOP853<T, V, D> {
136    fn init<F>(&mut self, ode: &F, t0: T, tf: T, y0: &V) -> Result<Evals, Error<T, V>>
137    where
138        F: ODE<T, V, D>,
139    {
140        let mut evals = Evals::new();
141
142        // Set Current State as Initial State
143        self.t = t0;
144        self.y = *y0;
145
146        // Calculate derivative at t0
147        ode.diff(t0, y0, &mut self.k[0]);
148        evals.fcn += 1; // Increment function evaluations for initial derivative calculation
149
150        // Initialize Previous State
151        self.t_old = self.t;
152        self.y_old = self.y;
153
154        // Calculate Initial Step
155        if self.h0 == T::zero() {
156            self.h0 = h_init(
157                ode, t0, tf, y0, 8, self.rtol, self.atol, self.h_min, self.h_max,
158            );
159            evals.fcn += 1; // Increment function evaluations for initial step size calculation
160
161            // Adjust h0 to be within bounds
162            let posneg = (tf - t0).signum();
163            if self.h0.abs() < self.h_min.abs() {
164                self.h0 = self.h_min.abs() * posneg;
165            } else if self.h0.abs() > self.h_max.abs() {
166                self.h0 = self.h_max.abs() * posneg;
167            }
168        }
169
170        // Check if h0 is within bounds, and h_min and h_max are valid
171        match validate_step_size_parameters::<T, V, D>(self.h0, self.h_min, self.h_max, t0, tf) {
172            Ok(h0) => self.h = h0,
173            Err(status) => return Err(status),
174        }
175
176        // Make sure iteration variables are reset
177        self.h_lamb = T::zero();
178        self.non_stiff_counter = 0;
179        self.stiffness_counter = 0;
180
181        // ODENumericalMethod is ready to go
182        self.status = Status::Initialized;
183
184        Ok(evals)
185    }
186
187    fn step<F>(&mut self, ode: &F) -> Result<Evals, Error<T, V>>
188    where
189        F: ODE<T, V, D>,
190    {
191        let mut evals = Evals::new();
192
193        // Check if Max Steps Reached
194        if self.steps >= self.max_steps {
195            self.status = Status::Error(Error::MaxSteps {
196                t: self.t,
197                y: self.y,
198            });
199            return Err(Error::MaxSteps {
200                t: self.t,
201                y: self.y,
202            });
203        }
204
205        // Check if Step Size is too smaller then machine default_epsilon
206        if self.h.abs() < T::default_epsilon() {
207            self.status = Status::Error(Error::StepSize {
208                t: self.t,
209                y: self.y,
210            });
211            return Err(Error::StepSize {
212                t: self.t,
213                y: self.y,
214            });
215        }
216
217        // The twelve stages
218        ode.diff(
219            self.t + self.c[1] * self.h,
220            &(self.y + self.k[0] * (self.a[1][0] * self.h)),
221            &mut self.k[1],
222        );
223        ode.diff(
224            self.t + self.c[2] * self.h,
225            &(self.y + self.k[0] * (self.a[2][0] * self.h) + self.k[1] * (self.a[2][1] * self.h)),
226            &mut self.k[2],
227        );
228        ode.diff(
229            self.t + self.c[3] * self.h,
230            &(self.y + self.k[0] * (self.a[3][0] * self.h) + self.k[2] * (self.a[3][2] * self.h)),
231            &mut self.k[3],
232        );
233        ode.diff(
234            self.t + self.c[4] * self.h,
235            &(self.y
236                + self.k[0] * (self.a[4][0] * self.h)
237                + self.k[2] * (self.a[4][2] * self.h)
238                + self.k[3] * (self.a[4][3] * self.h)),
239            &mut self.k[4],
240        );
241        ode.diff(
242            self.t + self.c[5] * self.h,
243            &(self.y
244                + self.k[0] * (self.a[5][0] * self.h)
245                + self.k[3] * (self.a[5][3] * self.h)
246                + self.k[4] * (self.a[5][4] * self.h)),
247            &mut self.k[5],
248        );
249        ode.diff(
250            self.t + self.c[6] * self.h,
251            &(self.y
252                + self.k[0] * (self.a[6][0] * self.h)
253                + self.k[3] * (self.a[6][3] * self.h)
254                + self.k[4] * (self.a[6][4] * self.h)
255                + self.k[5] * (self.a[6][5] * self.h)),
256            &mut self.k[6],
257        );
258        ode.diff(
259            self.t + self.c[7] * self.h,
260            &(self.y
261                + self.k[0] * (self.a[7][0] * self.h)
262                + self.k[3] * (self.a[7][3] * self.h)
263                + self.k[4] * (self.a[7][4] * self.h)
264                + self.k[5] * (self.a[7][5] * self.h)
265                + self.k[6] * (self.a[7][6] * self.h)),
266            &mut self.k[7],
267        );
268        ode.diff(
269            self.t + self.c[8] * self.h,
270            &(self.y
271                + self.k[0] * (self.a[8][0] * self.h)
272                + self.k[3] * (self.a[8][3] * self.h)
273                + self.k[4] * (self.a[8][4] * self.h)
274                + self.k[5] * (self.a[8][5] * self.h)
275                + self.k[6] * (self.a[8][6] * self.h)
276                + self.k[7] * (self.a[8][7] * self.h)),
277            &mut self.k[8],
278        );
279        ode.diff(
280            self.t + self.c[9] * self.h,
281            &(self.y
282                + self.k[0] * (self.a[9][0] * self.h)
283                + self.k[3] * (self.a[9][3] * self.h)
284                + self.k[4] * (self.a[9][4] * self.h)
285                + self.k[5] * (self.a[9][5] * self.h)
286                + self.k[6] * (self.a[9][6] * self.h)
287                + self.k[7] * (self.a[9][7] * self.h)
288                + self.k[8] * (self.a[9][8] * self.h)),
289            &mut self.k[9],
290        );
291        ode.diff(
292            self.t + self.c[10] * self.h,
293            &(self.y
294                + self.k[0] * (self.a[10][0] * self.h)
295                + self.k[3] * (self.a[10][3] * self.h)
296                + self.k[4] * (self.a[10][4] * self.h)
297                + self.k[5] * (self.a[10][5] * self.h)
298                + self.k[6] * (self.a[10][6] * self.h)
299                + self.k[7] * (self.a[10][7] * self.h)
300                + self.k[8] * (self.a[10][8] * self.h)
301                + self.k[9] * (self.a[10][9] * self.h)),
302            &mut self.k[1],
303        );
304        let t_new = self.t + self.h;
305        let yy1 = self.y
306            + self.k[0] * (self.a[11][0] * self.h)
307            + self.k[3] * (self.a[11][3] * self.h)
308            + self.k[4] * (self.a[11][4] * self.h)
309            + self.k[5] * (self.a[11][5] * self.h)
310            + self.k[6] * (self.a[11][6] * self.h)
311            + self.k[7] * (self.a[11][7] * self.h)
312            + self.k[8] * (self.a[11][8] * self.h)
313            + self.k[9] * (self.a[11][9] * self.h)
314            + self.k[1] * (self.a[11][10] * self.h);
315        ode.diff(t_new, &yy1, &mut self.k[2]);
316        self.k[3] = self.k[0] * self.b[0]
317            + self.k[5] * self.b[5]
318            + self.k[6] * self.b[6]
319            + self.k[7] * self.b[7]
320            + self.k[8] * self.b[8]
321            + self.k[9] * self.b[9]
322            + self.k[1] * self.b[10]
323            + self.k[2] * self.b[11];
324        self.k[4] = self.y + self.k[3] * self.h;
325
326        // Log evals
327        evals.fcn += 11; // Increment function evaluations for the 11 derivative calculations
328
329        // Error Estimation
330        let mut err = T::zero();
331        let mut err2 = T::zero();
332
333        let n = self.y.len();
334        for i in 0..n {
335            let sk = self.atol + self.rtol * self.y.get(i).abs().max(self.k[4].get(i).abs());
336            let erri = self.k[3].get(i)
337                - self.bhh[0] * self.k[0].get(i)
338                - self.bhh[1] * self.k[8].get(i)
339                - self.bhh[2] * self.k[2].get(i);
340            err2 += (erri / sk).powi(2);
341            let erri = self.er[0] * self.k[0].get(i)
342                + self.er[5] * self.k[5].get(i)
343                + self.er[6] * self.k[6].get(i)
344                + self.er[7] * self.k[7].get(i)
345                + self.er[8] * self.k[8].get(i)
346                + self.er[9] * self.k[9].get(i)
347                + self.er[10] * self.k[1].get(i)
348                + self.er[11] * self.k[2].get(i);
349            err += (erri / sk).powi(2);
350        }
351        let mut deno = err + T::from_f64(0.01).unwrap() * err2;
352        if deno <= T::zero() {
353            deno = T::one();
354        }
355        err = self.h.abs() * err * (T::one() / (deno * T::from_usize(n).unwrap())).sqrt();
356
357        // Computation of h_new
358        self.fac11 = err.powf(self.expo1);
359        // Lund-stabilization
360        self.fac = self.fac11 / self.facold.powf(self.beta);
361        // Requirement that fac1 <= h_new/h <= fac2
362        self.fac = self.facc2.max(self.facc1.min(self.fac / self.safe));
363        let mut h_new = self.h / self.fac;
364
365        if err <= T::one() {
366            // Step Accepted
367            self.facold = err.max(T::from_f64(1.0e-4).unwrap());
368            let y_new = self.k[4];
369            ode.diff(t_new, &y_new, &mut self.k[3]);
370            evals.fcn += 1; // Increment function evaluations for the new derivative calculation
371
372            // stiffness detection
373            if self.n_accepted % self.n_stiff == 0 {
374                let mut stdnum = T::zero();
375                let mut stden = T::zero();
376                let sqr = self.k[3] - self.k[2];
377                for i in 0..sqr.len() {
378                    stdnum += sqr.get(i).powi(2);
379                }
380                let sqr = self.k[4] - yy1;
381                for i in 0..sqr.len() {
382                    stden += sqr.get(i).powi(2);
383                }
384
385                if stden > T::zero() {
386                    self.h_lamb = self.h * (stdnum / stden).sqrt();
387                }
388                if self.h_lamb > T::from_f64(6.1).unwrap() {
389                    self.non_stiff_counter = 0;
390                    self.stiffness_counter += 1;
391                    if self.stiffness_counter == 15 {
392                        // Early Exit Stiffness Detected
393                        self.status = Status::Error(Error::Stiffness {
394                            t: self.t,
395                            y: self.y,
396                        });
397                        return Err(Error::Stiffness {
398                            t: self.t,
399                            y: self.y,
400                        });
401                    }
402                } else {
403                    self.non_stiff_counter += 1;
404                    if self.non_stiff_counter == 6 {
405                        self.stiffness_counter = 0;
406                    }
407                }
408            }
409
410            // Preperation for dense / continuous output
411            self.cont[0] = self.y;
412            let ydiff = self.k[4] - self.y;
413            self.cont[1] = ydiff;
414            let bspl = self.k[0] * self.h - ydiff;
415            self.cont[2] = bspl;
416            self.cont[3] = ydiff - self.k[3] * self.h - bspl;
417
418            self.cont[4] = self.k[0] * self.b_dense[0][0]
419                + self.k[5] * self.b_dense[0][5]
420                + self.k[6] * self.b_dense[0][6]
421                + self.k[7] * self.b_dense[0][7]
422                + self.k[8] * self.b_dense[0][8]
423                + self.k[9] * self.b_dense[0][9]
424                + self.k[1] * self.b_dense[0][10]
425                + self.k[2] * self.b_dense[0][11];
426
427            self.cont[5] = self.k[0] * self.b_dense[1][0]
428                + self.k[5] * self.b_dense[1][5]
429                + self.k[6] * self.b_dense[1][6]
430                + self.k[7] * self.b_dense[1][7]
431                + self.k[8] * self.b_dense[1][8]
432                + self.k[9] * self.b_dense[1][9]
433                + self.k[1] * self.b_dense[1][10]
434                + self.k[2] * self.b_dense[1][11];
435
436            self.cont[6] = self.k[0] * self.b_dense[2][0]
437                + self.k[5] * self.b_dense[2][5]
438                + self.k[6] * self.b_dense[2][6]
439                + self.k[7] * self.b_dense[2][7]
440                + self.k[8] * self.b_dense[2][8]
441                + self.k[9] * self.b_dense[2][9]
442                + self.k[1] * self.b_dense[2][10]
443                + self.k[2] * self.b_dense[2][11];
444
445            self.cont[7] = self.k[0] * self.b_dense[3][0]
446                + self.k[5] * self.b_dense[3][5]
447                + self.k[6] * self.b_dense[3][6]
448                + self.k[7] * self.b_dense[3][7]
449                + self.k[8] * self.b_dense[3][8]
450                + self.k[9] * self.b_dense[3][9]
451                + self.k[1] * self.b_dense[3][10]
452                + self.k[2] * self.b_dense[3][11];
453
454            ode.diff(
455                self.t + self.c_dense[0] * self.h,
456                &(self.y
457                    + (self.k[0] * self.a_dense[0][0]
458                        + self.k[6] * self.a_dense[0][6]
459                        + self.k[7] * self.a_dense[0][7]
460                        + self.k[8] * self.a_dense[0][8]
461                        + self.k[9] * self.a_dense[0][9]
462                        + self.k[1] * self.a_dense[0][10]
463                        + self.k[2] * self.a_dense[0][11]
464                        + self.k[3] * self.a_dense[0][12])
465                        * self.h),
466                &mut self.k[9],
467            );
468
469            ode.diff(
470                self.t + self.c_dense[1] * self.h,
471                &(self.y
472                    + (self.k[0] * self.a_dense[1][0]
473                        + self.k[5] * self.a_dense[1][5]
474                        + self.k[6] * self.a_dense[1][6]
475                        + self.k[7] * self.a_dense[1][7]
476                        + self.k[1] * self.a_dense[1][10]
477                        + self.k[2] * self.a_dense[1][11]
478                        + self.k[3] * self.a_dense[1][12]
479                        + self.k[9] * self.a_dense[1][13])
480                        * self.h),
481                &mut self.k[1],
482            );
483
484            ode.diff(
485                self.t + self.c_dense[2] * self.h,
486                &(self.y
487                    + (self.k[0] * self.a_dense[2][0]
488                        + self.k[5] * self.a_dense[2][5]
489                        + self.k[6] * self.a_dense[2][6]
490                        + self.k[7] * self.a_dense[2][7]
491                        + self.k[8] * self.a_dense[2][8]
492                        + self.k[3] * self.a_dense[2][12]
493                        + self.k[9] * self.a_dense[2][13]
494                        + self.k[1] * self.a_dense[2][14])
495                        * self.h),
496                &mut self.k[2],
497            );
498
499            // Log evals
500            evals.fcn += 3; // Increment function evaluations for the dense output calculations
501
502            // Final preparation - add contributions from the extra stages and scale
503            self.cont[4] = (self.cont[4]
504                + self.k[3] * self.b_dense[0][12]
505                + self.k[9] * self.b_dense[0][13]
506                + self.k[1] * self.b_dense[0][14]
507                + self.k[2] * self.b_dense[0][15])
508                * self.h;
509
510            self.cont[5] = (self.cont[5]
511                + self.k[3] * self.b_dense[1][12]
512                + self.k[9] * self.b_dense[1][13]
513                + self.k[1] * self.b_dense[1][14]
514                + self.k[2] * self.b_dense[1][15])
515                * self.h;
516
517            self.cont[6] = (self.cont[6]
518                + self.k[3] * self.b_dense[2][12]
519                + self.k[9] * self.b_dense[2][13]
520                + self.k[1] * self.b_dense[2][14]
521                + self.k[2] * self.b_dense[2][15])
522                * self.h;
523
524            self.cont[7] = (self.cont[7]
525                + self.k[3] * self.b_dense[3][12]
526                + self.k[9] * self.b_dense[3][13]
527                + self.k[1] * self.b_dense[3][14]
528                + self.k[2] * self.b_dense[3][15])
529                * self.h;
530
531            // For Interpolation
532            self.y_old = self.y;
533            self.t_old = self.t;
534            self.h_old = self.h;
535
536            // Update State
537            self.k[0] = self.k[3];
538            self.y = self.k[4];
539            self.t = t_new;
540
541            // Check if previous step rejected
542            if let Status::RejectedStep = self.status {
543                h_new = self.h.min(h_new);
544                self.status = Status::Solving;
545            }
546        } else {
547            // Step Rejected
548            h_new = self.h / self.facc1.min(self.fac11 / self.safe);
549            self.status = Status::RejectedStep;
550        }
551
552        // Step Complete
553        self.h = constrain_step_size(h_new, self.h_min, self.h_max);
554        Ok(evals)
555    }
556
557    fn t(&self) -> T {
558        self.t
559    }
560
561    fn y(&self) -> &V {
562        &self.y
563    }
564
565    fn t_prev(&self) -> T {
566        self.t_old
567    }
568
569    fn y_prev(&self) -> &V {
570        &self.y_old
571    }
572
573    fn h(&self) -> T {
574        self.h
575    }
576
577    fn set_h(&mut self, h: T) {
578        self.h = h;
579    }
580
581    fn status(&self) -> &Status<T, V, D> {
582        &self.status
583    }
584
585    fn set_status(&mut self, status: Status<T, V, D>) {
586        self.status = status;
587    }
588}
589
590impl<T: Real, V: State<T>, D: CallBackData> Interpolation<T, V> for DOP853<T, V, D> {
591    fn interpolate(&mut self, t_interp: T) -> Result<V, Error<T, V>> {
592        // Check if interpolation is out of bounds
593        if t_interp < self.t_old || t_interp > self.t {
594            return Err(Error::OutOfBounds {
595                t_interp,
596                t_prev: self.t_old,
597                t_curr: self.t,
598            });
599        }
600
601        // Evaluate the interpolation polynomial at the requested time
602        let s = (t_interp - self.t_old) / self.h_old;
603        let s1 = T::one() - s;
604
605        // Compute the interpolated value using nested polynomial evaluation
606        let conpar = self.cont[4] + (self.cont[5] + (self.cont[6] + self.cont[7] * s) * s1) * s;
607
608        let y_interp = self.cont[0]
609            + (self.cont[1] + (self.cont[2] + (self.cont[3] + conpar * s1) * s) * s1) * s;
610
611        Ok(y_interp)
612    }
613}
614
615impl<T: Real, V: State<T>, D: CallBackData> DOP853<T, V, D> {
616    /// Creates a new DOP853 ODENumericalMethod.
617    ///
618    /// # Returns
619    /// * `system` - Function that defines the ordinary differential equation dy/dt = f(t, y).
620    /// # Returns
621    /// * DOP853 Struct ready to go for solving.
622    ///  
623    pub fn new() -> Self {
624        DOP853 {
625            ..Default::default()
626        }
627    }
628
629    // Builder Functions
630    pub fn rtol(mut self, rtol: T) -> Self {
631        self.rtol = rtol;
632        self
633    }
634
635    pub fn atol(mut self, atol: T) -> Self {
636        self.atol = atol;
637        self
638    }
639
640    pub fn h0(mut self, h0: T) -> Self {
641        self.h0 = h0;
642        self
643    }
644
645    pub fn h_max(mut self, h_max: T) -> Self {
646        self.h_max = h_max;
647        self
648    }
649
650    pub fn h_min(mut self, h_min: T) -> Self {
651        self.h_min = h_min;
652        self
653    }
654
655    pub fn max_steps(mut self, max_steps: usize) -> Self {
656        self.max_steps = max_steps;
657        self
658    }
659
660    pub fn n_stiff(mut self, n_stiff: usize) -> Self {
661        self.n_stiff = n_stiff;
662        self
663    }
664
665    pub fn safe(mut self, safe: T) -> Self {
666        self.safe = safe;
667        self
668    }
669
670    pub fn beta(mut self, beta: T) -> Self {
671        self.beta = beta;
672        self
673    }
674
675    pub fn fac1(mut self, fac1: T) -> Self {
676        self.fac1 = fac1;
677        self
678    }
679
680    pub fn fac2(mut self, fac2: T) -> Self {
681        self.fac2 = fac2;
682        self
683    }
684
685    pub fn expo1(mut self, expo1: T) -> Self {
686        self.expo1 = expo1;
687        self
688    }
689
690    pub fn facc1(mut self, facc1: T) -> Self {
691        self.facc1 = facc1;
692        self
693    }
694
695    pub fn facc2(mut self, facc2: T) -> Self {
696        self.facc2 = facc2;
697        self
698    }
699}
700
701impl<T: Real, V: State<T>, D: CallBackData> Default for DOP853<T, V, D> {
702    fn default() -> Self {
703        // Convert coefficient arrays from f64 to type T
704        let a = DOP853_A.map(|row| row.map(|x| T::from_f64(x).unwrap()));
705        let b = DOP853_B.map(|x| T::from_f64(x).unwrap());
706        let c = DOP853_C.map(|x| T::from_f64(x).unwrap());
707        let er = DOP853_ER.map(|x| T::from_f64(x).unwrap());
708        let bhh = DOP853_BHH.map(|x| T::from_f64(x).unwrap());
709
710        let a_dense = DOP853_A_DENSE.map(|row| row.map(|x| T::from_f64(x).unwrap()));
711        let c_dense = DOP853_C_DENSE.map(|x| T::from_f64(x).unwrap());
712        let b_dense = DOP853_DENSE.map(|row| row.map(|x| T::from_f64(x).unwrap()));
713
714        // Create arrays of zeros for k and cont matrices
715        let k_zeros = [V::zeros(); 12];
716        let cont_zeros = [V::zeros(); 8];
717
718        DOP853 {
719            // State Variables
720            t: T::zero(),
721            y: V::zeros(),
722            h: T::zero(),
723
724            // Settings
725            h0: T::zero(),
726            rtol: T::from_f64(1e-3).unwrap(),
727            atol: T::from_f64(1e-6).unwrap(),
728            h_max: T::infinity(),
729            h_min: T::zero(),
730            max_steps: 1_000_000,
731            n_stiff: 100,
732            safe: T::from_f64(0.9).unwrap(),
733            fac1: T::from_f64(0.33).unwrap(),
734            fac2: T::from_f64(6.0).unwrap(),
735            beta: T::from_f64(0.0).unwrap(),
736            expo1: T::from_f64(1.0 / 8.0).unwrap(),
737            facc1: T::from_f64(1.0 / 0.33).unwrap(),
738            facc2: T::from_f64(1.0 / 6.0).unwrap(),
739            facold: T::from_f64(1.0e-4).unwrap(),
740            fac11: T::zero(),
741            fac: T::zero(),
742
743            // Butcher Tableau Coefficients
744            a,
745            b,
746            c,
747            er,
748            bhh,
749            a_dense,
750            c_dense,
751            b_dense,
752
753            // Status and Counters
754            status: Status::Uninitialized,
755            h_lamb: T::zero(),
756            non_stiff_counter: 0,
757            stiffness_counter: 0,
758            steps: 0,
759            n_accepted: 0,
760
761            // Coefficents and temporary storage
762            k: k_zeros,
763            y_old: V::zeros(),
764            t_old: T::zero(),
765            h_old: T::zero(),
766            cont: cont_zeros,
767        }
768    }
769}
770
771// DOP853 Butcher Tableau
772
773// 12 Stage Core
774
775// A matrix (12x12, lower triangular)
776const DOP853_A: [[f64; 12]; 12] = [
777    [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
778    [
779        5.260_015_195_876_773E-2,
780        0.0,
781        0.0,
782        0.0,
783        0.0,
784        0.0,
785        0.0,
786        0.0,
787        0.0,
788        0.0,
789        0.0,
790        0.0,
791    ],
792    [
793        1.972_505_698_453_79E-2,
794        5.917_517_095_361_37E-2,
795        0.0,
796        0.0,
797        0.0,
798        0.0,
799        0.0,
800        0.0,
801        0.0,
802        0.0,
803        0.0,
804        0.0,
805    ],
806    [
807        2.958_758_547_680_685E-2,
808        0.0,
809        8.876_275_643_042_054E-2,
810        0.0,
811        0.0,
812        0.0,
813        0.0,
814        0.0,
815        0.0,
816        0.0,
817        0.0,
818        0.0,
819    ],
820    [
821        2.413_651_341_592_667E-1,
822        0.0,
823        -8.845_494_793_282_861E-1,
824        9.248_340_032_617_92E-1,
825        0.0,
826        0.0,
827        0.0,
828        0.0,
829        0.0,
830        0.0,
831        0.0,
832        0.0,
833    ],
834    [
835        3.703_703_703_703_703_5E-2,
836        0.0,
837        0.0,
838        1.708_286_087_294_738_6E-1,
839        1.254_676_875_668_224_2E-1,
840        0.0,
841        0.0,
842        0.0,
843        0.0,
844        0.0,
845        0.0,
846        0.0,
847    ],
848    [
849        3.7109375E-2,
850        0.0,
851        0.0,
852        1.702_522_110_195_440_5E-1,
853        6.021_653_898_045_596E-2,
854        -1.7578125E-2,
855        0.0,
856        0.0,
857        0.0,
858        0.0,
859        0.0,
860        0.0,
861    ],
862    [
863        3.709_200_011_850_479E-2,
864        0.0,
865        0.0,
866        1.703_839_257_122_399_8E-1,
867        1.072_620_304_463_732_8E-1,
868        -1.531_943_774_862_440_2E-2,
869        8.273_789_163_814_023E-3,
870        0.0,
871        0.0,
872        0.0,
873        0.0,
874        0.0,
875    ],
876    [
877        6.241_109_587_160_757E-1,
878        0.0,
879        0.0,
880        -3.360_892_629_446_941_4,
881        -8.682_193_468_417_26E-1,
882        2.759_209_969_944_671E1,
883        2.015_406_755_047_789_4E1,
884        -4.348_988_418_106_996E1,
885        0.0,
886        0.0,
887        0.0,
888        0.0,
889    ],
890    [
891        4.776_625_364_382_643_4E-1,
892        0.0,
893        0.0,
894        -2.488_114_619_971_667_7,
895        -5.902_908_268_368_43E-1,
896        2.123_005_144_818_119_3E1,
897        1.527_923_363_288_242_3E1,
898        -3.328_821_096_898_486E1,
899        -2.033_120_170_850_862_7E-2,
900        0.0,
901        0.0,
902        0.0,
903    ],
904    [
905        -9.371_424_300_859_873E-1,
906        0.0,
907        0.0,
908        5.186_372_428_844_064,
909        1.091_437_348_996_729_5,
910        -8.149_787_010_746_927,
911        -1.852_006_565_999_696E1,
912        2.273_948_709_935_050_5E1,
913        2.493_605_552_679_652_3,
914        -3.046_764_471_898_219_6,
915        0.0,
916        0.0,
917    ],
918    [
919        2.273_310_147_516_538,
920        0.0,
921        0.0,
922        -1.053_449_546_673_725E1,
923        -2.000_872_058_224_862_5,
924        -1.795_893_186_311_88E1,
925        2.794_888_452_941_996E1,
926        -2.858_998_277_135_023_5,
927        -8.872_856_933_530_63,
928        1.236_056_717_579_430_3E1,
929        6.433_927_460_157_636E-1,
930        0.0,
931    ],
932];
933
934// C coefficients (nodes)
935const DOP853_C: [f64; 12] = [
936    0.0,                      // C1 (not given in constants, but must be 0)
937    5.260_015_195_876_773E-2, // C2
938    7.890_022_793_815_16E-2,  // C3
939    1.183_503_419_072_274E-1, // C4
940    2.816_496_580_927_726E-1, // C5
941    3.333_333_333_333_333E-1, // C6
942    0.25E+00,                 // C7
943    3.076_923_076_923_077E-1, // C8
944    6.512_820_512_820_513E-1, // C9
945    0.6E+00,                  // C10
946    8.571_428_571_428_571E-1, // C11
947    1.0,                      // C12 (final point, not explicitly given)
948];
949
950// B coefficients (weights for main method)
951const DOP853_B: [f64; 12] = [
952    5.429_373_411_656_876_5E-2, // B1
953    0.0,                        // B2-B5 are zero (not explicitly listed)
954    0.0,
955    0.0,
956    0.0,
957    4.450_312_892_752_409,      // B6
958    1.891_517_899_314_500_3,    // B7
959    -5.801_203_960_010_585,     // B8
960    3.111_643_669_578_199E-1,   // B9
961    -1.521_609_496_625_161E-1,  // B10
962    2.013_654_008_040_303_4E-1, // B11
963    4.471_061_572_777_259E-2,   // B12
964];
965
966// Error estimation coefficients
967
968// Error estimation coefficients (constructed from ER values)
969const DOP853_ER: [f64; 12] = [
970    1.312_004_499_419_488E-2, // ER1
971    0.0,                      // ER2-ER5 are zero
972    0.0,
973    0.0,
974    0.0,
975    -1.225_156_446_376_204_4,    // ER6
976    -4.957_589_496_572_502E-1,   // ER7
977    1.664_377_182_454_986_4,     // ER8
978    -3.503_288_487_499_736_6E-1, // ER9
979    3.341_791_187_130_175E-1,    // ER10
980    8.192_320_648_511_571E-2,    // ER11
981    -2.235_530_786_388_629_4E-2, // ER12
982];
983
984const DOP853_BHH: [f64; 3] = [
985    2.440_944_881_889_764E-1,   // BHH1
986    7.338_466_882_816_118E-1,   // BHH2
987    2.205_882_352_941_176_6E-2, // BHH3
988];
989
990// Dense output Coefficients
991
992// Dense output A coefficients (for the 3 extra stages used in interpolation)
993const DOP853_A_DENSE: [[f64; 16]; 3] = [
994    // Stage 14 coefficients (C14 = 0.1)
995    [
996        5.616_750_228_304_795_4E-2, // A141
997        0.0,
998        0.0,
999        0.0,
1000        0.0,
1001        0.0,                         // A142-A146 (zero)
1002        2.535_002_102_166_248_3E-1,  // A147
1003        -2.462_390_374_708_025E-1,   // A148
1004        -1.241_914_232_638_163_7E-1, // A149
1005        1.532_917_982_787_656_8E-1,  // A1410
1006        8.201_052_295_634_69E-3,     // A1411
1007        7.567_897_660_545_699E-3,    // A1412
1008        -8.298E-3,                   // A1413
1009        0.0,
1010        0.0,
1011        0.0, // A1414-A1416 (zero/not used)
1012    ],
1013    // Stage 15 coefficients (C15 = 0.2)
1014    [
1015        3.183_464_816_350_214E-2, // A151
1016        0.0,
1017        0.0,
1018        0.0,
1019        0.0,                        // A152-A155 (zero)
1020        2.830_090_967_236_677_6E-2, // A156
1021        5.354_198_830_743_856_6E-2, // A157
1022        -5.492_374_857_139_099E-2,  // A158
1023        0.0,
1024        0.0,                         // A159-A1510 (zero)
1025        -1.083_473_286_972_493_2E-4, // A1511
1026        3.825_710_908_356_584E-4,    // A1512
1027        -3.404_650_086_874_045_6E-4, // A1513
1028        1.413_124_436_746_325E-1,    // A1514
1029        0.0,
1030        0.0, // A1515-A1516 (zero/not used)
1031    ],
1032    // Stage 16 coefficients (C16 = 0.777...)
1033    [
1034        -4.288_963_015_837_919_4E-1, // A161
1035        0.0,
1036        0.0,
1037        0.0,
1038        0.0,                      // A162-A165 (zero)
1039        -4.697_621_415_361_164,   // A166
1040        7.683_421_196_062_599,    // A167
1041        4.068_989_818_397_11,     // A168
1042        3.567_271_874_552_811E-1, // A169
1043        0.0,
1044        0.0,
1045        0.0,                         // A1610-A1612 (zero)
1046        -1.399_024_165_159_014_5E-3, // A1613
1047        2.947_514_789_152_772_4,     // A1614
1048        -9.150_958_472_179_87,       // A1615
1049        0.0,                         // A1616 (not used)
1050    ],
1051];
1052
1053const DOP853_C_DENSE: [f64; 3] = [
1054    0.1E+00,                  // C14
1055    0.2E+00,                  // C15
1056    7.777_777_777_777_778E-1, // C16
1057];
1058
1059// Dense output coefficients for stage 4
1060const DOP853_D4: [f64; 16] = [
1061    -8.428_938_276_109_013, // D41
1062    0.0,
1063    0.0,
1064    0.0,
1065    0.0,                       // D42-D45 are zero
1066    5.667_149_535_193_777E-1,  // D46
1067    -3.068_949_945_949_891_7,  // D47
1068    2.384_667_656_512_07,      // D48
1069    2.117_034_582_445_028,     // D49
1070    -8.713_915_837_779_73E-1,  // D410
1071    2.240_437_430_260_788_3,   // D411
1072    6.315_787_787_694_688E-1,  // D412
1073    -8.899_033_645_133_331E-2, // D413
1074    1.814_850_552_085_472_7E1, // D414
1075    -9.194_632_392_478_356,    // D415
1076    -4.436_036_387_594_894,    // D416
1077];
1078
1079// Dense output coefficients for stages 5, 6, and 7 follow same pattern
1080const DOP853_D5: [f64; 16] = [
1081    1.042_750_864_257_913_4E1, // D51
1082    0.0,
1083    0.0,
1084    0.0,
1085    0.0,                        // D52-D55 are zero
1086    2.422_834_917_752_581_7E2,  // D56
1087    1.652_004_517_172_702_8E2,  // D57
1088    -3.745_467_547_226_902E2,   // D58
1089    -2.211_366_685_312_530_6E1, // D59
1090    7.733_432_668_472_264,      // D510
1091    -3.067_408_473_108_939_8E1, // D511
1092    -9.332_130_526_430_229,     // D512
1093    1.569_723_812_177_084_5E1,  // D513
1094    -3.113_940_321_956_517_8E1, // D514
1095    -9.352_924_358_844_48,      // D515
1096    3.581_684_148_639_408E1,    // D516
1097];
1098
1099const DOP853_D6: [f64; 16] = [
1100    1.998_505_324_200_243_3E1, // D61
1101    0.0,
1102    0.0,
1103    0.0,
1104    0.0,                        // D62-D65 are zero
1105    -3.870_373_087_493_518E2,   // D66
1106    -1.891_781_381_951_675_8E2, // D67
1107    5.278_081_592_054_236E2,    // D68
1108    -1.157_390_253_995_963E1,   // D69
1109    6.881_232_694_696_3,        // D610
1110    -1.000_605_096_691_083_8,   // D611
1111    7.777_137_798_053_443E-1,   // D612
1112    -2.778_205_752_353_508,     // D613
1113    -6.019_669_523_126_412E1,   // D614
1114    8.432_040_550_667_716E1,    // D615
1115    1.199_229_113_618_279E1,    // D616
1116];
1117
1118const DOP853_D7: [f64; 16] = [
1119    -2.569_393_346_270_375E1, // D71
1120    0.0,
1121    0.0,
1122    0.0,
1123    0.0,                        // D72-D75 are zero
1124    -1.541_897_486_902_364_3E2, // D76
1125    -2.315_293_791_760_455E2,   // D77
1126    3.576_391_179_106_141E2,    // D78
1127    9.340_532_418_362_432E1,    // D79
1128    -3.745_832_313_645_163E1,   // D710
1129    1.040_996_495_089_623E2,    // D711
1130    2.984_029_342_666_05E1,     // D712
1131    -4.353_345_659_001_114E1,   // D713
1132    9.632_455_395_918_828E1,    // D714
1133    -3.917_726_167_561_544E1,   // D715
1134    -1.497_268_362_579_856_4E2, // D716
1135];
1136
1137// Dense output coefficients as a 3D array [stage][coefficient_index]
1138const DOP853_DENSE: [[f64; 16]; 4] = [DOP853_D4, DOP853_D5, DOP853_D6, DOP853_D7];
differential_equations/ode/methods/runge_kutta/explicit/dop853.rs

differential_equations/ode/methods/runge_kutta/explicit/
dop853.rs
