gardesk/garcalc / 839c390

Browse files

add linear substitution and u-substitution to integrator

Authored by mfwolffe <wolffemf@dukes.jmu.edu>
SHA
839c3903f17e4594afa8c22e3095fefc8a6b5f5a
Parents
25b1d04
Tree
9058590

1 changed file

StatusFile+-
M garcalc-cas/src/symbolic.rs 281 1
garcalc-cas/src/symbolic.rsmodified
@@ -445,6 +445,46 @@ impl Integrator {
445445
                     ));
446446
                 }
447447
 
448
+                // Linear substitution: ∫ (ax+b)^n dx = (ax+b)^(n+1) / (a*(n+1))
449
+                if !exp.contains_var(var) {
450
+                    if let Some((a, _b)) = Self::extract_linear(base, var) {
451
+                        if exp.is_negative_one() {
452
+                            // ∫ (ax+b)^(-1) dx = ln|ax+b| / a
453
+                            return Ok(Expr::div(
454
+                                Expr::func("ln", vec![Expr::func("abs", vec![(**base).clone()])]),
455
+                                a,
456
+                            ));
457
+                        }
458
+                        let n_plus_1 = Expr::add(vec![(**exp).clone(), Expr::Integer(1)]);
459
+                        return Ok(Expr::div(
460
+                            Expr::pow((**base).clone(), n_plus_1.clone()),
461
+                            Expr::mul(vec![a, n_plus_1]),
462
+                        ));
463
+                    }
464
+                }
465
+
466
+                // Linear substitution for exponentials: ∫ e^(ax+b) dx = e^(ax+b) / a
467
+                if let Expr::Symbol(s) = &**base {
468
+                    if s.as_str() == "e" {
469
+                        if let Some((a, _b)) = Self::extract_linear(exp, var) {
470
+                            return Ok(Expr::div(
471
+                                Expr::pow(Expr::symbol("e"), (**exp).clone()),
472
+                                a,
473
+                            ));
474
+                        }
475
+                    }
476
+                }
477
+
478
+                // ∫ a^(cx+d) dx = a^(cx+d) / (c * ln(a))
479
+                if !base.contains_var(var) {
480
+                    if let Some((c, _d)) = Self::extract_linear(exp, var) {
481
+                        return Ok(Expr::div(
482
+                            expr.clone(),
483
+                            Expr::mul(vec![c, Expr::func("ln", vec![(**base).clone()])]),
484
+                        ));
485
+                    }
486
+                }
487
+
448488
                 // Return unevaluated
