diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index edd686f29..90cb47a9c 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -154,6 +154,21 @@ integrand is then `x^(2n + 1) (c + d x)^n (2 c + 3 d x)`, the derivative of `x^( | `"x^(2*n)*(c+d*x)^n*(2*c*x+3*d*x^2)".ToEntity().Integrate("x")` | `integral(...)` | 61 characters | | `"x^(n-1)*(c+d*x^2)^n*(c*x+3*d*x^3)".ToEntity().Integrate("x")` | `integral(...)` | 65 characters | +### A power of a quotient of linears over a half-odd power of their product is integrated + +**Answers where there were none.** `e^(n acoth(a x))/(c - a^2 c x^2)^(3/2)` was declined, and so were its `^(5/2)` to +`^(9/2)` and the same beside `x`, `x^2` and `x^3`, Rubi's 7.4.2. The parser folds `e^(n acoth(a x))` into +`((a x + 1)/(a x - 1))^(n/2)`, a power `T^alpha` of a quotient of linears, and over a half-odd power of a multiple `Q` of +their product the antiderivative is `T^alpha P(x) Q^(1/2 - K)` for a polynomial `P` whose coefficients a three-term +recurrence gives; where the recurrence leaves an equation it cannot meet, the antiderivative is not of that shape and +it is still declined ([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"e^(n*acoth(a*x))/(c-a^2*c*x^2)^(3/2)".ToEntity().Integrate("x")` | `integral(...)` | 140 characters | +| `"e^(n*acoth(a*x))*x/(c-a^2*c*x^2)^(3/2)".ToEntity().Integrate("x")` | `integral(...)` | 155 characters | +| `"e^(n*acoth(a*x))*x^3/(c-a^2*c*x^2)^(5/2)".ToEntity().Integrate("x")` | `integral(...)` | 397 characters | + ### A power of a quotient of linears beside whole powers of them is integrated under the quotient **Answers where there were none.** `e^(n acoth(a x))/(c - a^2 c x^2)^2` was declined, and so were its `^3` and `^4`, Rubi's diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index c28cedfb6..d07e5a3f1 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -1002,6 +1002,109 @@ private static bool IsADerivativeOfAFactorBelowTheBar(Entity expr, Entity.Variab } } + /// + /// x^j T^alpha/Q^(K + 1/2), a power of a quotient of linears T = (p + q x)/(r + s x) + /// over a half-odd power of a multiple of their product, Q = kappa (p + q x)(r + s x), + /// closed by Liouville's ansatz F = T^alpha P(x) Q^(1/2 - K). With D = q r - p s, + /// T'/T = D/((p + q x)(r + s x)), so F' = T^alpha Q^(-K - 1/2) (alpha kappa D P + Q P' - (K - 1/2) Q' P), + /// and asking the bracket to be x^j is, with Q = A x^2 + B x + C, the three-term + /// system A (m - 2K) p_(m-1) + (alpha kappa D + B (m - K + 1/2)) p_m + C (m + 1) p_(m+1) = [m = j], + /// solved from p_0 up, the last equation fixing p_0 and the one past it checked. + /// e^(n acoth(a x)) arrives folded as ((a x + 1)/(a x - 1))^(n/2): Rubi's 7.4.2 + /// e^(n acoth(a x))/(c - a^2 c x^2)^(3/2) is e^(n acoth(a x)) (a x - n)/(a c (1 - n^2) sqrt(c - a^2 c x^2)). + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + /// The principal powers are differentiated as such wherever T and Q are not 0, + /// on a cut as off it, so the answer holds where the integrand is defined. Where the check + /// past the last equation fails, the antiderivative is not of this shape -- a hypergeometric + /// one -- and the rule declines. Only for a power of T that is not whole, and + /// Q without a root at 0, which the system divides by. + /// + internal static Entity? SolveAPowerOfAQuotientOfLinearsOverAHalfOddPowerOfTheirProduct(Entity expr, Entity.Variable x) + { + Powf? power = null; + Entity constant = Number.Integer.One; + var j = 0; + Entity? quadratic = null; + ERational? quadraticPower = null; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (!underneath && power is null && factor is Powf(Divf(var above, var below), var alpha) + && !alpha.ContainsNode(x) && alpha.Evaled is not Number.Integer && above.ContainsNode(x) && below.ContainsNode(x)) + { + power = (Powf)factor; + continue; + } + if (!underneath && APowerOfTheVariable(factor, x) is { } k && k > 0) + { + j += k; + continue; + } + if (quadratic is null && factor is Powf(var q, Number.Rational r) && r.ERational.Denominator.Equals(EInteger.FromInt32(2))) + { + quadratic = q; + quadraticPower = underneath ? r.ERational.Negate() : r.ERational; + continue; + } + return null; + } + if (power is null || quadratic is null || quadraticPower is null || j > 12 || power.Base is not Divf(var top, var bottom) + || !TreeAnalyzer.TryGetPolyLinear(top, x, out var qs, out var ps) || TreeAnalyzer.IsZero(qs) + || !TreeAnalyzer.TryGetPolyLinear(bottom, x, out var ss, out var rs) || TreeAnalyzer.IsZero(ss) + || !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var a2, out var b2, out var c2) || TreeAnalyzer.IsZero(a2)) + return null; + // Q a multiple of the linears' product, and its power -(K + 1/2). + var kappa = Functions.PartialFractions.Bare((a2 / (qs * ss)).Simplify()); + if (!VanishesIdentically(b2 - kappa * (ps * ss + qs * rs)) || !VanishesIdentically(c2 - kappa * ps * rs) || VanishesIdentically(c2)) + return null; + var half = quadraticPower.Negate().Subtract(ERational.FromInt32(1).Divide(ERational.FromInt32(2))).ToLowestTerms(); + if (!half.Denominator.Equals(EInteger.One) || !half.Numerator.CanFitInInt32() || half.Numerator.Abs().CompareTo(EInteger.FromInt32(12)) > 0) + return null; + var bigK = half.Numerator.ToInt32Checked(); + var d = System.Math.Max(2 * bigK - 1, j - 1); + if (d < 0 || d > 12) + return null; + var determinant = qs * rs - ps * ss; + var alphaKappaD = (power.Exponent * kappa * determinant).InnerSimplified; + Entity Beta(int m) => alphaKappaD + b2 * (Number.Integer.Create(m - bigK) + Number.Rational.Create(1, 2)); + Entity Rhs(int m) => m == j ? Number.Integer.One : Number.Integer.Zero; + // p_m = u_m + v_m lambda, lambda = p_0. + var u = new Entity[d + 2]; + var v = new Entity[d + 2]; + (u[0], v[0]) = (Number.Integer.Zero, Number.Integer.One); + Entity Lower(Entity[] w, int m) => m - 1 < 0 ? Number.Integer.Zero : a2 * (m - 2 * bigK) * w[m - 1]; + for (var m = 0; m < d; m++) + { + var divisor = c2 * (m + 1); + u[m + 1] = Functions.PartialFractions.Bare(((Rhs(m) - Lower(u, m) - Beta(m) * u[m]) / divisor).Simplify()); + v[m + 1] = Functions.PartialFractions.Bare(((Number.Integer.Zero - Lower(v, m) - Beta(m) * v[m]) / divisor).Simplify()); + } + // The degree-d equation, with p_(d+1) = 0, fixes lambda. + var constantPart = (Lower(u, d) + Beta(d) * u[d]).InnerSimplified; + var lambdaPart = (Lower(v, d) + Beta(d) * v[d]).InnerSimplified; + if (VanishesIdentically(lambdaPart)) + return null; + var lambda = Functions.PartialFractions.Bare(((Rhs(d) - constantPart) / lambdaPart).Simplify()); + var coefficients = Enumerable.Range(0, d + 1).Select(m => Functions.PartialFractions.Bare((u[m] + v[m] * lambda).Simplify())).ToArray(); + // The one past it: A (d + 1 - 2K) p_d = [d + 1 = j]. + if (!VanishesIdentically(a2 * (d + 1 - 2 * bigK) * coefficients[d] - Rhs(d + 1))) + return null; + Entity polynomial = Number.Integer.Zero; + for (var m = 0; m <= d; m++) + if (!VanishesIdentically(coefficients[m])) + polynomial += coefficients[m] * (m == 0 ? Number.Integer.One : m == 1 ? x : MathS.Pow(x, m)); + if (polynomial == Number.Integer.Zero) + return null; + var raised = Number.Rational.Create(1 - 2 * bigK, 2); + return constant * power * polynomial * MathS.Pow(quadratic, raised); + } + /// /// A quotient by the square of a base, by parts against the base's reciprocal, where the /// base's derivative makes the rest exact: N/(M v^2) is g v'/v^2 with diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 278c548e5..e43b5088c 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -791,6 +791,8 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) if ((answer = IndefiniteIntegralSolver.SolveByPartsAgainstALogarithmOfARadicalSum(expr, x)) is { }) return answer; // A function of tanh(L) or coth(L) alone that is not rational in it, under u = tanh(L). if ((answer = IndefiniteIntegralSolver.SolveAFunctionOfTheHyperbolicTangent(expr, x, integrateByParts)) is { }) return answer; + // A power of a quotient of linears over a half-odd power of their product, by Liouville's ansatz. + if ((answer = IndefiniteIntegralSolver.SolveAPowerOfAQuotientOfLinearsOverAHalfOddPowerOfTheirProduct(expr, x)) is { }) return answer; // A power of an exponential with a positive base is the exponential of the product, // exactly, and only that spelling is one the exponential rules read. if ((answer = IndefiniteIntegralSolver.SolveByFlatteningAPowerOfAnExponential(expr, x, integrateByParts)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/APowerOfAQuotientOfLinearsOverAHalfOddPowerOfTheirProductIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfAQuotientOfLinearsOverAHalfOddPowerOfTheirProductIntegralTest.cs new file mode 100644 index 000000000..8915f309e --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerOfAQuotientOfLinearsOverAHalfOddPowerOfTheirProductIntegralTest.cs @@ -0,0 +1,66 @@ +// +// 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 +{ + /// + /// x^j e^(n acoth(a x))/(c - a^2 c x^2)^(K + 1/2), closed by Liouville's ansatz + /// e^(n acoth(a x)) P(x)/(c - a^2 c x^2)^(K - 1/2): the first is Rubi's + /// e^(n acoth(a x)) (a x - n)/(a c (1 - n^2) sqrt(c - a^2 c x^2)). They were declined. + /// Rubi's 7.4.2. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APowerOfAQuotientOfLinearsOverAHalfOddPowerOfTheirProductIntegralTest + { + [Fact] + public void RubisAnswer() + { + var integral = "e^(n*acoth(a*x))/(c - a^2*c*x^2)^(3/2)".ToEntity().Integrate("x").Substitute("C", 0); + var rubi = "e^(n*acoth(a*x))*(a*x - n)/(a*c*(1 - n^2)*sqrt(c - a^2*c*x^2))".ToEntity(); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("c", 0.7).Substitute("n", 0.4); + foreach (var at in new[] { -2.9, -1.4, 1.3, 2.9 }) + { + var difference = (Pinned(integral) - Pinned(rubi)).Substitute("x", at).EvalNumerical(); + Assert.True(Math.Abs((double)difference.RealPart) + Math.Abs((double)difference.ImaginaryPart) < 1e-9, $"{difference} at {at}"); + } + } + + [Theory] + [InlineData("e^(n*acoth(a*x))/(c - a^2*c*x^2)^(5/2)")] + [InlineData("e^(n*acoth(a*x))*x/(c - a^2*c*x^2)^(3/2)")] + [InlineData("e^(n*acoth(a*x))*x^3/(c - a^2*c*x^2)^(5/2)")] + [InlineData("e^(n*atanh(a*x))/(c - a^2*c*x^2)^(3/2)")] + public void ByTheAnsatz(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("c", 0.7).Substitute("n", 0.4); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + // Inside and outside the roots at -1/1.3 and 1/1.3. + foreach (var at in new[] { -2.9, -1.4, -0.4, 0.3, 0.6, 1.3, 2.9 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + if (want.IsNaN) + continue; + compared++; + Assert.True(Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart) + < 1e-9 * Math.Max(1, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)), + $"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}"); + } + } +}