1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126

//! Butcher tableau
use crate::analysis::differential_equation::ordinary::ExplicitODE;
use crate::algebra::linear::Vector;
use crate::algebra::abstr::Real;

use std::fmt::Debug;

#[cfg(feature = "serde")]
use serde::{Deserialize, Serialize};


/// ```math
/// y_{n+1} = y_n + h \sum_{i = 1}^{s}b_i k_i
/// ```
/// ```math
/// k_i = f(t_n + c_ih, y_n + h \sum_{j=1}^{s}a_{ij}k_j), \quad i = 1, \dots , s
/// ```
#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
#[derive(Clone, Debug)]
pub struct ButcherFixedStepSize<T>
{
    a: Vec<T>,
    b: Vec<T>,
    c: Vec<T>,
}

impl<T> ButcherFixedStepSize<T>
{
    pub fn new(a: Vec<T>, b: Vec<T>, c: Vec<T>) -> ButcherFixedStepSize<T>
    {
        ButcherFixedStepSize
        {
            a,
            b,
            c
        }
    }
}

impl<T> ButcherFixedStepSize<T>
    where T: Real
{
    pub fn do_step<F>(self: &Self, prob: &F, t_n: &T, x_n: &Vector<T>, h: &T) -> Vector<T>
        where F: ExplicitODE<T>,
    {
        let mut k: Vec<Vector<T>> = Vec::with_capacity(self.b.len());
        let (rows, _columns): (usize, usize) = x_n.dim();

        k.push(prob.func(t_n, x_n));

        for j in 1..self.b.len()
        {
            let i_b = (j - 1) * j / 2;
            let i_e = i_b + j;

            let sum= self.a[i_b..i_e].iter().zip(k.iter()).map(|(a_jl, k_l)| k_l * a_jl).fold(Vector::zero(rows),  | a, b| a + b);

            k.push(prob.func(&(*t_n + self.c[j - 1] * *h), &(x_n + &(&sum * h))));
        }

        let sum: Vector<T> = self.b.iter().zip(k.iter()).map(|(b, k_j)|  k_j * b).fold(Vector::zero(rows), |a, b| a + b);
        x_n + &(&sum * h)
    }
}

#[cfg_attr(feature = "serde", derive(Serialize, Deserialize))]
#[derive(Clone, Debug)]
pub struct ButcherAdaptiveStepSize<T>
{
    a: Vec<T>,
    b: Vec<T>,
    b_s: Vec<T>,
    c: Vec<T>,
}

///
/// ```math
/// y_{n+1}^* = y_n 1 h \sum_{i=1}^{s} b_{i}^{s}$k_{i}
/// ```
///
impl<T: Debug> ButcherAdaptiveStepSize<T>
{
    pub fn new(a: Vec<T>, b: Vec<T>, b_s: Vec<T>, c: Vec<T>) -> ButcherAdaptiveStepSize<T>
    {
        ButcherAdaptiveStepSize
        {
            a,
            b,
            b_s,
            c
        }
    }
}

impl<T> ButcherAdaptiveStepSize<T>
    where T: Real + Debug
{
    pub fn do_step<F>(self: &Self, prob: &F, t_n: &T, x_n: &Vector<T>, h: &T) -> (Vector<T>, Vector<T>)
        where F: ExplicitODE<T>,
    {
        let mut k: Vec<Vector<T>> = Vec::with_capacity(self.b.len());
        let (rows, _columns): (usize, usize) = x_n.dim();

        k.push(prob.func(t_n, x_n));

        for j in 1..self.b.len()
        {
            let i_b = (j - 1) * j / 2;
            let i_e = i_b + j;

            let sum= self.a[i_b..i_e].iter().zip(k.iter()).map(|(a_jl, k_l)| k_l * a_jl).fold(Vector::zero(rows),  | a, b| a + b);

            let k_i = prob.func(&(*t_n + self.c[j - 1] * *h), &(x_n + &(&sum * h))) ;

            k.push(k_i);
        }

        let sum: Vector<T> = self.b.iter().zip(k.iter()).map(|(b, k_j)| k_j * b).fold(Vector::zero(rows), |a, b| a + b);
        let x_n_1 = x_n + &(&sum * h);

        let sum_s: Vector<T> = self.b_s.iter().zip(k.iter()).map(|(b, k_j)|  k_j * b).fold(Vector::zero(rows), |a, b| a + b);
        let x_s_n_1 = x_n + &(&sum_s * h);

        (x_n_1, x_s_n_1)
    }
}