qualia_core_db/specialized_libs/
symbolic_limits.rs1use super::symbolic_algebra::{differentiate, simplify, Expr};
7
8const TOL: f64 = 1e-9;
9const MAX_HOPITAL: usize = 8;
10
11fn eval_at(e: &Expr, x: &str, v: f64) -> Option<f64> {
12 let mut env = std::collections::HashMap::new();
13 env.insert(x.to_string(), v);
14 e.eval(&env).filter(|r| r.is_finite())
15}
16
17pub fn limit(f: &Expr, x: &str, a: f64) -> Option<f64> {
21 if let Some(v) = eval_at(f, x, a) {
23 return Some(v);
24 }
25 if let Expr::Div(num, den) = f {
27 let (mut n, mut d) = ((**num).clone(), (**den).clone());
28 for _ in 0..MAX_HOPITAL {
29 let nv = eval_at(&n, x, a);
30 let dv = eval_at(&d, x, a);
31 match (nv, dv) {
32 (Some(nn), Some(dd)) => {
33 if dd.abs() > TOL {
34 return Some(nn / dd); }
36 if nn.abs() > TOL {
37 return None; }
39 n = simplify(&differentiate(&n, x));
41 d = simplify(&differentiate(&d, x));
42 }
43 _ => return None,
44 }
45 }
46 }
47 None
48}
49
50pub fn limit_at_infinity(f: &Expr, x: &str) -> Option<f64> {
53 let mut prev = eval_at(f, x, 1e3)?;
54 for &t in &[1e4, 1e5, 1e6, 1e7] {
55 let cur = eval_at(f, x, t)?;
56 if (cur - prev).abs() < 1e-6 {
57 return Some(cur);
58 }
59 prev = cur;
60 }
61 Some(prev)
63}
64
65#[cfg(test)]
66mod tests {
67 use super::super::symbolic_algebra::{add, c, div, mul, pow, sub, var};
68 use super::*;
69
70 #[test]
71 fn lhopital_on_zero_over_zero() {
72 let f = div(sub(pow(var("x"), 2), c(1.0)), sub(var("x"), c(1.0)));
74 assert!((limit(&f, "x", 1.0).unwrap() - 2.0).abs() < 1e-7);
75 let g = div(sub(pow(var("x"), 2), c(4.0)), sub(var("x"), c(2.0)));
77 assert!((limit(&g, "x", 2.0).unwrap() - 4.0).abs() < 1e-7);
78 }
79
80 #[test]
81 fn direct_substitution_when_defined() {
82 let f = add(pow(var("x"), 2), c(1.0));
84 assert!((limit(&f, "x", 3.0).unwrap() - 10.0).abs() < 1e-9);
85 }
86
87 #[test]
88 fn divergent_fails_closed() {
89 let f = div(c(1.0), var("x"));
91 assert!(limit(&f, "x", 0.0).is_none());
92 }
93
94 #[test]
95 fn rational_limit_at_infinity() {
96 let f = div(
98 add(mul(c(2.0), pow(var("x"), 2)), c(3.0)),
99 sub(pow(var("x"), 2), c(1.0)),
100 );
101 assert!((limit_at_infinity(&f, "x").unwrap() - 2.0).abs() < 1e-3);
102 }
103}