bacon-sci 0.16.2

Scientific computing in Rust
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

use crate::polynomial::Polynomial;
use crate::roots::{muller_polynomial, newton_polynomial};
use num_complex::Complex64;

#[test]
fn newton() {
    let mut poly = Polynomial::new();
    poly.set_coefficient(2, 1.0);
    poly.set_coefficient(0, -1.0);

    let solution = newton_polynomial(1.5, &poly, 0.0001, 1000).unwrap();
    assert!(approx_eq!(f64, solution, 1.0, epsilon = 0.0001));
    let solution = newton_polynomial(-1.5, &poly, 0.0001, 1000).unwrap();
    assert!(approx_eq!(f64, solution, -1.0, epsilon = 0.0001));

    poly.set_coefficient(2, 2.0);
    poly.set_coefficient(3, 5.0);
    poly.set_coefficient(1, -2.0);
    let solution = newton_polynomial(2.0, &poly, 0.0001, 1000).unwrap();
    assert!(approx_eq!(f64, solution, 0.66157, epsilon = 0.00001));
}

#[test]
fn muller() {
    let mut poly = polynomial![1.0, 0.0, -1.0];

    let solution = muller_polynomial((0.0, 0.5, 1.5), &poly, 0.0001, 1000).unwrap();
    assert!(approx_eq!(f64, solution.re, 1.0, epsilon = 0.0001));
    assert!(approx_eq!(f64, solution.im, 0.0, epsilon = 0.0001));
    let solution = muller_polynomial((0.0, -0.5, -1.5), &poly, 0.0001, 1000).unwrap();
    assert!(approx_eq!(f64, solution.re, -1.0, epsilon = 0.0001));
    assert!(approx_eq!(f64, solution.im, 0.0, epsilon = 0.0001));

    poly = polynomial![5.0, 2.0, -2.0, -1.0];
    let solution = muller_polynomial((0.0, 1.0, 2.0), &poly, 0.0001, 1000).unwrap();
    assert!(approx_eq!(f64, solution.re, 0.66157, epsilon = 0.00001));
    assert!(approx_eq!(f64, solution.im, 0.0, epsilon = 0.00001));
}

#[test]
fn polynomial_roots() {
    let poly = polynomial![1.0, -1.0];

    let roots = poly.roots(1e-10, 1000).unwrap();
    assert_eq!(roots.len(), 1);
    assert!(approx_eq!(f64, roots[0].re, 1.0, epsilon = 1e-10));
    assert!(approx_eq!(f64, roots[0].im, 0.0, epsilon = 1e-10));

    let poly = polynomial![0.0];
    let roots = poly.roots(1e-10, 1000).unwrap();
    assert_eq!(roots.len(), 1);
    assert!(approx_eq!(f64, roots[0].re, 0.0, epsilon = 1e-10));
    assert!(approx_eq!(f64, roots[0].im, 0.0, epsilon = 1e-10));

    let poly = polynomial![1.0];
    let roots = poly.roots(1e-10, 1000);
    assert!(roots.is_err(), "Should have not found a root");

    let poly = polynomial![1.0, 0.0, -1.0];
    let cmplx = poly.make_complex();
    let roots = poly.roots(1e-10, 1000).unwrap();
    assert_eq!(roots.len(), 2);
    for (ind, root) in roots.iter().enumerate() {
        let val = cmplx.evaluate(*root);
        assert!(approx_eq!(f64, val.re, 0.0, epsilon = 1e-10));
        assert!(approx_eq!(f64, val.im, 0.0, epsilon = 1e-10));
        for (j, r) in roots.iter().enumerate() {
            if j == ind {
                continue;
            }
            assert!(
                !(approx_eq!(f64, root.re, r.re, epsilon = 0.0001)
                    && approx_eq!(f64, root.im, r.im, epsilon = 0.0001))
            );
        }
    }

    let poly = polynomial![1.0, -2.0, -2.0, 1.0];
    let cmplx = poly.make_complex();
    let roots = poly.roots(1e-10, 1000).unwrap();
    assert_eq!(roots.len(), 3);
    for (ind, root) in roots.iter().enumerate() {
        let val = cmplx.evaluate(*root);
        assert!(approx_eq!(f64, val.re, 0.0, epsilon = 1e-10));
        assert!(approx_eq!(f64, val.im, 0.0, epsilon = 1e-10));
        for (j, r) in roots.iter().enumerate() {
            if j == ind {
                continue;
            }
            assert!(
                !(approx_eq!(f64, root.re, r.re, epsilon = 0.0001)
                    && approx_eq!(f64, root.im, r.im, epsilon = 0.0001))
            );
        }
    }

    let poly = polynomial![
        Complex64::new(5.0, 0.0),
        Complex64::new(-3.0, 0.0),
        Complex64::new(0.0, 0.0),
        Complex64::new(2.0, 0.0),
        Complex64::new(1.0, 0.0)
    ];
    let roots = poly.roots(1e-10, 1000).unwrap();
    assert_eq!(roots.len(), 4);
    for (ind, root) in roots.iter().enumerate() {
        let val = poly.evaluate(*root);
        assert!(approx_eq!(f64, val.re, 0.0, epsilon = 1e-10));
        assert!(approx_eq!(f64, val.im, 0.0, epsilon = 1e-10));
        for (j, r) in roots.iter().enumerate() {
            if j == ind {
                continue;
            }
            assert!(
                !(approx_eq!(f64, root.re, r.re, epsilon = 0.0001)
                    && approx_eq!(f64, root.im, r.im, epsilon = 0.0001))
            );
        }
    }
}