gardesk/garcalc / dcd4f0a

Browse files

Improve integral evaluation and symbolic simplification

Authored by mfwolffe <wolffemf@dukes.jmu.edu>
SHA
dcd4f0a7669096b59b5525c44c657ad2ab2b10e1
Parents
99835d3
Tree
96abcf6

2 changed files

StatusFile+-
M garcalc-cas/src/eval.rs 149 5
M garcalc-cas/src/symbolic.rs 48 0
garcalc-cas/src/eval.rsmodified
@@ -7,7 +7,7 @@ use std::collections::HashMap;
77
 use std::f64::consts::{E, PI};
88
99
 use crate::error::{CasError, Result};
10
-use crate::expr::{Expr, Rational, Sign};
10
+use crate::expr::{Expr, Rational, Sign, Symbol};
1111
 use crate::symbolic::{Differentiator, Integrator, Limits, Simplifier, Solver};
1212
1313
 /// Variable bindings for evaluation
@@ -185,10 +185,18 @@ impl Evaluator {
185185
                 upper,
186186
             } => {
187187
                 if let (Some(l), Some(u)) = (lower, upper) {
188
-                    // Definite integral - evaluate to number
188
+                    // Definite integral - try symbolic antiderivative first.
189189
                     let result = Integrator::integrate_definite(inner, var, l, u)?;
190190
                     let simplified = Simplifier::simplify(&result);
191
-                    self.eval(&simplified)
191
+                    if Self::is_unevaluated_definite_integral(&simplified) {
192
+                        // If no closed form is available, fall back to numerical quadrature.
193
+                        match self.eval_definite_integral_numeric(inner, var, l, u) {
194
+                            Ok(value) => Ok(value),
195
+                            Err(_) => Ok(simplified),
196
+                        }
197
+                    } else {
198
+                        self.eval(&simplified)
199
+                    }
192200
                 } else {
193201
                     // Indefinite integral - return symbolic result
194202
                     let result = Integrator::integrate(inner, var)?;
@@ -545,8 +553,16 @@ impl Evaluator {
545553
                 // integrate(expr, var, lower, upper)
546554
                 if let Expr::Symbol(var) = &args[1] {
547555
                     let result = Integrator::integrate_definite(&args[0], var, &args[2], &args[3])?;
548
-                    // Try to evaluate the result numerically
549
-                    self.eval(&Simplifier::simplify(&result))
556
+                    let simplified = Simplifier::simplify(&result);
557
+                    if Self::is_unevaluated_definite_integral(&simplified) {
558
+                        match self.eval_definite_integral_numeric(&args[0], var, &args[2], &args[3])
559
+                        {
560
+                            Ok(value) => Ok(value),
561
+                            Err(_) => Ok(simplified),
562
+                        }
563
+                    } else {
564
+                        self.eval(&simplified)
565
+                    }
550566
                 } else {
551567
                     Err(CasError::Type(
552568
                         "integrate requires variable as second argument".to_string(),
@@ -996,6 +1012,115 @@ impl Evaluator {
9961012
         }
9971013
     }
9981014
1015
+    fn is_unevaluated_definite_integral(expr: &Expr) -> bool {
1016
+        matches!(
1017
+            expr,
1018
+            Expr::Integral {
1019
+                lower: Some(_),
1020
+                upper: Some(_),
1021
+                ..
1022
+            }
1023
+        )
1024
+    }
1025
+
1026
+    fn eval_definite_integral_numeric(
1027
+        &self,
1028
+        body: &Expr,
1029
+        var: &Symbol,
1030
+        lower: &Expr,
1031
+        upper: &Expr,
1032
+    ) -> Result<Expr> {
1033
+        let lower_eval = self.eval(lower)?;
1034
+        let upper_eval = self.eval(upper)?;
1035
+        let mut a = self.to_f64(&lower_eval)?;
1036
+        let mut b = self.to_f64(&upper_eval)?;
1037
+
1038
+        if !a.is_finite() || !b.is_finite() {
1039
+            return Err(CasError::Type(
1040
+                "integral bounds must be finite numbers".to_string(),
1041
+            ));
1042
+        }
1043
+
1044
+        if (a - b).abs() < 1e-14 {
1045
+            return Ok(Expr::Integer(0));
1046
+        }
1047
+
1048
+        let mut sign = 1.0;
1049
+        if a > b {
1050
+            std::mem::swap(&mut a, &mut b);
1051
+            sign = -1.0;
1052
+        }
1053
+
1054
+        let mut slices = 64usize;
1055
+        let mut estimate = self.simpson_integral(body, var, a, b, slices)?;
1056
+        for _ in 0..8 {
1057
+            slices *= 2;
1058
+            let refined = self.simpson_integral(body, var, a, b, slices)?;
1059
+            if (refined - estimate).abs() <= 1e-10 * (1.0 + refined.abs()) {
1060
+                return Ok(Self::float_to_expr(sign * refined));
1061
+            }
1062
+            estimate = refined;
1063
+        }
1064
+
1065
+        Ok(Self::float_to_expr(sign * estimate))
1066
+    }
1067
+
1068
+    fn simpson_integral(
1069
+        &self,
1070
+        body: &Expr,
1071
+        var: &Symbol,
1072
+        a: f64,
1073
+        b: f64,
1074
+        slices: usize,
1075
+    ) -> Result<f64> {
1076
+        if slices == 0 || slices % 2 != 0 {
1077
+            return Err(CasError::EvaluationError(
1078
+                "simpson integration requires a positive even number of slices".to_string(),
1079
+            ));
1080
+        }
1081
+
1082
+        let h = (b - a) / slices as f64;
1083
+        let mut acc =
1084
+            self.eval_integrand_point(body, var, a)? + self.eval_integrand_point(body, var, b)?;
1085
+
1086
+        for i in 1..slices {
1087
+            let x = a + i as f64 * h;
1088
+            let fx = self.eval_integrand_point(body, var, x)?;
1089
+            if i % 2 == 0 {
1090
+                acc += 2.0 * fx;
1091
+            } else {
1092
+                acc += 4.0 * fx;
1093
+            }
1094
+        }
1095
+
1096
+        Ok(acc * h / 3.0)
1097
+    }
1098
+
1099
+    fn eval_integrand_point(&self, body: &Expr, var: &Symbol, x: f64) -> Result<f64> {
1100
+        let substituted = Simplifier::substitute(body, var, &Expr::Float(x));
1101
+        let evaluated = self.eval(&substituted)?;
1102
+        let value = self.to_f64(&evaluated)?;
1103
+        if value.is_finite() {
1104
+            Ok(value)
1105
+        } else {
1106
+            Err(CasError::EvaluationError(format!(
1107
+                "integrand is not finite at {x}"
1108
+            )))
1109
+        }
1110
+    }
1111
+
1112
+    fn float_to_expr(x: f64) -> Expr {
1113
+        if !x.is_finite() {
1114
+            return Expr::Float(x);
1115
+        }
1116
+        let rounded = x.round();
1117
+        if (x - rounded).abs() < 1e-10 && rounded >= i64::MIN as f64 && rounded <= i64::MAX as f64 {
1118
+            Expr::Integer(rounded as i64)
1119
+        } else {
1120
+            Expr::Float(x)
1121
+        }
1122
+    }
1123
+
9991124
     fn eval_sum(
10001125
         &self,
10011126
         body: &Expr,
@@ -1647,6 +1772,25 @@ mod tests {
16471772
         assert_eq!(eval("product(n, n, 1, 4)").unwrap(), Expr::Integer(24));
16481773
     }
16491774
1775
+    #[test]
1776
+    fn test_definite_integral_numeric_fallback_for_non_elementary_antiderivative() {
1777
+        let result = eval_to_f64("integrate(x^(x+2), x, 0, 1)");
1778
+        assert!((result - 0.2781176122).abs() < 1e-8);
1779
+    }
1780
+
1781
+    #[test]
1782
+    fn test_definite_integral_with_symbolic_bounds_stays_symbolic() {
1783
+        let symbolic = eval("integrate(x^(x+2), x, 0, n)").unwrap();
1784
+        assert!(matches!(
1785
+            symbolic,
1786
+            Expr::Integral {
1787
+                lower: Some(_),
1788
+                upper: Some(_),
1789
+                ..
1790
+            }
1791
+        ));
1792
+    }
1793
+
16501794
     #[test]
16511795
     fn test_empty_discrete_range() {
16521796
         assert_eq!(eval("sum(n, n, 5, 1)").unwrap(), Expr::Integer(0));
garcalc-cas/src/symbolic.rsmodified
@@ -4,6 +4,7 @@
44
55
 use crate::error::{CasError, Result};
66
 use crate::expr::{Expr, Symbol};
7
+use std::f64::consts::E;
78
89
 /// Symbolic differentiator
910
 pub struct Differentiator;
@@ -329,6 +330,25 @@ impl Integrator {
329330
         // Get antiderivative
330331
         let antideriv = Self::integrate(expr, var)?;
331332
333
+        // If we could not find an antiderivative, keep the definite integral
334
+        // unevaluated instead of substituting and incorrectly collapsing to zero.
335
+        if matches!(
336
+            &antideriv,
337
+            Expr::Integral {
338
+                expr: inner,
339
+                var: v,
340
+                lower: None,
341
+                upper: None
342
+            } if **inner == *expr && v == var
343
+        ) {
344
+            return Ok(Expr::Integral {
345
+                expr: Box::new(expr.clone()),
346
+                var: var.clone(),
347
+                lower: Some(Box::new(lower.clone())),
348
+                upper: Some(Box::new(upper.clone())),
349
+            });
350
+        }
351
+
332352
         // F(upper) - F(lower)
333353
         let f_upper = Simplifier::substitute(&antideriv, var, upper);
334354
         let f_lower = Simplifier::substitute(&antideriv, var, lower);
@@ -1229,6 +1249,19 @@ impl Simplifier {
12291249
             Expr::Func(name, args) => {
12301250
                 let sargs: Vec<_> = args.iter().map(Self::simplify_impl).collect();
12311251
1252
+                if name == "ln" && sargs.len() == 1 {
1253
+                    match &sargs[0] {
1254
+                        Expr::Integer(1) => return Expr::Integer(0),
1255
+                        Expr::Rational(r) if r.den != 0 && r.num == r.den => {
1256
+                            return Expr::Integer(0);
1257
+                        }
1258
+                        Expr::Float(x) if (*x - 1.0).abs() < 1e-12 => return Expr::Integer(0),
1259
+                        Expr::Symbol(sym) if sym.as_str() == "e" => return Expr::Integer(1),
1260
+                        Expr::Float(x) if (*x - E).abs() < 1e-12 => return Expr::Integer(1),
1261
+                        _ => {}
1262
+                    }
1263
+                }
1264
+
12321265
                 if name == "factorial" && sargs.len() == 1 {
12331266
                     if let Expr::Integer(n) = sargs[0] {
12341267
                         if n >= 0 {
@@ -1979,6 +2012,21 @@ mod tests {
19792012
         assert_eq!(simplified, Expr::Integer(1));
19802013
     }
19812014
2015
+    #[test]
2016
+    fn test_simplify_ln_of_e() {
2017
+        let expr = Expr::func("ln", vec![Expr::symbol("e")]);
2018
+        let simplified = Simplifier::simplify(&expr);
2019
+        assert_eq!(simplified, Expr::Integer(1));
2020
+    }
2021
+
2022
+    #[test]
2023
+    fn test_diff_e_to_x_simplifies_to_e_to_x() {
2024
+        let expr = Expr::pow(Expr::symbol("e"), Expr::symbol("x"));
2025
+        let result = Differentiator::diff(&expr, &Symbol::new("x")).unwrap();
2026
+        let simplified = Simplifier::simplify(&result);
2027
+        assert_eq!(simplified, Expr::pow(Expr::symbol("e"), Expr::symbol("x")));
2028
+    }
2029
+
19822030
     #[test]
19832031
     fn test_simplify_product_residual_factorial_inverse() {
19842032
         let expr = Expr::mul(vec![