numlib 0.1.0

A simple numerical library implementing common numberical algorithms 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
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166

// implementation of integrating techniques

/* Composite Trapezoidal Quadrature
 * Function: Computes an integral estimation using the composite trapezoidal rule
 * Parameters: 
 * f - the function,
 * n - the number of subintervals, 
 * a - the intial endpoint, 
 * b - the final endpoint
 * Returns:
 * The approximate solution to the integral
 */

pub fn composite_trapezoid(a: f64, b: f64, n: i64, f: fn(f64) -> f64) -> f64{
    let h : f64 = (b-a)/(n as f64).abs();
    let mut sum : f64 = 0.0;
    for i in 1..n {
        sum += f(a+(i as f64)*h);
    }
    return h*(sum + f(a)/(2 as f64) + f(b)/(2 as f64));
}

/* Simpsons Rule 
 * Function: Computes an integral estimation using simpsons rule. Preferable
 * to the trapezoidal rule in terms of error except on periodically smooth functions
 * Parameters: 
 * f - the function,
 * n - the number of subintervals,
 * a - the initial endpoint, 
 * b - the final endpoint
 * Returns:
 * The approximate solution to the integral
 *
 */

pub fn simpson_quadrature(a: f64, b: f64, n:i64, f: fn(f64) -> f64) -> f64{
    let h : f64 = (b-a)/(n as f64).abs();
    let mut sum :f64 = 0.0;
    for i in 1..n{
        if i % 2 == 1 {
            sum += 4.0 * f(a+(i as f64)*h);
        } else {
            sum += 2.0 * f(a+(i as f64)*h);
        }
    }
    return h*(sum + f(a) + f(b))/(3.0);
}

/* Adaptive Simpson's rule
 * Function: Computes an integral estimation by using the simspson rule and adapting to some user
 * specified tolerence.
 * Parameters:
 * f - the function
 * a - the inital endpoint
 * b - the final endpoint
 * eps - the tolerence
 * Returns:
 * The result of the integral
 */
pub fn adaptive_simpson(a: f64, b: f64, f: fn(f64) -> f64, eps: f64) -> f64 {
    let fa : f64 = f(a);
    let fb : f64 = f(b);
    let (m, fm, whole) = simpson_single_eval(f, a, b, fa, fb);
    return adaptive_simpson_quadrature(f, a, b, fa, fb, eps, m, fm, whole);
}

/* Adaptive_Simpson_Quadrature
 * Function: Performs the recursive nature of calculating adaptive simpsons rule to a tolerence
 * Parameters:
 * f - the function
 * a - the initial endpoint
 * b - the final endpoint
 * fa - the function evaluated at the initial endpoint
 * fb - the function evaluated at the final endpoint
 * eps - the tolerence
 * m - the midpoint
 * fm - the function evaluated at the midpoint
 * whole - the result of a single eval of simpsons quadrature
 * Returns: 
 * Either a recursive call to itself or the solution to the integral
 */
fn adaptive_simpson_quadrature(f: fn(f64) -> f64, a: f64, b: f64, fa: f64, fb: f64, eps: f64, m: f64, fm: f64, whole: f64) -> f64 {
    let (lm, flm, left) = simpson_single_eval(f, a, m, fa, fm);
    let (rm, frm, right) = simpson_single_eval(f, m, b, fm, fb);
    let delta : f64 = left+right-whole;
    if delta.abs() <= 15.0*eps {
        return left+right+delta/15.0;
    }
    return adaptive_simpson_quadrature(f, a, m, fa, fm, eps/2.0, lm, flm, left) + adaptive_simpson_quadrature(f, m, b, fm, fb, eps/2.0, rm, frm, right);
}

/* Simpson_Single_eval
 * Private Function: Computes simpsons rule at its midpoint and returns the value
 * Parameters:
 * f - the function
 * a - the inital endpoint
 * b - the final endpoint
 * fa - the function evaluated at the initial endpoint
 * fb - the function evaluated at the final endpoint
 * Returns: 
 * The result of a single eval of simpson's rule
 */
fn simpson_single_eval(f: fn(f64) -> f64, a: f64, b: f64, fa: f64, fb: f64) -> (f64, f64, f64) {
    let m : f64 = (a+b)/2.0;
    let fm : f64 = f(m);
    return (m, fm, (b-a).abs() / 6.0 * (fa + 4.0*fm  + fb));
}

/* 3/8 Simpson's Rule
 * Function: Computes an integral estimation using the 3/8s simpsons rule. About twice as accurate
 * as simpson's original rule. 
 * Parameters:
 * f - the function, 
 * n - the number of subintervals,
 * a - the initial endpoint, 
 * b - the final endpoint
 * Returns: The approximate solution to the integral. 
 *
 */

pub fn three_eighths_simpson(a: f64, b: f64, n: i64, f: fn(f64) -> f64) -> f64{
    let h : f64 = (b-a)/(n as f64).abs();
    let mut sum : f64 = 0.0;
    for i in 1..n{
        if i % 3 == 0 {
            sum += 2.0*f(a+(i as f64)*h);
        } else{
            sum += 3.0*f(a+(i as f64)*h);
        }
    }
    return 3.0*h*sum/8.0;
}

/* Unit tests
 * Each method is tested in the integrate class below
 */

#[cfg(test)]
mod tests{
    // imports outer scope integrate class 
    use super::*;
    use assert_approx_eq::assert_approx_eq;

    // square function to use in tests
    fn square(x : f64) -> f64 {
        return x*x;
    }

    #[test]
    fn trapezoid_quad() {
        let a : f64 =  0.66679;
        assert_approx_eq!(a, composite_trapezoid(-1.0, 1.0, 100, square), 2 as f64);
    }

    #[test]
    fn simpson_test() {
        let a : f64 = 0.66679;
        assert_approx_eq!(a, simpson_quadrature(-1.0, 1.0, 100, square), 2 as f64);
    }

    #[test] 
    fn three_eighth_simpson_test() { 
        let a: f64 = 0.66679;
        assert_approx_eq!(a, three_eighths_simpson(-1.0, 1.0, 100, square), 2 as f64);
    }
}