449489
                 Ok(Expr::Integral {
450490
                     expr: Box::new(expr.clone()),
@@ -473,6 +513,18 @@ impl Integrator {
473513
             return Self::integrate(factors[0], var);
474514
         }
475515
 
516
+        // U-substitution: look for f(g(x)) * g'(x) patterns
517
+        // For each pair of factors, check if one is the derivative of the inner
518
+        // function of the other (up to a constant multiple).
519
+        if factors.len() == 2 {
520
+            // Try both orderings: factor[0] as f(g(x)), factor[1] as g'(x) and vice versa
521
+            for (fi, fj) in [(0, 1), (1, 0)] {
522
+                if let Some(result) = Self::try_u_substitution(factors[fi], factors[fj], var) {
523
+                    return Ok(result);
524
+                }
525
+            }
526
+        }
527
+
476528
         // Return unevaluated
477529
         Ok(Expr::Integral {
478530
             expr: Box::new(Expr::mul(factors.iter().map(|f| (*f).clone()).collect())),
@@ -482,6 +534,223 @@ impl Integrator {
482534
         })
483535
     }
484536
 
537
+    /// Try u-substitution: given f_expr (containing g(x)) and candidate g'(x),
538
+    /// check if candidate equals g'(x) up to a constant, and if so compute the integral.
539
+    fn try_u_substitution(f_expr: &Expr, candidate_deriv: &Expr, var: &Symbol) -> Option<Expr> {
540
+        // Extract inner function g(x) from f_expr
541
+        let inner = Self::extract_inner_function(f_expr, var)?;
542
+
543
+        // Compute g'(x)
544
+        let g_prime = Differentiator::diff(&inner, var).ok()?;
545
+        let g_prime_simplified = Simplifier::simplify(&g_prime);
546
+
547
+        // Check if candidate_deriv = c * g'(x) for some constant c
548
+        let constant_multiple = Self::is_constant_multiple(candidate_deriv, &g_prime_simplified, var)?;
549
+
550
+        // Build f(u) by replacing g(x) with u in f_expr
551
+        let u = Symbol::new("_u_sub_");
552
+        let f_of_u = Self::replace_subexpr(f_expr, &inner, &Expr::Symbol(u.clone()));
553
+
554
+        // If substitution didn't fully replace var, this pattern doesn't work
555
+        if f_of_u.contains_var(var) {
556
+            return None;
557
+        }
558
+
559
+        let antideriv = Self::integrate(&f_of_u, &u).ok()?;
560
+        // Check it's not unevaluated
561
+        if matches!(&antideriv, Expr::Integral { .. }) {
562
+            return None;
563
+        }
564
+
565
+        // Substitute g(x) back for u
566
+        let result = Simplifier::substitute(&antideriv, &u, &inner);
567
+
568
+        // Multiply by constant_multiple
569
+        if constant_multiple.is_one() {
570
+            Some(result)
571
+        } else {
572
+            Some(Expr::mul(vec![constant_multiple, result]))
573
+        }
574
+    }
575
+
576
+    /// Replace occurrences of `target` subexpression with `replacement` in `expr`.
577
+    fn replace_subexpr(expr: &Expr, target: &Expr, replacement: &Expr) -> Expr {
578
+        if Self::exprs_structurally_equal(expr, target) {
579
+            return replacement.clone();
580
+        }
581
+        match expr {
582
+            Expr::Neg(e) => Expr::neg(Self::replace_subexpr(e, target, replacement)),
583
+            Expr::Add(terms) => Expr::add(
584
+                terms.iter().map(|t| Self::replace_subexpr(t, target, replacement)).collect(),
585
+            ),
586
+            Expr::Mul(factors) => Expr::mul(
587
+                factors.iter().map(|f| Self::replace_subexpr(f, target, replacement)).collect(),
588
+            ),
589
+            Expr::Pow(base, exp) => Expr::pow(
590
+                Self::replace_subexpr(base, target, replacement),
591
+                Self::replace_subexpr(exp, target, replacement),
592
+            ),
593
+            Expr::Func(name, args) => Expr::func(
594
+                name,
595
+                args.iter().map(|a| Self::replace_subexpr(a, target, replacement)).collect(),
596
+            ),
597
+            _ => expr.clone(),
598
+        }
599
+    }
600
+
601
+    /// Extract the inner function g(x) from an expression that looks like f(g(x)).
602
+    /// Returns the innermost composite argument that contains `var`.
603
+    fn extract_inner_function(expr: &Expr, var: &Symbol) -> Option<Expr> {
604
+        match expr {
605
+            // f(g(x)) — the inner function is g(x)
606
+            Expr::Func(_name, args) if args.len() == 1 => {
607
+                let arg = &args[0];
608
+                if *arg != Expr::Symbol(var.clone()) && arg.contains_var(var) {
609
+                    Some(arg.clone())
610
+                } else {
611
+                    None
612
+                }
613
+            }
614
+            // (g(x))^n — the inner function is g(x) if g(x) is not just x
615
+            Expr::Pow(base, exp) if !exp.contains_var(var) && base.contains_var(var) => {
616
+                if **base != Expr::Symbol(var.clone()) {
617
+                    Some((**base).clone())
618
+                } else {
619
+                    None
620
+                }
621
+            }
622
+            _ => None,
623
+        }
624
+    }
625
+
626
+    /// Check if `expr` is a constant multiple of `target` w.r.t. `var`.
627
+    /// Returns the constant factor if so.
628
+    fn is_constant_multiple(expr: &Expr, target: &Expr, var: &Symbol) -> Option<Expr> {
629
+        // Simplify both for comparison
630
+        let expr_s = Simplifier::simplify(expr);
631
+        let target_s = Simplifier::simplify(target);
632
+
633
+        // Direct equality: multiple is 1
634
+        if Self::exprs_structurally_equal(&expr_s, &target_s) {
635
+            return Some(Expr::Integer(1));
636
+        }
637
+
638
+        // Check if expr = c * target where c is constant
639
+        // Try dividing: if expr/target simplifies to a constant, that's our multiple
640
+        let ratio = Expr::div(expr_s.clone(), target_s.clone());
641
+        let ratio_simplified = Simplifier::simplify(&ratio);
642
+        if !ratio_simplified.contains_var(var) {
643
+            return Some(ratio_simplified);
644
+        }
645
+
646
+        None
647
+    }
648
+
649
+    /// Simple structural equality check (after simplification)
650
+    fn exprs_structurally_equal(a: &Expr, b: &Expr) -> bool {
651
+        match (a, b) {
652
+            (Expr::Integer(x), Expr::Integer(y)) => x == y,
653
+            (Expr::Float(x), Expr::Float(y)) => (x - y).abs() < 1e-12,
654
+            (Expr::Symbol(x), Expr::Symbol(y)) => x == y,
655
+            (Expr::Neg(x), Expr::Neg(y)) => Self::exprs_structurally_equal(x, y),
656
+            (Expr::Add(xs), Expr::Add(ys)) | (Expr::Mul(xs), Expr::Mul(ys)) => {
657
+                xs.len() == ys.len()
658
+                    && xs
659
+                        .iter()
660
+                        .zip(ys.iter())
661
+                        .all(|(x, y)| Self::exprs_structurally_equal(x, y))
662
+            }
663
+            (Expr::Pow(xb, xe), Expr::Pow(yb, ye)) => {
664
+                Self::exprs_structurally_equal(xb, yb) && Self::exprs_structurally_equal(xe, ye)
665
+            }
666
+            (Expr::Func(xn, xa), Expr::Func(yn, ya)) => {
667
+                xn == yn
668
+                    && xa.len() == ya.len()
669
+                    && xa
670
+                        .iter()
671
+                        .zip(ya.iter())
672
+                        .all(|(x, y)| Self::exprs_structurally_equal(x, y))
673
+            }
674
+            (Expr::Rational(x), Expr::Rational(y)) => x == y,
675
+            _ => false,
676
+        }
677
+    }
678
+
679
+    /// Extract the linear form ax + b from an expression, returning (a, b).
680
+    /// Returns None if the expression is not linear in `var`.
681
+    fn extract_linear(expr: &Expr, var: &Symbol) -> Option<(Expr, Expr)> {
682
+        match expr {
683
+            // Just x → (1, 0)
684
+            Expr::Symbol(s) if s == var => Some((Expr::Integer(1), Expr::Integer(0))),
685
+
686
+            // a*x or x*a → (a, 0)
687
+            Expr::Mul(factors) => {
688
+                let (constants, var_factors): (Vec<_>, Vec<_>) =
689
+                    factors.iter().partition(|f| !f.contains_var(var));
690
+                // Need exactly one var factor that is x
691
+                if var_factors.len() == 1 && *var_factors[0] == Expr::Symbol(var.clone()) {
692
+                    let a = if constants.is_empty() {
693
+                        Expr::Integer(1)
694
+                    } else {
695
+                        Expr::mul(constants.into_iter().cloned().collect())
696
+                    };
697
+                    Some((a, Expr::Integer(0)))
698
+                } else {
699
+                    None
700
+                }
701
+            }
702
+
703
+            // ax + b or b + ax
704
+            Expr::Add(terms) => {
705
+                let mut a = Expr::Integer(0);
706
+                let mut b_parts = Vec::new();
707
+                for term in terms {
708
+                    if !term.contains_var(var) {
709
+                        b_parts.push(term.clone());
710
+                    } else if *term == Expr::Symbol(var.clone()) {
711
+                        a = Expr::add(vec![a, Expr::Integer(1)]);
712
+                    } else if let Expr::Mul(factors) = term {
713
+                        let (constants, var_factors): (Vec<_>, Vec<_>) =
714
+                            factors.iter().partition(|f| !f.contains_var(var));
715
+                        if var_factors.len() == 1 && *var_factors[0] == Expr::Symbol(var.clone()) {
716
+                            let coeff = if constants.is_empty() {
717
+                                Expr::Integer(1)
718
+                            } else {
719
+                                Expr::mul(constants.into_iter().cloned().collect())
720
+                            };
721
+                            a = Expr::add(vec![a, coeff]);
722
+                        } else {
723
+                            return None; // Non-linear term
724
+                        }
725
+                    } else {
726
+                        return None; // Non-linear term containing var
727
+                    }
728
+                }
729
+                let a = Simplifier::simplify(&a);
730
+                if a.is_zero() {
731
+                    return None; // No linear term
732
+                }
733
+                let b = if b_parts.is_empty() {
734
+                    Expr::Integer(0)
735
+                } else {
736
+                    Expr::add(b_parts)
737
+                };
738
+                Some((a, b))
739
+            }
740
+
741
+            // -x → (-1, 0)
742
+            Expr::Neg(inner) => {
743
+                if **inner == Expr::Symbol(var.clone()) {
744
+                    Some((Expr::Integer(-1), Expr::Integer(0)))
745
+                } else {
746
+                    None
747
+                }
748
+            }
749
+
750
+            _ => None,
751
+        }
752
+    }
753
+
485754
     fn integrate_func(name: &str, args: &[Expr], var: &Symbol) -> Result<Expr> {
486755
         if args.is_empty() {
487756
             return Ok(Expr::Integral {
@@ -494,8 +763,19 @@ impl Integrator {
494763
 
495764
         let arg = &args[0];
496765
 
497
-        // Only handle simple case where arg = var
766
+        // Try linear substitution: ∫ f(ax+b) dx = (1/a)*F(ax+b)
498767
         if *arg != Expr::Symbol(var.clone()) {
768
+            if let Some((a, _b)) = Self::extract_linear(arg, var) {
769
+                // Integrate as if arg = var, then divide by the linear coefficient
770
+                let simple_var = Expr::Symbol(var.clone());
771
+                let simple_integral = Self::integrate_func(name, &[simple_var], var)?;
772
+                // Check we actually got an antiderivative (not unevaluated)
773
+                if !matches!(&simple_integral, Expr::Integral { .. }) {
774
+                    // Substitute ax+b for x in the result, then divide by a
775
+                    let substituted = Simplifier::substitute(&simple_integral, var, arg);
776
+                    return Ok(Expr::div(substituted, a));
777
+                }
778
+            }
499779
             return Ok(Expr::Integral {
500780
                 expr: Box::new(Expr::func(name, args.to_vec())),
501781
                 var: var.clone(),