qualia_core_db/solvers/vector_calculus/
integrals.rs

1//! Numeric line and surface integrals over parametric curves/surfaces given as
2//! closures. Derivatives of the parametrization are taken by central finite difference,
3//! so the caller supplies only the position map. Trapezoidal quadrature.
4
5const FD: f64 = 1e-6;
6
7fn dot(a: &[f64], b: &[f64]) -> f64 {
8    a.iter().zip(b).map(|(x, y)| x * y).sum()
9}
10fn norm(a: &[f64]) -> f64 {
11    dot(a, a).sqrt()
12}
13fn cross3(a: &[f64], b: &[f64]) -> [f64; 3] {
14    [
15        a[1] * b[2] - a[2] * b[1],
16        a[2] * b[0] - a[0] * b[2],
17        a[0] * b[1] - a[1] * b[0],
18    ]
19}
20
21/// Central-difference derivative of a curve `r(t)`.
22fn dcurve<C: Fn(f64) -> Vec<f64>>(curve: &C, t: f64) -> Vec<f64> {
23    let p = curve(t + FD);
24    let m = curve(t - FD);
25    p.iter()
26        .zip(&m)
27        .map(|(a, b)| (a - b) / (2.0 * FD))
28        .collect()
29}
30
31/// Scalar line integral `∫_C f ds = ∫ f(r(t)) |r'(t)| dt`.
32pub fn line_integral_scalar<F, C>(f: F, curve: C, t0: f64, t1: f64, steps: usize) -> f64
33where
34    F: Fn(&[f64]) -> f64,
35    C: Fn(f64) -> Vec<f64>,
36{
37    let h = (t1 - t0) / steps as f64;
38    let g = |t: f64| f(&curve(t)) * norm(&dcurve(&curve, t));
39    let mut sum = 0.5 * (g(t0) + g(t1));
40    for i in 1..steps {
41        sum += g(t0 + i as f64 * h);
42    }
43    sum * h
44}
45
46/// Work / vector line integral `∫_C F·dr = ∫ F(r(t))·r'(t) dt`.
47pub fn line_integral_work<F, C>(field: F, curve: C, t0: f64, t1: f64, steps: usize) -> f64
48where
49    F: Fn(&[f64]) -> Vec<f64>,
50    C: Fn(f64) -> Vec<f64>,
51{
52    let h = (t1 - t0) / steps as f64;
53    let g = |t: f64| dot(&field(&curve(t)), &dcurve(&curve, t));
54    let mut sum = 0.5 * (g(t0) + g(t1));
55    for i in 1..steps {
56        sum += g(t0 + i as f64 * h);
57    }
58    sum * h
59}
60
61/// Flux of a 3-D field through a parametric surface `r(u,v)`:
62/// `∫∫ F·(r_u × r_v) du dv`. The orientation is that of `r_u × r_v`.
63pub fn surface_flux<F, S>(
64    field: F,
65    surf: S,
66    u0: f64,
67    u1: f64,
68    v0: f64,
69    v1: f64,
70    steps: usize,
71) -> f64
72where
73    F: Fn(&[f64]) -> Vec<f64>,
74    S: Fn(f64, f64) -> Vec<f64>,
75{
76    let hu = (u1 - u0) / steps as f64;
77    let hv = (v1 - v0) / steps as f64;
78    let integrand = |u: f64, v: f64| -> f64 {
79        let ru: Vec<f64> = {
80            let p = surf(u + FD, v);
81            let m = surf(u - FD, v);
82            p.iter()
83                .zip(&m)
84                .map(|(a, b)| (a - b) / (2.0 * FD))
85                .collect()
86        };
87        let rv: Vec<f64> = {
88            let p = surf(u, v + FD);
89            let m = surf(u, v - FD);
90            p.iter()
91                .zip(&m)
92                .map(|(a, b)| (a - b) / (2.0 * FD))
93                .collect()
94        };
95        let n = cross3(&ru, &rv);
96        dot(&field(&surf(u, v)), &n)
97    };
98    // Trapezoidal over the 2-D grid.
99    let mut sum = 0.0;
100    for i in 0..=steps {
101        for j in 0..=steps {
102            let w = {
103                let wu = if i == 0 || i == steps { 0.5 } else { 1.0 };
104                let wv = if j == 0 || j == steps { 0.5 } else { 1.0 };
105                wu * wv
106            };
107            sum += w * integrand(u0 + i as f64 * hu, v0 + j as f64 * hv);
108        }
109    }
110    sum * hu * hv
111}
112
113#[cfg(test)]
114mod tests {
115    use super::*;
116    use core::f64::consts::PI;
117
118    #[test]
119    fn conservative_work_is_potential_difference() {
120        // F = (y, x) = ∇(xy). Work from (0,0)→(1,1) along the diagonal = xy|=1.
121        let f = |p: &[f64]| vec![p[1], p[0]];
122        let curve = |t: f64| vec![t, t];
123        let w = line_integral_work(f, curve, 0.0, 1.0, 400);
124        assert!((w - 1.0).abs() < 1e-6);
125    }
126
127    #[test]
128    fn greens_theorem_on_the_unit_circle() {
129        // ∮ F·dr with F = (−y, x) around the unit circle = 2·area = 2π.
130        let f = |p: &[f64]| vec![-p[1], p[0]];
131        let circle = |t: f64| vec![t.cos(), t.sin()];
132        let w = line_integral_work(f, circle, 0.0, 2.0 * PI, 2000);
133        assert!((w - 2.0 * PI).abs() < 1e-4);
134    }
135
136    #[test]
137    fn arc_length_via_scalar_line_integral() {
138        // ∫_C 1 ds over the unit circle = circumference = 2π.
139        let len = line_integral_scalar(
140            |_| 1.0,
141            |t: f64| vec![t.cos(), t.sin()],
142            0.0,
143            2.0 * PI,
144            2000,
145        );
146        assert!((len - 2.0 * PI).abs() < 1e-4);
147    }
148
149    #[test]
150    fn divergence_theorem_flux_through_sphere() {
151        // Flux of F = (x,y,z) through the unit sphere = ∫∫∫ div F dV = 3·(4/3 π) = 4π.
152        // Parametrize with polar angle u∈[0,π] first, azimuth v∈[0,2π] second, so
153        // r_u × r_v is the *outward* normal (the divergence theorem's orientation).
154        let f = |p: &[f64]| vec![p[0], p[1], p[2]];
155        let sphere = |u: f64, v: f64| vec![u.sin() * v.cos(), u.sin() * v.sin(), u.cos()];
156        let flux = surface_flux(f, sphere, 0.0, PI, 0.0, 2.0 * PI, 120);
157        assert!((flux - 4.0 * PI).abs() < 1e-2, "flux {flux} vs 4π");
158    }
159}
qualia_core_db/solvers/vector_calculus/integrals.rs

qualia_core_db/solvers/vector_calculus/
integrals.rs
