qualia_core_db/specialized_libs/
symbolic_ode.rs

1//! **Symbolic differential equations** (Gap analysis §3.4) — closed-form solutions for the
2//! standard *solvable* classes of ODE, plus first-order-linear PDE (method of
3//! characteristics) and second-order-linear PDE type classification.
4//!
5//! Honest scope (everything outside it returns [`OdeError::NotSupported`] — never a
6//! fabricated solution):
7//!
8//! **ODE**
9//! - **Separable** `y' = g(x)·h(y)` → implicit `∫dy/h(y) = ∫g(x)dx + C`, using the CAS
10//!   integrator ([`crate::specialized_libs::symbolic_integration`]); fails closed when
11//!   either integral is outside the integrator's table.
12//! - **Linear first-order, constant coefficients** `y' + a·y = b` → explicit.
13//! - **Linear second-order, constant coefficients** `a·y'' + b·y' + c·y = 0` → explicit, via
14//!   the characteristic equation (distinct-real / repeated / complex roots).
15//!
16//! **PDE**
17//! - **First-order linear homogeneous** `a·uₓ + b·u_y = 0` → `u = F(b·x − a·y)` (an arbitrary
18//!   differentiable `F`; the characteristic invariant is returned).
19//! - **Second-order linear** `A·uₓₓ + B·u_xy + C·u_yy + … ` → elliptic / parabolic / hyperbolic
20//!   classification by the discriminant `B² − 4AC`.
21//!
22//! A general nonlinear/variable-coefficient PDE solver is *not* attempted — that is genuinely
23//! beyond a bounded module, and the contract here is "solve the supported classes exactly,
24//! refuse the rest", not "pretend".
25
26use super::symbolic_algebra::{add, c, cos, div, exp, mul, sin, var, Expr};
27use super::symbolic_integration::{integrate, IntegrationError};
28
29#[derive(Debug, Clone, PartialEq, Eq)]
30pub enum OdeError {
31    /// The equation is outside the supported solvable classes.
32    NotSupported,
33    /// A required antiderivative is outside the CAS integrator's table.
34    NotIntegrable,
35}
36
37impl From<IntegrationError> for OdeError {
38    fn from(_: IntegrationError) -> Self {
39        OdeError::NotIntegrable
40    }
41}
42
43/// The solution of an ODE.
44#[derive(Debug, Clone, PartialEq)]
45pub enum OdeSolution {
46    /// `y(x)` given explicitly; integration constants appear as variables `C`, `C1`, `C2`.
47    Explicit(Expr),
48    /// An implicit relation `F(y) = G(x)` (with `G` already including `+ C`).
49    Implicit { f_y: Expr, g_x: Expr },
50}
51
52/// Solve the **separable** ODE `y' = g(x)·h(y)` as `∫ dy/h(y) = ∫ g(x) dx + C`. The two
53/// antiderivatives are computed by the CAS integrator; either being non-integrable yields
54/// [`OdeError::NotIntegrable`].
55pub fn solve_separable(
56    g_x: &Expr,
57    h_y: &Expr,
58    xvar: &str,
59    yvar: &str,
60) -> Result<OdeSolution, OdeError> {
61    let f_y = integrate(&div(c(1.0), h_y.clone()), yvar)?; // ∫ dy / h(y)
62    let g_int = integrate(g_x, xvar)?; // ∫ g(x) dx
63    Ok(OdeSolution::Implicit {
64        f_y,
65        g_x: add(g_int, var("C")),
66    })
67}
68
69/// Solve the **linear first-order constant-coefficient** ODE `y' + a·y = b`.
70/// - `a ≠ 0` → `y = b/a + C·e^{−a x}`.
71/// - `a = 0` → `y = b·x + C`.
72pub fn solve_linear_first_order(a: f64, b: f64, xvar: &str) -> OdeSolution {
73    if a == 0.0 {
74        OdeSolution::Explicit(add(mul(c(b), var(xvar)), var("C")))
75    } else {
76        let homogeneous = mul(var("C"), exp(mul(c(-a), var(xvar))));
77        OdeSolution::Explicit(add(c(b / a), homogeneous))
78    }
79}
80
81/// Solve the **linear second-order homogeneous constant-coefficient** ODE
82/// `a·y'' + b·y' + c·y = 0` via its characteristic equation `a·r² + b·r + c = 0`. `a = 0`
83/// is not second-order → [`OdeError::NotSupported`].
84pub fn solve_linear_second_order(
85    a: f64,
86    b: f64,
87    cc: f64,
88    xvar: &str,
89) -> Result<OdeSolution, OdeError> {
90    if a == 0.0 {
91        return Err(OdeError::NotSupported);
92    }
93    let disc = b * b - 4.0 * a * cc;
94    let x = var(xvar);
95    let sol = if disc > 1e-12 {
96        // Distinct real roots: y = C1·e^{r1 x} + C2·e^{r2 x}.
97        let s = disc.sqrt();
98        let r1 = (-b + s) / (2.0 * a);
99        let r2 = (-b - s) / (2.0 * a);
100        add(
101            mul(var("C1"), exp(mul(c(r1), x.clone()))),
102            mul(var("C2"), exp(mul(c(r2), x))),
103        )
104    } else if disc.abs() <= 1e-12 {
105        // Repeated root r: y = (C1 + C2·x)·e^{r x}.
106        let r = -b / (2.0 * a);
107        mul(add(var("C1"), mul(var("C2"), x.clone())), exp(mul(c(r), x)))
108    } else {
109        // Complex roots α ± βi: y = e^{α x}·(C1·cos(β x) + C2·sin(β x)).
110        let alpha = -b / (2.0 * a);
111        let beta = (-disc).sqrt() / (2.0 * a);
112        mul(
113            exp(mul(c(alpha), x.clone())),
114            add(
115                mul(var("C1"), cos(mul(c(beta), x.clone()))),
116                mul(var("C2"), sin(mul(c(beta), x))),
117            ),
118        )
119    };
120    Ok(OdeSolution::Explicit(sol))
121}
122
123// ── PDE ──────────────────────────────────────────────────────────────────────────
124
125/// The solution of a (supported) PDE.
126#[derive(Debug, Clone, PartialEq)]
127pub enum PdeSolution {
128    /// `u(x, y) = F(invariant)` for an arbitrary differentiable `F` (method of
129    /// characteristics for a first-order linear homogeneous PDE).
130    GeneralFunctionOf { invariant: Expr },
131}
132
133/// The type of a second-order linear PDE.
134#[derive(Debug, Clone, Copy, PartialEq, Eq)]
135pub enum PdeClass {
136    /// `B² − 4AC < 0` (e.g. Laplace `uₓₓ + u_yy = 0`).
137    Elliptic,
138    /// `B² − 4AC = 0` (e.g. the heat equation).
139    Parabolic,
140    /// `B² − 4AC > 0` (e.g. the wave equation `uₓₓ − u_yy = 0`).
141    Hyperbolic,
142}
143
144/// Solve `a·uₓ + b·u_y = 0` by the method of characteristics: `u` is an arbitrary function
145/// of the invariant `b·x − a·y`. Requires `(a, b) ≠ (0, 0)`.
146pub fn solve_first_order_linear_pde(
147    a: f64,
148    b: f64,
149    xvar: &str,
150    yvar: &str,
151) -> Result<PdeSolution, OdeError> {
152    if a == 0.0 && b == 0.0 {
153        return Err(OdeError::NotSupported);
154    }
155    // Characteristic invariant ξ = b·x − a·y (constant along characteristics).
156    let invariant = add(mul(c(b), var(xvar)), mul(c(-a), var(yvar)));
157    Ok(PdeSolution::GeneralFunctionOf { invariant })
158}
159
160/// Classify the second-order linear PDE `A·uₓₓ + B·u_xy + C·u_yy + … = …` by the
161/// discriminant `B² − 4AC`.
162pub fn classify_second_order_pde(a_xx: f64, b_xy: f64, c_yy: f64) -> PdeClass {
163    let disc = b_xy * b_xy - 4.0 * a_xx * c_yy;
164    if disc < -1e-12 {
165        PdeClass::Elliptic
166    } else if disc <= 1e-12 {
167        PdeClass::Parabolic
168    } else {
169        PdeClass::Hyperbolic
170    }
171}
172
173#[cfg(test)]
174mod tests {
175    use super::super::symbolic_algebra::{differentiate, pow, simplify, var};
176    use super::*;
177    use std::collections::HashMap;
178
179    fn env(pairs: &[(&str, f64)]) -> HashMap<String, f64> {
180        pairs.iter().map(|&(k, v)| (k.to_string(), v)).collect()
181    }
182
183    #[test]
184    fn separable_growth() {
185        // y' = 1·y  (g = 1, h = y)  →  ∫dy/y = ∫dx + C  →  ln(y) = x + C.
186        let sol = solve_separable(&c(1.0), &var("y"), "x", "y").unwrap();
187        match sol {
188            OdeSolution::Implicit { f_y, g_x } => {
189                // f_y = ln(y) ; evaluate F(y=e) = 1.
190                assert!(
191                    (f_y.eval(&env(&[("y", std::f64::consts::E)])).unwrap() - 1.0).abs() < 1e-9
192                );
193                // g_x = x + C ; at x=2, C=3 → 5.
194                assert!((g_x.eval(&env(&[("x", 2.0), ("C", 3.0)])).unwrap() - 5.0).abs() < 1e-9);
195            }
196            _ => panic!("expected implicit solution"),
197        }
198    }
199
200    #[test]
201    fn separable_fails_closed() {
202        // h(y) = sin(y²) → ∫dy/sin(y²) is not in the table → NotIntegrable, not fabricated.
203        let h = super::super::symbolic_algebra::sin(pow(var("y"), 2));
204        assert_eq!(
205            solve_separable(&c(1.0), &h, "x", "y").unwrap_err(),
206            OdeError::NotIntegrable
207        );
208    }
209
210    #[test]
211    fn linear_first_order_satisfies_the_ode() {
212        // y' + 2y = 6  →  y = 3 + C·e^{−2x}. Check y' + 2y = 6 at samples with C set.
213        let OdeSolution::Explicit(y) = solve_linear_first_order(2.0, 6.0, "x") else {
214            panic!()
215        };
216        let yp = simplify(&differentiate(&y, "x"));
217        for &(x, cval) in &[(0.0, 1.0), (0.7, -2.0), (1.5, 4.0)] {
218            let e = env(&[("x", x), ("C", cval)]);
219            let residual = yp.eval(&e).unwrap() + 2.0 * y.eval(&e).unwrap();
220            assert!(
221                (residual - 6.0).abs() < 1e-7,
222                "residual {residual} at x={x}"
223            );
224        }
225    }
226
227    /// Substitute an explicit solution into `a·y'' + b·y' + c·y` and assert ≈ 0.
228    fn verify_second_order(y: &Expr, a: f64, b: f64, cc: f64) {
229        let yp = simplify(&differentiate(y, "x"));
230        let ypp = simplify(&differentiate(&yp, "x"));
231        for &(x, c1, c2) in &[(0.0, 1.0, 0.5), (0.8, -1.0, 2.0), (1.7, 3.0, -1.5)] {
232            let e = env(&[("x", x), ("C1", c1), ("C2", c2)]);
233            let r = a * ypp.eval(&e).unwrap() + b * yp.eval(&e).unwrap() + cc * y.eval(&e).unwrap();
234            assert!(r.abs() < 1e-6, "residual {r} at x={x}");
235        }
236    }
237
238    #[test]
239    fn second_order_distinct_real_roots() {
240        // y'' − 3y' + 2y = 0 → roots 1, 2 → C1 e^x + C2 e^{2x}.
241        let OdeSolution::Explicit(y) = solve_linear_second_order(1.0, -3.0, 2.0, "x").unwrap()
242        else {
243            panic!()
244        };
245        verify_second_order(&y, 1.0, -3.0, 2.0);
246    }
247
248    #[test]
249    fn second_order_repeated_root() {
250        // y'' − 2y' + y = 0 → repeated root 1 → (C1 + C2 x) e^x.
251        let OdeSolution::Explicit(y) = solve_linear_second_order(1.0, -2.0, 1.0, "x").unwrap()
252        else {
253            panic!()
254        };
255        verify_second_order(&y, 1.0, -2.0, 1.0);
256    }
257
258    #[test]
259    fn second_order_complex_roots() {
260        // y'' + y = 0 → roots ±i → C1 cos x + C2 sin x.
261        let OdeSolution::Explicit(y) = solve_linear_second_order(1.0, 0.0, 1.0, "x").unwrap()
262        else {
263            panic!()
264        };
265        verify_second_order(&y, 1.0, 0.0, 1.0);
266    }
267
268    #[test]
269    fn second_order_rejects_non_second_order() {
270        assert_eq!(
271            solve_linear_second_order(0.0, 1.0, 1.0, "x").unwrap_err(),
272            OdeError::NotSupported
273        );
274    }
275
276    #[test]
277    fn transport_pde_invariant_satisfies_equation() {
278        // a uₓ + b u_y = 0 with a=2,b=3 → u = F(3x − 2y). Test F=(·)² : u=(3x−2y)².
279        let PdeSolution::GeneralFunctionOf { invariant } =
280            solve_first_order_linear_pde(2.0, 3.0, "x", "y").unwrap();
281        let u = pow(invariant, 2);
282        let ux = simplify(&differentiate(&u, "x"));
283        let uy = simplify(&differentiate(&u, "y"));
284        for &(x, y) in &[(0.0, 0.0), (1.0, 2.0), (-1.0, 0.5)] {
285            let e = env(&[("x", x), ("y", y)]);
286            let r = 2.0 * ux.eval(&e).unwrap() + 3.0 * uy.eval(&e).unwrap();
287            assert!(r.abs() < 1e-7, "transport residual {r}");
288        }
289    }
290
291    #[test]
292    fn second_order_pde_classification() {
293        // Laplace uₓₓ + u_yy: A=1,B=0,C=1 → elliptic.
294        assert_eq!(classify_second_order_pde(1.0, 0.0, 1.0), PdeClass::Elliptic);
295        // Wave uₓₓ − u_yy: A=1,B=0,C=−1 → hyperbolic.
296        assert_eq!(
297            classify_second_order_pde(1.0, 0.0, -1.0),
298            PdeClass::Hyperbolic
299        );
300        // Heat-like uₓₓ (no u_yy): A=1,B=0,C=0 → parabolic.
301        assert_eq!(
302            classify_second_order_pde(1.0, 0.0, 0.0),
303            PdeClass::Parabolic
304        );
305    }
306}
qualia_core_db/specialized_libs/symbolic_ode.rs

qualia_core_db/specialized_libs/
symbolic_ode.rs
