qualia_core_db/solvers/transforms/
laplace.rs

1//! Laplace transform `L{f}(s) = ∫₀^∞ e^{−st} f(t) dt`.
2//!
3//! Two paths: a **numerical** transform by Simpson quadrature for any closure (general,
4//! always available), and a **symbolic** table transform over the CAS expression type.
5//! The current `Expr` algebra has no exp/trig variants, so the symbolic table covers the
6//! cases it *can* represent — constants, integer powers `tⁿ`, and their linear
7//! combinations — and **fails closed** (`NotTransformable`) on anything else rather than
8//! returning a wrong transform.
9
10use crate::specialized_libs::symbolic_algebra::{add, c, div, neg, pow, sub, var, Expr};
11
12#[derive(Debug, Clone, PartialEq, Eq)]
13pub enum LaplaceError {
14    /// The expression is outside the symbolic table (e.g. needs exp/trig the CAS lacks).
15    NotTransformable,
16    /// Domain error in the numeric transform (`s ≤ 0`, non-positive horizon, …).
17    OutOfDomain,
18}
19
20/// Numerical Laplace transform of `f` at `s`, integrating to `t_max` with `steps` (even)
21/// Simpson panels. Requires `s > 0`, `t_max > 0`, `steps ≥ 2` even.
22pub fn laplace_numeric<F: Fn(f64) -> f64>(
23    f: F,
24    s: f64,
25    t_max: f64,
26    steps: usize,
27) -> Result<f64, LaplaceError> {
28    if s <= 0.0 || t_max <= 0.0 || steps < 2 || steps % 2 != 0 {
29        return Err(LaplaceError::OutOfDomain);
30    }
31    let h = t_max / steps as f64;
32    let g = |t: f64| (-s * t).exp() * f(t);
33    let mut sum = g(0.0) + g(t_max);
34    for i in 1..steps {
35        let t = i as f64 * h;
36        sum += if i % 2 == 1 { 4.0 } else { 2.0 } * g(t);
37    }
38    Ok(sum * h / 3.0)
39}
40
41fn factorial(n: i32) -> f64 {
42    (1..=n).fold(1.0, |a, k| a * k as f64)
43}
44
45/// Symbolic Laplace transform of `expr` in the time variable `t`, returning an `Expr` in
46/// the complex frequency variable `s`. Handles constants (`c → c/s`), powers
47/// (`tⁿ → n!/s^{n+1}`), negation, sums/differences, and scalar multiples; everything
48/// else fails closed.
49pub fn laplace_table(expr: &Expr) -> Result<Expr, LaplaceError> {
50    transform(expr, "t")
51}
52
53fn transform(expr: &Expr, t: &str) -> Result<Expr, LaplaceError> {
54    match expr {
55        Expr::Const(k) => Ok(div(c(*k), var("s"))), // L{k} = k/s
56        Expr::Var(name) if name == t => Ok(div(c(1.0), pow(var("s"), 2))), // L{t} = 1/s²
57        Expr::Pow(base, n) if is_var(base, t) && *n >= 0 => {
58            // L{tⁿ} = n!/s^{n+1}
59            Ok(div(c(factorial(*n)), pow(var("s"), n + 1)))
60        }
61        Expr::Neg(a) => Ok(neg(transform(a, t)?)),
62        Expr::Add(a, b) => Ok(add(transform(a, t)?, transform(b, t)?)),
63        Expr::Sub(a, b) => Ok(sub(transform(a, t)?, transform(b, t)?)),
64        Expr::Mul(a, b) => {
65            // Scalar · f(t): L is the scalar times L{f}.
66            if let Expr::Const(_) = **a {
67                Ok(crate::specialized_libs::symbolic_algebra::mul(
68                    (**a).clone(),
69                    transform(b, t)?,
70                ))
71            } else if let Expr::Const(_) = **b {
72                Ok(crate::specialized_libs::symbolic_algebra::mul(
73                    (**b).clone(),
74                    transform(a, t)?,
75                ))
76            } else {
77                Err(LaplaceError::NotTransformable)
78            }
79        }
80        _ => Err(LaplaceError::NotTransformable),
81    }
82}
83
84fn is_var(e: &Expr, name: &str) -> bool {
85    matches!(e, Expr::Var(v) if v == name)
86}
87
88#[cfg(test)]
89mod tests {
90    use super::*;
91    use std::collections::HashMap;
92
93    fn eval_at_s(e: &Expr, s: f64) -> f64 {
94        let mut env = HashMap::new();
95        env.insert("s".to_string(), s);
96        e.eval(&env).unwrap()
97    }
98
99    #[test]
100    fn numeric_transforms_match_closed_forms() {
101        // L{1}(s) = 1/s
102        assert!((laplace_numeric(|_| 1.0, 2.0, 60.0, 4000).unwrap() - 0.5).abs() < 1e-4);
103        // L{e^{-t}}(s=1) = 1/(s+1) = 1/2
104        assert!((laplace_numeric(|t| (-t).exp(), 1.0, 60.0, 4000).unwrap() - 0.5).abs() < 1e-4);
105        // L{t}(s=1) = 1/s² = 1
106        assert!((laplace_numeric(|t| t, 1.0, 80.0, 8000).unwrap() - 1.0).abs() < 1e-3);
107    }
108
109    #[test]
110    fn symbolic_table_powers_and_linearity() {
111        // L{t²} = 2/s³ ; at s=2 → 2/8 = 0.25
112        let l = laplace_table(&pow(var("t"), 2)).unwrap();
113        assert!((eval_at_s(&l, 2.0) - 0.25).abs() < 1e-12);
114        // L{3} = 3/s ; at s=3 → 1
115        let l2 = laplace_table(&c(3.0)).unwrap();
116        assert!((eval_at_s(&l2, 3.0) - 1.0).abs() < 1e-12);
117        // L{t + 5} = 1/s² + 5/s ; at s=1 → 6
118        let l3 = laplace_table(&add(var("t"), c(5.0))).unwrap();
119        assert!((eval_at_s(&l3, 1.0) - 6.0).abs() < 1e-12);
120    }
121
122    #[test]
123    fn fails_closed_on_unrepresentable() {
124        // sqrt(t) is not in the polynomial table.
125        let e = crate::specialized_libs::symbolic_algebra::sqrt(var("t"));
126        assert_eq!(
127            laplace_table(&e).unwrap_err(),
128            LaplaceError::NotTransformable
129        );
130        assert_eq!(
131            laplace_numeric(|_| 1.0, -1.0, 10.0, 100).unwrap_err(),
132            LaplaceError::OutOfDomain
133        );
134    }
135}
qualia_core_db/solvers/transforms/laplace.rs

qualia_core_db/solvers/transforms/
laplace.rs
