russell_ode 2.7.0

Solvers for ordinary differential equations and differential algebraic equations
Documentation 
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
127
128
129
130

use russell_lab::{Vector, approx_eq};
use russell_ode::{Method, OdeSolver, Params, Samples};

#[test]
fn test_bweuler_hairer_wanner_eq1() {
    // get ODE system
    let (system, x0, mut y0, mut args, y_fn_x) = Samples::hairer_wanner_eq1();
    let ndim = system.get_ndim();

    // final x
    let x1 = 1.5;

    // set configuration parameters
    let params = Params::new(Method::BwEuler);
    let mut solver = OdeSolver::new(params, system).unwrap();

    // solve the ODE system
    let h_equal = Some(1.875 / 50.0);
    solver.solve(&mut y0, x0, x1, h_equal, &mut args, None).unwrap();

    // get statistics
    let stat = solver.stats();

    // compare with a previous implementation
    approx_eq(y0[0], 0.09060476604187756, 1e-15);
    assert_eq!(stat.h_accepted, h_equal.unwrap());

    // compare with the analytical solution
    let mut y1_correct = Vector::new(ndim);
    y_fn_x(&mut y1_correct, x1, &mut args);
    approx_eq(y0[0], y1_correct[0], 5e-5);

    // print and check statistics
    println!("{}", stat);
    assert_eq!(stat.n_function, 80);
    assert_eq!(stat.n_jacobian, 40);
    assert_eq!(stat.n_factor, 40);
    assert_eq!(stat.n_lin_sol, 40);
    assert_eq!(stat.n_steps, 40);
    assert_eq!(stat.n_accepted, 40);
    assert_eq!(stat.n_rejected, 0);
    assert_eq!(stat.n_iterations, 2);
    assert_eq!(stat.n_iterations_max, 2);
}

#[test]
fn test_bweuler_hairer_wanner_eq1_num_jac() {
    // get ODE system
    let (system, x0, mut y0, mut args, y_fn_x) = Samples::hairer_wanner_eq1();
    let ndim = system.get_ndim();

    // final x
    let x1 = 1.5;

    // set configuration parameters
    let mut params = Params::new(Method::BwEuler);
    params.newton.use_numerical_jacobian = true;

    // solve the ODE system
    let mut solver = OdeSolver::new(params, system).unwrap();
    let h_equal = Some(1.875 / 50.0);
    solver.solve(&mut y0, x0, x1, h_equal, &mut args, None).unwrap();

    // get statistics
    let stat = solver.stats();

    // compare with a previous implementation
    approx_eq(y0[0], 0.09060476598021044, 1e-11);
    assert_eq!(stat.h_accepted, h_equal.unwrap());

    // compare with the analytical solution
    let mut y1_correct = Vector::new(ndim);
    y_fn_x(&mut y1_correct, x1, &mut args);
    approx_eq(y0[0], y1_correct[0], 5e-5);

    // print and check statistics
    println!("{}", stat);
    assert_eq!(stat.n_function, 120);
    assert_eq!(stat.n_jacobian, 40);
    assert_eq!(stat.n_factor, 40);
    assert_eq!(stat.n_lin_sol, 40);
    assert_eq!(stat.n_steps, 40);
    assert_eq!(stat.n_accepted, 40);
    assert_eq!(stat.n_rejected, 0);
    assert_eq!(stat.n_iterations, 2);
    assert_eq!(stat.n_iterations_max, 2);
}

#[test]
fn test_bweuler_hairer_wanner_eq1_modified_newton() {
    // get ODE system
    let (system, x0, mut y0, mut args, y_fn_x) = Samples::hairer_wanner_eq1();
    let ndim = system.get_ndim();

    // final x
    let x1 = 1.5;

    // set configuration parameters
    let mut params = Params::new(Method::BwEuler);
    params.bweuler.use_modified_newton = true;

    // solve the ODE system
    let mut solver = OdeSolver::new(params, system).unwrap();
    let h_equal = Some(1.875 / 50.0);
    solver.solve(&mut y0, x0, x1, h_equal, &mut args, None).unwrap();

    // get statistics
    let stat = solver.stats();

    // compare with a previous implementation
    approx_eq(y0[0], 0.09060476604187756, 1e-15);
    assert_eq!(stat.h_accepted, h_equal.unwrap());

    // compare with the analytical solution
    let mut y1_correct = Vector::new(ndim);
    y_fn_x(&mut y1_correct, x1, &mut args);
    approx_eq(y0[0], y1_correct[0], 5e-5);

    // print and check statistics
    println!("{}", stat);
    assert_eq!(stat.n_function, 80);
    assert_eq!(stat.n_jacobian, 1);
    assert_eq!(stat.n_factor, 1);
    assert_eq!(stat.n_lin_sol, 40);
    assert_eq!(stat.n_steps, 40);
    assert_eq!(stat.n_accepted, 40);
    assert_eq!(stat.n_rejected, 0);
    assert_eq!(stat.n_iterations, 2);
    assert_eq!(stat.n_iterations_max, 2);
}