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

/// Backward Differentiation Formula
/// Author: Axect

use structure::matrix::*;
use structure::vector::*;
use structure::dual::*;
use numerical::utils::*;
use util::non_macro::*;


// =============================================================================
// BDF1
// =============================================================================

/// BDF1 (Backward Euler Method)
///
/// t_init: initial param
/// ys_init: initial values
/// h: Step size
/// rtol: Relative tolerance (1e-15 is recommended)
pub fn bdf1<F>(t_init: f64, ys_init: Vec<f64>, f: F, h: f64, rtol: f64, num: usize) -> Matrix
    where F: Fn(Dual, Vec<Dual>) -> Vec<Dual> + Copy
{
    let mut ys = ys_init.clone();
    let mut t = t_init.clone();
    let mut records = zeros(num + 1, ys_init.len() + 1);
    for i in 0 .. (num+1) {
        records.subs_row(i, concat(vec![t], ys.clone()));
        t += h;
        ys = one_step_bdf1(t, ys.clone(), f, h, rtol);
    }

    records
}


/// One step for BDF1 (Backward Euler)
pub fn one_step_bdf1<F>(t: f64, ys: Vec<f64>, f: F, h: f64, rtol: f64) -> Vec<f64>
    where F: Fn(Dual, Vec<Dual>) -> Vec<Dual> + Copy
{
    // new t = t + h
    let new_t = t + h;

    // One step forward euler for initial guess
    // y_{n+1}^{(0)} = y_n + hf(t_n, y_n)
    let new_ys_0 = ys.add(
        &(f(dual(t, 0.), ys.conv_dual()).values().fmap(|y| h*y))
    );

    let mut iter_prev_ys = new_ys_0;
    let mut iter_next_ys = non_auto_update(new_t, ys.clone(), iter_prev_ys.clone(), f, h);
    let mut err = iter_next_ys.sub(&iter_prev_ys).norm();
    let mut max_iter = 0;

    while err >= rtol && max_iter <= 10 {
        iter_prev_ys = iter_next_ys.clone();
        iter_next_ys = non_auto_update(new_t, ys.clone(), iter_prev_ys.clone(), f, h);
        err = iter_next_ys.sub(&iter_prev_ys).norm();
        max_iter += 1;
    }

    iter_next_ys
}

#[allow(non_snake_case)]
fn non_auto_update<F>(t: f64, yn: Vec<f64>, ys: Vec<f64>, f: F, h: f64) -> Vec<f64>
    where F: Fn(Dual, Vec<Dual>) -> Vec<Dual> + Copy
{
    let n = ys.len();
    let f_vec = |_xs: Vec<Dual>| f(dual(t, 0.), ys.conv_dual()); // n x 1

    // Df_y_{ij} = ∂f_i/∂y_j where y = f(t, y)
    let Df = jacobian(ys.clone(), f_vec); // n x n
    let I = eye(n);
    let DF = I - h * Df;
    let f_xs = f(dual(t, 0.), ys.conv_dual()).values();

    // F(y_{n+1}^i) = y_{n+1}^i - y_n - hf(t_n+h, y_{n+1}^i) 
    let F_ys = ys.sub(&yn).sub(&(f_xs.fmap(|x| h*x)));

    ys.sub(&(DF.inv().unwrap() % F_ys).col(0))
}