diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 1b2b09dbe..96bad4a36 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -3145,6 +3145,23 @@ Both columns measured on a build, `v2.5.0` against this change. | `"x*sin(ln(x))/ln(x)".Integrate("x")` | `integral(x * sin(ln(x)) / ln(x), x)` | `(-1/2 * i) * Ei((2 + i) * ln(x)) + 1/2 * i * Ei((2 - i) * ln(x)) + C` | | `"cos(a+b*ln(c*x^n))^2".Integrate("x")` | `integral(cos(a + b * ln(c * x ^ n)) ^ 2, x)` | an antiderivative in `sin` and `cos` of `2 (a + b ln(c x^n))`, beside `x (c x^n)^(-1/n)` | +### Exponentials below the bar beside a power of a linear are integrated to the exponential integral + +**Answers where there were none.** `1/(x e^(2x))` was declined where `e^(-2x)/x` was answered, and so +were `1/((c + d x)(a + a tanh(e + f x)))` and its powers, Rubi's 6.3.1, and the same with `coth`, 6.4.1: +the rules for an exponential over a linear read it above the bar. Written in the exponential every +other one is a whole power of, exponentials whose denominator is then a power of that one alone are a +sum of its powers -- `1/(1 + tanh(z))` is `(1 + e^(-2z))/2` -- and each term over the linear is an +exponential integral. Where anything else stays below the bar, `1/(x (1 + e^x))`, the integral is not +of this kind and is declined as before ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"1/(x*exp(2*x))".ToEntity().Integrate("x")` | `integral(...)` | `Ei(-2 x)` | +| `"1/(x*(1 + tanh(x)))".ToEntity().Integrate("x")` | `integral(...)` | `ln(x)/2 + Ei(-2 x)/2` | +| `"1/((c + d*x)*(a + a*tanh(e + f*x)))".ToEntity().Integrate("x")` | `integral(...)` | a logarithm and an exponential integral | +| `"1/((c + d*x)^2*(a + a*coth(e + f*x))^3)".ToEntity().Integrate("x")` | `integral(...)` | powers of `c + d x` and exponential integrals | + ### An exponential or a hyperbolic function over several linears is split into partial fractions over them `e^x/(x (x + 1))` was left unevaluated, where `e^x/x` and `e^x/(x + 1)` were each answered with diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index ab8942487..30c20076f 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -5309,6 +5309,141 @@ void AddALinear(Entity linear, int power) return Functions.PartialFractions.Bare((constant * sum).InnerSimplified); } + /// + /// Exponentials of one linear below the bar, beside polynomials of x with one below + /// the bar: in w, the exponential every other is a whole power of, the exponentials + /// are a rational function of w, and where that is a sum of whole powers of w + /// -- its denominator in lowest terms a power of w alone -- each term is an + /// exponential of a linear over the polynomials, which + /// answers onto the exponential + /// integral. 1 + tanh(z) is 2 e^(2z)/(e^(2z) + 1), so + /// 1/((c + d x)(a + a tanh(e + f x))) is (1 + e^(-2(e + f x)))/(2 a (c + d x)). + /// Rubi's 6.3.1 and 6.4.1, (c + d x)^m (a + a tanh(e + f x))^n and the same with + /// coth, for a negative m. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + /// + /// The rules for an exponential over a linear read it above the bar: e^(-2x)/x was + /// answered and 1/(x e^(2x)) declined. Where the denominator keeps a factor other + /// than w after cancelling, 1/(x (1 + e^x)), the integral is neither elementary + /// nor an exponential integral, and this declines. + /// + internal static Entity? SolveAnExponentialBelowTheBarBesideAPowerOfALinear(Entity expr, Entity.Variable x) + { + bool IsAnExponential(Entity node) => node is Powf(var b, var p) && !b.ContainsNode(x) && p.ContainsNode(x); + if (!expr.Nodes.Any(node => node is Divf(_, var below) && below.Nodes.Any(IsAnExponential) + || node is Powf(var raised, Number.Integer { IsNegative: true }) && raised.Nodes.Any(IsAnExponential))) + return null; + Entity constant = Number.Integer.One; + Entity aboveX = Number.Integer.One; + Entity belowX = Number.Integer.One; + Entity exponentials = Number.Integer.One; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (factor.Nodes.Any(IsAnExponential)) + { + exponentials = underneath ? exponentials / factor : exponentials * factor; + continue; + } + // A polynomial, or a whole power of one, above or below the bar. + var (@base, power) = factor is Powf(var raised, Number.Integer whole) && whole.EInteger.CanFitInInt32() && !whole.EInteger.IsZero + ? (raised, whole.EInteger.ToInt32Checked()) : (factor, 1); + if (underneath) + power = -power; + if (!TreeAnalyzer.TryGetPolynomial(@base, x, out var read) || read.Keys.Any(degree => degree.Sign < 0) || read.Values.Any(coefficient => coefficient.ContainsNode(x))) + return null; + if (power > 0) + aboveX = aboveX == Number.Integer.One ? MathS.Pow(@base, power) : aboveX * MathS.Pow(@base, power); + else + belowX = belowX == Number.Integer.One ? MathS.Pow(@base, -power) : belowX * MathS.Pow(@base, -power); + } + // Nothing of x below the bar: a polynomial beside the exponentials is the rules' by parts. + if (belowX == Number.Integer.One || exponentials.Complexity > 160) + return null; + + // One base, and every exponent a rational multiple of the first, its offset included, + // so that w is a power of the base at that exponent and no constant is left beside it. + Entity? commonBase = null; + Entity? unit = null, unitSlope = null, unitOffset = null; + var multiples = new Dictionary(); + foreach (var node in exponentials.Nodes) + { + if (!IsAnExponential(node) || multiples.ContainsKey(node) || node is not Powf(var @base, var exponent)) + continue; + if (!TreeAnalyzer.TryGetPolyLinear(exponent, x, out var slope, out var offset) || TreeAnalyzer.IsZero(slope)) + return null; + if (commonBase is null) + { + (commonBase, unit, unitSlope, unitOffset) = (@base, exponent, slope, offset); + multiples[node] = ERational.One; + continue; + } + if (@base != commonBase + || Functions.PartialFractions.Bare((slope / unitSlope!).Simplify()).Evaled is not Number.Rational ratio || ratio.ERational.IsZero + || !IsTheZeroPolynomial((offset - ratio * unitOffset!).InnerSimplified)) + return null; + multiples[node] = ratio.ERational; + } + if (commonBase is null || unit is null) + return null; + if (commonBase != MathS.e && (commonBase.Evaled is Number.Complex and not Number.Real || commonBase.Evaled is Number.Real { IsNegative: true } || TreeAnalyzer.IsZero(commonBase))) + return null; + var numerators = EInteger.Zero; + var denominators = EInteger.One; + foreach (var multiple in multiples.Values) + { + numerators = numerators.Gcd(multiple.Numerator.Abs()); + denominators = denominators.Multiply(multiple.Denominator).Divide(denominators.Gcd(multiple.Denominator)); + } + var step = ERational.Create(numerators, denominators); + if (multiples.Values.Any(multiple => multiple.Divide(step).ToLowestTerms().Numerator.Abs().CompareTo(EInteger.FromInt32(12)) > 0)) + return null; + var w = Variable.CreateUnique(expr, "w_exp"); + var inW = exponentials.Replace(node => + multiples.TryGetValue(node, out var multiple) ? MathS.Pow(w, Number.Integer.Create(multiple.Divide(step).ToLowestTerms().Numerator)) : node); + if (inW.ContainsNode(x)) + return null; + + // A sum of whole powers of w: its denominator in lowest terms a power of w alone. + var (top, bottom) = Functions.SingleQuotient.Of(inW); + if (!TreeAnalyzer.TryGetPolynomial(bottom, w, out var below) || below.Count != 1) + { + if (!Functions.PolynomialGcd.TryCancel(top, bottom, out var cancelled, maxComplexity: 1024)) + return null; + (top, bottom) = Functions.SingleQuotient.Of(cancelled is Providedf(var inner, _) ? inner : cancelled); + if (!TreeAnalyzer.TryGetPolynomial(bottom, w, out below) || below.Count != 1) + return null; + } + var lowest = below.Keys.Single(); + var leading = below[lowest]; + if (!TreeAnalyzer.TryGetPolynomial(top, w, out var above) || above.Count > 16 || TreeAnalyzer.IsZero(leading)) + return null; + + Entity sum = Number.Integer.Zero; + foreach (var term in above) + { + if (TreeAnalyzer.IsZero(term.Value)) + continue; + var power = term.Key.Subtract(lowest); + var exponential = power.IsZero ? Number.Integer.One + : MathS.Pow(commonBase, (Number.Rational.Create(step.Multiply(ERational.FromEInteger(power))) * unit).InnerSimplified); + var question = constant * term.Value / leading * aboveX * exponential / belowX; + var integral = power.IsZero + ? Integration.ComputeIndefiniteIntegral(question, x, false) + : SolveAnExponentialOfALinearOverAPowerOfALinear(question, x) ?? SolveAnExponentialOverSeveralLinears(question, x) + ?? Integration.ComputeIndefiniteIntegral(question, x, false); + if (integral is null || integral.Nodes.Any(node => node is Integralf)) + return null; + sum += integral; + } + return sum == Number.Integer.Zero ? null : Functions.PartialFractions.Bare(sum.InnerSimplified); + } + /// /// An exponential of a linear times sines and cosines of linears, and a polynomial, over /// linears: each sine and cosine is written as exponentials, sin(c x) = (e^(i c x) - e^(-i c x))/(2i), diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 3086fd9c7..087ac67fb 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -656,6 +656,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // And either over several linears, split into partial fractions over them first, as the // trigonometric rule below splits. if ((answer = IndefiniteIntegralSolver.SolveAnExponentialOverSeveralLinears(expr, x)) is { }) return answer; + // And exponentials below the bar, where they are whole powers of one exponential and + // its powers alone are left below: `1/((c + d x)(a + a tanh(e + f x)))`. + if ((answer = IndefiniteIntegralSolver.SolveAnExponentialBelowTheBarBesideAPowerOfALinear(expr, x)) is { }) return answer; // And with sines and cosines beside the exponential, written as exponentials: each term // is then the exponential's, with a complex rate. if ((answer = IndefiniteIntegralSolver.SolveAnExponentialTimesATrigonometricOverLinears(expr, x)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/ExponentialBelowTheBarIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ExponentialBelowTheBarIntegralTest.cs new file mode 100644 index 000000000..fd36483b3 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ExponentialBelowTheBarIntegralTest.cs @@ -0,0 +1,57 @@ +// +// Copyright (c) 2019-2026 Angouri. +// AngouriMath is licensed under MIT. +// Details: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md. +// Website: https://am.angouri.org. +// + +using System; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// Exponentials of one linear below the bar, beside a power of a linear: written in the + /// exponential the others are whole powers of, they are a sum of whole powers of it, and + /// each term over the linear is an exponential integral. 1/(x e^(2x)) was declined + /// where e^(-2x)/x was answered, and so was 1/((c + d x)(a + a tanh(k + f x))), + /// Rubi's 6.3.1, which is (1 + e^(-2(k + f x)))/(2 a (c + d x)). + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class ExponentialBelowTheBarIntegralTest + { + [Theory] + [InlineData("1/(x*exp(2*x))")] + [InlineData("(exp(2*x) + 1)/(2*x*exp(2*x))")] + [InlineData("1/(x*(1 + tanh(x)))")] + [InlineData("1/((c + d*x)*(a + a*tanh(k + f*x)))")] + [InlineData("1/((c + d*x)^2*(a + a*tanh(k + f*x))^3)")] + [InlineData("1/((c + d*x)*(a + a*coth(k + f*x))^2)")] + [InlineData("1/((c + d*x)*(a - a*tanh(e + f*x)))")] + [InlineData("1/(x*(x + 1)*exp(x))")] + public void ExpandedOverTheLinear(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.DoesNotContain("NaN", integral.Stringize()); + Entity Pinned(Entity e) => e.Substitute("a", 2.3).Substitute("c", 1.3).Substitute("d", 1.7).Substitute("f", 0.9).Substitute("k", 0.4); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in new[] { -1.7, -0.9, 0.3, 0.8, 1.6, 2.9 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + if (want.IsNaN || Math.Abs((double)want.ImaginaryPart) > 1e-12) + continue; + compared++; + var got = derivative.Substitute("x", at).EvalNumerical(); + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand}"); + } + } +}