qualia_core_db/solvers/vector_calculus/
differential.rs1use 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
10pub fn gradient(f: &Expr, vars: &[&str]) -> Vec<Expr> {
12 vars.iter().map(|v| partial(f, v)).collect()
13}
14
15pub 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
27pub 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
42pub 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 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 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 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 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}