Integral implemented using a trapezoidal Riemann sum

Question

I tried to make this a general purpose integral function but I want to know if it's efficient and idiomatic Rust.

use std::mem;


/// calculates the signed area between the function f and the x axis from
/// x = a to b using a trapezoidal Riemann sum. precision is the number of
/// trapezoids to calculate
pub fn integral<F>(a: f64, b: f64, f: F, precision: u32) -> f64
    where F: Fn(f64) -> f64 {

    let mut a = a;
    let mut b = b;

    let mut sign = 1.0;

    if a > b {
        mem::swap(&mut a, &mut b);
        sign = -1.0;
    }

    let delta = (b - a).abs() / precision as f64;
    let mut result = 0.0;

    for trapezoid in 0..precision {
        let left_side = a + (delta * trapezoid as f64);
        let right_size = left_side + delta;

        result += 0.5 * (f(left_side) + f(right_size)) * delta;
    }

    result * sign
}

And here's a test case

fn f(x: f64) -> f64 {
    (3.0 * x * x) + (4.0 * x) + 7.0
}

fn main() {
    let a = integral(0.0, 11.5, f, 1000000);
    println!("{}", a); // expect approx 1865.88
}

Answer 1

1. Embrace even more static analysis tools, such as clippy. It provides warnings like:

   warning: `if _ { .. } else { .. }` is an expression
     --> src/main.rs:12:5
      |
   12 |       let mut sign = 1.0;
      |  _____^ starting here...
   13 | |
   14 | |     if a > b {
   15 | |         mem::swap(&mut a, &mut b);
   16 | |         sign = -1.0;
   17 | |     }
      | |_____^ ...ending here
      |
   

2. You don't need to rebind a variable to make it mutable. you can just add mut in the argument list.

3. When a where clause is used, the { moves to the next line.

4. Why is abs used? Isn't it guaranteed that b > a?

5. "side" vs "size" seems likely to cause confusion. Either make them the same or distinct.

6. Using map and sum avoids needing to make result mutable.

7. Consider removing all mutability; either by reassigning all the variables at once:

   let (a, b, sign) = if a > b {
       (b, a, -1.0)
   } else {
       (a, b, 1.0)
   };
   

   Or by adding a small shim that reverses the order and applies the negation; I like the latter solution because it ties the argument reversal and final negation closer together.

/// calculates the signed area between the function f and the x axis from
/// x = a to b using a trapezoidal Riemann sum. precision is the number of
/// trapezoids to calculate
pub fn integral<F>(a: f64, b: f64, f: F, precision: u32) -> f64
    where F: Fn(f64) -> f64
{
    fn core<F>(a: f64, b: f64, f: F, precision: u32) -> f64
        where F: Fn(f64) -> f64
    {
        let delta = (b - a) / precision as f64;

        (0..precision).map(|trapezoid| {
            let left_side = a + (delta * trapezoid as f64);
            let right_size = left_side + delta;

            0.5 * (f(left_side) + f(right_size)) * delta
        }).sum()
    }

    if a > b {
        -core(b, a, f, precision)
    } else {
        core(a, b, f, precision)
    }
}