Skip to main content

qualia_core_db/solvers/vector_calculus/
differential.rs

1//! Symbolic differential operators over the CAS, built on `differentiate`/`simplify`.
2
3use super::VecCalcError;
4use crate::specialized_libs::symbolic_algebra::{add, differentiate, simplify, sub, Expr};
5
6fn partial(f: &Expr, v: &str) -> Expr {
7    simplify(&differentiate(f, v))
8}
9
10/// `∇f = [∂f/∂x₁, …, ∂f/∂xₙ]`.
11pub fn gradient(f: &Expr, vars: &[&str]) -> Vec<Expr> {
12    vars.iter().map(|v| partial(f, v)).collect()
13}
14
15/// Divergence `∇·F = Σ ∂Fᵢ/∂xᵢ`. Requires `field.len() == vars.len()`.
16pub fn divergence(field: &[Expr], vars: &[&str]) -> Result<Expr, VecCalcError> {
17    if field.len() != vars.len() || field.is_empty() {
18        return Err(VecCalcError::DimensionMismatch);
19    }
20    let mut acc = partial(&field[0], vars[0]);
21    for i in 1..field.len() {
22        acc = add(acc, partial(&field[i], vars[i]));
23    }
24    Ok(simplify(&acc))
25}
26
27/// Curl `∇×F` of a 3-D field. Requires exactly 3 components and 3 variables `[x,y,z]`.
28/// Returns `[ ∂Fz/∂y−∂Fy/∂z , ∂Fx/∂z−∂Fz/∂x , ∂Fy/∂x−∂Fx/∂y ]`.
29pub fn curl(field: &[Expr], vars: &[&str]) -> Result<[Expr; 3], VecCalcError> {
30    if field.len() != 3 || vars.len() != 3 {
31        return Err(VecCalcError::DimensionMismatch);
32    }
33    let (fx, fy, fz) = (&field[0], &field[1], &field[2]);
34    let (x, y, z) = (vars[0], vars[1], vars[2]);
35    Ok([
36        simplify(&sub(partial(fz, y), partial(fy, z))),
37        simplify(&sub(partial(fx, z), partial(fz, x))),
38        simplify(&sub(partial(fy, x), partial(fx, y))),
39    ])
40}
41
42/// Laplacian `∇²f = Σ ∂²f/∂xᵢ²`.
43pub fn laplacian(f: &Expr, vars: &[&str]) -> Result<Expr, VecCalcError> {
44    if vars.is_empty() {
45        return Err(VecCalcError::DimensionMismatch);
46    }
47    let second = |v: &str| partial(&partial(f, v), v);
48    let mut acc = second(vars[0]);
49    for v in &vars[1..] {
50        acc = add(acc, second(v));
51    }
52    Ok(simplify(&acc))
53}
54
55#[cfg(test)]
56mod tests {
57    use super::*;
58    use crate::specialized_libs::symbolic_algebra::{add as eadd, mul, pow, var};
59    use std::collections::HashMap;
60
61    fn at(e: &Expr, p: &[(&str, f64)]) -> f64 {
62        let env: HashMap<String, f64> = p.iter().map(|&(k, v)| (k.to_string(), v)).collect();
63        e.eval(&env).unwrap()
64    }
65
66    #[test]
67    fn divergence_of_position_field_is_three() {
68        // F = (x, y, z) → ∇·F = 3
69        let field = [var("x"), var("y"), var("z")];
70        let d = divergence(&field, &["x", "y", "z"]).unwrap();
71        assert!((at(&d, &[("x", 1.0), ("y", 2.0), ("z", 3.0)]) - 3.0).abs() < 1e-9);
72    }
73
74    #[test]
75    fn curl_of_a_gradient_is_zero() {
76        // f = x²y + z ; ∇f then curl(∇f) = 0.
77        let f = eadd(mul(pow(var("x"), 2), var("y")), var("z"));
78        let g = gradient(&f, &["x", "y", "z"]);
79        let cc = curl(&g, &["x", "y", "z"]).unwrap();
80        for comp in &cc {
81            assert!(at(comp, &[("x", 1.3), ("y", -0.7), ("z", 2.0)]).abs() < 1e-9);
82        }
83    }
84
85    #[test]
86    fn divergence_of_a_curl_is_zero() {
87        // F = (x²z, x y², y z²) ; ∇·(∇×F) = 0.
88        let field = [
89            mul(pow(var("x"), 2), var("z")),
90            mul(var("x"), pow(var("y"), 2)),
91            mul(var("y"), pow(var("z"), 2)),
92        ];
93        let c = curl(&field, &["x", "y", "z"]).unwrap();
94        let d = divergence(&c, &["x", "y", "z"]).unwrap();
95        assert!(at(&d, &[("x", 0.9), ("y", 1.1), ("z", -0.4)]).abs() < 1e-8);
96    }
97
98    #[test]
99    fn laplacian_of_r_squared_is_six() {
100        // f = x²+y²+z² → ∇²f = 6
101        let f = eadd(eadd(pow(var("x"), 2), pow(var("y"), 2)), pow(var("z"), 2));
102        let l = laplacian(&f, &["x", "y", "z"]).unwrap();
103        assert!((at(&l, &[("x", 5.0), ("y", -2.0), ("z", 1.0)]) - 6.0).abs() < 1e-9);
104    }
105
106    #[test]
107    fn curl_dimension_guard() {
108        assert_eq!(
109            curl(&[var("x"), var("y")], &["x", "y"]).unwrap_err(),
110            VecCalcError::DimensionMismatch
111        );
112    }
113}