40_trapezoidal_integration/
40_trapezoidal_integration.rs

1//! # Example: Trapezoidal Integration
2//!
3//! Run: cargo run --example 40_trapezoidal_integration
4//!
5//! ## Problem
6//! Approximate the area under a curve (a definite integral) from sampled function
7//! values, and compare against the known exact value.
8//!
9//! ## Math idea
10//! The composite trapezoidal rule approximates `∫ f` by summing trapezoids between
11//! adjacent samples:
12//! `∫ ≈ h · (f₀/2 + f₁ + … + f_{n-1} + fₙ/2)`,
13//! which equals `h · (Σ f − (f₀ + fₙ)/2)`. Here `f(x) = x²` on `[0, 1]`, whose exact
14//! integral is `1/3`.
15//!
16//! ## Tensor representation
17//! The sample grid is a 1-D `Tensor` from `linspace`; the function values are a
18//! `Tensor` computed elementwise, and the running total uses a `Tensor` reduction.
19//!
20//! ## What this demonstrates
21//! - `Tensor::linspace` for an evenly spaced grid;
22//! - elementwise squaring (`&x * &x`) to evaluate `f`;
23//! - a whole-tensor `sum` reduction in the trapezoidal formula.
24//!
25//! ## Expected output
26//! ```text
27//! trapezoidal estimate: 0.335000
28//! exact integral (1/3): 0.333333
29//! abs error:            0.001667
30//! Trapezoidal integration: OK
31//! ```
32//!
33//! This is a numerical approximation, not an integration library.
34
35use matten::Tensor;
36
37fn main() {
38    // Evenly spaced sample grid on [0, 1].
39    let x = Tensor::linspace(0.0, 1.0, 11);
40    let xs = x.as_slice();
41    let h = xs[1] - xs[0];
42
43    // f(x) = x^2, evaluated with elementwise multiplication.
44    let f = &x * &x;
45    let fs = f.as_slice();
46    let n = fs.len();
47
48    // Composite trapezoidal rule: h * (Σ f - (f_first + f_last)/2).
49    let estimate = h * (f.sum() - (fs[0] + fs[n - 1]) / 2.0);
50    let exact = 1.0 / 3.0; // ∫_0^1 x^2 dx = 1/3
51    let err = (estimate - exact).abs();
52
53    println!("trapezoidal estimate: {estimate:.6}");
54    println!("exact integral (1/3): {exact:.6}");
55    println!("abs error:            {err:.6}");
56
57    assert!(err < 0.01);
58    println!("Trapezoidal integration: OK");
59}
40_trapezoidal_integration/40_trapezoidal_integration.rs

40_trapezoidal_integration/
40_trapezoidal_integration.rs
