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
//use numerical::utils::*;
//use structure::dual::*;
///// Backward Differentiation Formula
///// Author: Axect
//use structure::matrix::*;
//use structure::vector::*;
//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))
//}