scirs2-integrate 0.4.1

Numerical integration module for SciRS2 (scirs2-integrate)
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

use scirs2_core::ndarray::{array, ArrayView1};
use scirs2_integrate::ode::{solve_ivp, ODEMethod, ODEOptions};

#[allow(dead_code)]
fn main() {
    println!("LSODA Solver Debug Example");
    println!("-------------------------");
    println!("This example is for debugging the LSODA implementation");
    println!("We'll try different step sizes and parameters to understand the failures");

    // Use a very simple test problem: y' = -y (exponential decay)
    let decay_system = |_t: f64, y: ArrayView1<f64>| array![-y[0]];

    // Try with different initial step sizes and min step thresholds
    let step_sizes = [0.1, 0.01, 0.001];
    let min_steps = [1e-3, 1e-4, 1e-5, 1e-6];

    for &initial_step in &step_sizes {
        for &min_step in &min_steps {
            println!("\nTrying with initial_step={initial_step}, min_step={min_step}");

            let result = solve_ivp(
                decay_system,
                [0.0, 1.0], // Short time interval for testing
                array![1.0],
                Some(ODEOptions {
                    method: ODEMethod::LSODA,
                    rtol: 1e-4,
                    atol: 1e-6,
                    max_steps: 100,
                    h0: Some(initial_step),
                    min_step: Some(min_step),
                    ..Default::default()
                }),
            );

            match result {
                Ok(res) => {
                    println!(
                        "  Success! Steps: {}, Function evaluations: {}",
                        res.n_steps, res.n_eval
                    );
                    println!(
                        "  Final value: {}, Exact: {}",
                        res.y.last().expect("Operation failed")[0],
                        (-1.0f64).exp()
                    );
                    println!(
                        "  Error: {:.2e}",
                        (res.y.last().expect("Operation failed")[0] - (-1.0f64).exp()).abs()
                    );

                    if let Some(msg) = res.message {
                        println!("  Message: {msg}");
                    }
                }
                Err(e) => {
                    println!("  Failed: {e}");
                }
            }
        }
    }

    // Try a stiff problem with various parameters
    println!("\n\nTesting with stiff Van der Pol oscillator");
    let mu = 1000.0; // Very stiff
    let van_der_pol =
        |_t: f64, y: ArrayView1<f64>| array![y[1], mu * (1.0 - y[0].powi(2)) * y[1] - y[0]];

    // Try with different parameters
    let configs = [(0.01, 1e-4, 500), (0.001, 1e-5, 1000), (0.0001, 1e-6, 2000)];

    for &(initial_step, min_step, max_steps) in &configs {
        println!(
            "\nTrying Van der Pol with: initial_step={initial_step}, min_step={min_step}, max_steps={max_steps}"
        );

        let result = solve_ivp(
            van_der_pol,
            [0.0, 2.0], // Short interval for testing
            array![2.0, 0.0],
            Some(ODEOptions {
                method: ODEMethod::LSODA,
                rtol: 1e-4,
                atol: 1e-6,
                max_steps,
                h0: Some(initial_step),
                min_step: Some(min_step),
                ..Default::default()
            }),
        );

        match result {
            Ok(res) => {
                println!(
                    "  Success! Steps: {}, Function evaluations: {}",
                    res.n_steps, res.n_eval
                );
                println!(
                    "  Final state: [{:.4}, {:.4}]",
                    res.y.last().expect("Operation failed")[0],
                    res.y.last().expect("Operation failed")[1]
                );

                if let Some(msg) = res.message {
                    println!("  Message: {msg}");
                }
            }
            Err(e) => {
                println!("  Failed: {e}");
            }
        }
    }
}