# Integral implemented using a trapezoidal Riemann sum

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
}


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)
}
}

• You know, I didn't even notice left_side vs right_size. Editor filled it in for me. Apr 21, 2017 at 3:50