r2rs_mass/func/
area.rs

1/// # Adaptive Numerical Integration
2///
3/// ## Description:
4///
5/// Integrate a function of one variable over a finite range using a
6/// recursive adaptive method.  This function is mainly for
7/// demonstration purposes.
8///
9/// ## Usage:
10///
11/// area(f, a, b, ..., fa = f(a, ...), fb = f(b, ...),
12///  limit = 10, eps = 1e-05)
13///
14/// ## Arguments:
15///
16/// * f: The integrand as an ‘S’ function object.  The variable of
17///  integration must be the first argument.
18/// * a: Lower limit of integration.
19/// * b: Upper limit of integration.
20/// * ...: Additional arguments needed by the integrand.
21/// * fa: Function value at the lower limit.
22/// * fb: Function value at the upper limit.
23/// * limit: Limit on the depth to which recursion is allowed to go.
24/// * eps: Error tolerance to control the process.
25///
26/// ## Details:
27///
28/// The method divides the interval in two and compares the values
29/// given by Simpson's rule and the trapezium rule.  If these are
30/// within eps of each other the Simpson's rule result is given,
31/// otherwise the process is applied separately to each half of the
32/// interval and the results added together.
33///
34/// ## Value:
35///
36/// The integral from ‘a’ to ‘b’ of ‘f(x)’.
37///
38/// ## References:
39///
40/// Venables, W. N. and Ripley, B. D. (1994) _Modern Applied
41/// Statistics with S-Plus._ Springer.  pp. 105-110.
42///
43/// ## Examples:
44///
45/// ```r
46/// area(sin, 0, pi)  # integrate the sin function from 0 to pi.
47/// ```
48pub fn area(
49    f: &dyn Fn(f64) -> f64,
50    lower_limit: f64,
51    upper_limit: f64,
52    max_iterations: usize,
53    epsilon: f64,
54) -> f64 {
55    let h = upper_limit - lower_limit;
56    let d = (lower_limit + upper_limit) / 2.0;
57    let fa = f(lower_limit);
58    let fb = f(upper_limit);
59    let fd = f(d);
60    let a1 = ((fa + fb) * h) / 2.0;
61    let a2 = ((fa + 4.0 * fd + fb) * h) / 6.0;
62    if (a1 - a2).abs() < epsilon || max_iterations == 0 {
63        a2
64    } else {
65        area(f, lower_limit, d, max_iterations - 1, epsilon)
66            + area(f, d, upper_limit, max_iterations - 1, epsilon)
67    }
68}
r2rs_mass/func/area.rs

r2rs_mass/func/
area.rs
