diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index a7c295ac9..da2ea9c7c 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -213,6 +213,23 @@ is a polynomial in `t` times a power of it ([#233](https://github.com/asc-commun | `"e^(n*acoth(a*x))/(c-a^2*c*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | three powers of `(a x + 1)/(a x - 1)`, `n/2 - 1` to `n/2 + 1` | | `"e^(n*acoth(a*x))/(c-a^2*c*x^2)^4".ToEntity().Integrate("x")` | `integral(...)` | seven powers of it, `n/2 - 3` to `n/2 + 3` | +### A power of x beside powers of a quotient of linears and of their product is integrated as one term + +**Answers where there were none.** `e^(n artanh(a x)) x^2 (c - a^2 c x^2)^(-1 - n^2/2)` was declined, and so were +`e^(6 artanh(a x)) x^2/(c - a^2 c x^2)^19`, `e^(5 artanh(a x)) x^2/(c - a^2 c x^2)^(27/2)` and its reciprocal twin, and +`e^(artanh(a x)/2) x^2/(c - a^2 c x^2)^(9/8)`, Rubi's 7.3.6. With `T = (1 + a x)/(1 - a x)` and `Q = c - a^2 c x^2`, +each is `x^j T^alpha Q^beta`, and its antiderivative is `T^alpha Q^(beta + 1)` times a polynomial of degree `j - 1`, +found by one linear solve whose last equation is checked; where it fails the rule declines +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"e^(n*atanh(a*x))*x^2*(c-a^2*c*x^2)^(-1-1/2*n^2)".ToEntity().Integrate("x")` | `integral(...)` | 167 characters | +| `"e^(6*atanh(a*x))*x^2/(c-a^2*c*x^2)^19".ToEntity().Integrate("x")` | `integral(...)` | `T^3 Q^(-18) (6 a x - 1)/(210 a^3 c)`, 119 characters | +| `"e^(5*atanh(a*x))*x^2/(c-a^2*c*x^2)^(27/2)".ToEntity().Integrate("x")` | `integral(...)` | 125 characters | +| `"e^(1/2*atanh(a*x))*x^2/(c-a^2*c*x^2)^(9/8)".ToEntity().Integrate("x")` | `integral(...)` | 121 characters | +| `"e^(2*atanh(a*x))*x^2/(c-a^2*c*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 109 characters | + ### The last digits of an evaluated sine, cosine, tangent, arcsine or arctangent **Silent, and in the hundredth digit.** `EvalNumerical` at the default hundred digits gives the @@ -1387,7 +1404,10 @@ values, and too few points could be compared. The answer is now refused only whe a cosine `g`, so `x/g` is constant on each side of 0; the powers `-1` and `-2` are then sine and cosine integrals of `(a + b u)/(2 b)`, hyperbolic ones for the inverse hyperbolic cosine, and the recurrence run backwards gives the rest. The answers hold on the real domain of the integrand -([#718](https://github.com/asc-community/AngouriMath/issues/718)). +([#718](https://github.com/asc-community/AngouriMath/issues/718)). The six of Rubi's 7.1.5 +`1/(a + b asinh(+-i + d x^2))^n` are answered the same way, with `g = cosh(u/2 + c pi/4)`: `cosh(u/2) + c sinh(u/2)` is +`sqrt(2)` times it, and its square is `1 + c sinh(u) = c d x^2`. They are complex at every `x`, and are checked there +([#233](https://github.com/asc-community/AngouriMath/issues/233)). | Input | Was (2.5.0) | Now | |---|---|---| @@ -1397,6 +1417,8 @@ rest. The answers hold on the real domain of the integrand | `"1/(a+b*asin(-1+d*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 224 characters | | `"1/(a+b*acosh(1+d*x^2))^2".ToEntity().Integrate("x")` | `integral(...)` | 608 characters | | `"1/(a+b*acosh(-1+d*x^2))^3".ToEntity().Integrate("x")` | `integral(...)` | 633 characters | +| `"1/(a+b*asinh(i+d*x^2))".ToEntity().Integrate("x")` | `integral(...)` | 482 characters | +| `"1/(a+b*asinh(-i+d*x^2))^3".ToEntity().Integrate("x")` | `integral(...)` | 735 characters | ### The reciprocal of a binomial in a high power of a sine, cosine or hyperbolic one is split over its roots diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index a1d92a989..ce42563f6 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -1057,6 +1057,124 @@ private static bool IsADerivativeOfAFactorBelowTheBar(Entity expr, Entity.Variab return ((constant * MathS.Pow(determinant, Number.Integer.Create(j + k + EInteger.One))).Simplify() * sum).InnerSimplified; } + /// + /// x^j T^alpha Q^beta, for a power of a quotient of linears T = (p + q x)/(r + s x) + /// and one of a multiple of their product, Q = kappa (p + q x)(r + s x), where the + /// antiderivative is algebraic: T^alpha Q^(beta + 1) P(x) with P of degree + /// j - 1. With D = q r - p s, T'/T = kappa D/Q, so the derivative is + /// T^alpha Q^beta (alpha kappa D P + (beta + 1) Q' P + 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 + 2 beta + 1) p_(m-1) + (alpha kappa D + B (m + beta + 1)) p_m + C (m + 1) p_(m+1) = [m = j] + /// for m from 0 to j: solved from the top, p_(j-1) = 1/(A (j + 2 beta + 1)), + /// down to p_0, and the equation at m = 0 checked. Rubi's 7.3.6 + /// e^(n artanh(a x)) x^2 (c - a^2 c x^2)^(-1 - n^2/2) is + /// e^(n artanh(a x)) (1 - a n x)/(a^3 c n (1 - n^2) (c - a^2 c x^2)^(n^2/2)). + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + /// + /// One equation more than unknowns: where the check fails the antiderivative is not this + /// one term -- it has a logarithm, or is hypergeometric -- and the rule declines. Where it + /// holds, the rules after it ran past the budget: e^(6 artanh(a x)) x^2/(c - a^2 c x^2)^19 + /// was written as a rational function over a denominator of degree 38 and handed to the + /// partial fractions, and is (6 a x - 1)/(210 a^3 c^19 (1 - a x)^21 (1 + a x)^15); + /// e^(5 artanh(a x)) x^2/(c - a^2 c x^2)^(27/2) the same, and + /// e^(artanh(a x)/2) x^2/(c - a^2 c x^2)^(9/8) was declined. Where the denominator + /// is small the partial fractions answer it too, as a sum over its two roots: + /// e^(2 artanh(a x)) x^2/(c - a^2 c x^2)^3 is here (2 a x - 1) T Q^(-2)/(6 a^3 c). + /// The Moebius substitution after it takes a power of Q that fixes the degree of + /// P; this takes the degree from x^j. + /// + /// + /// The powers are differentiated as principal powers wherever T and Q are not 0, + /// on a cut as off it, so the answer holds where the integrand is defined, and is checked + /// at sampled points as well. Only for x above the bar: without it P has no + /// degree to take. + /// + /// + internal static Entity? SolveAPowerOfXBesidePowersOfAQuotientOfLinearsAndTheirProductAlgebraically(Entity expr, Entity.Variable x) + { + Entity? quotient = null, alpha = null, quadratic = null, beta = null; + Entity constant = Number.Integer.One; + var j = 0; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + if (!underneath && APowerOfTheVariable(factor, x) is { } k && k > 0) + { + j += k; + continue; + } + if (quotient is null && AsAPowerOfAQuotient(factor, x) is var (written, exponent)) + { + (quotient, alpha) = (written, underneath ? -exponent : exponent); + continue; + } + var (@base, power) = factor is Powf(var b, var e) ? (b, e) : (factor, (Entity)Number.Integer.One); + if (quadratic is null && !power.ContainsNode(x) && @base is not Divf) + { + (quadratic, beta) = (@base, underneath ? -power : power); + continue; + } + return null; + } + if (quotient is not Divf(var top, var bottom) || alpha is null || quadratic is null || beta is null || j < 1 || j > 12 + || !TreeAnalyzer.TryGetPolyLinear(top, x, out var q, out var p) || TreeAnalyzer.IsZero(q) + || !TreeAnalyzer.TryGetPolyLinear(bottom, x, out var s, out var r) || TreeAnalyzer.IsZero(s) + || !TreeAnalyzer.TryGetPolyQuadratic(quadratic, x, out var a2, out var b2, out var c2) || TreeAnalyzer.IsZero(a2)) + return null; + var determinant = (q * r - p * s).InnerSimplified; + var kappa = Functions.PartialFractions.Bare((a2 / (q * s)).Simplify()); + if (VanishesIdentically(determinant) + || !VanishesIdentically(b2 - kappa * (p * s + q * r)) || !VanishesIdentically(c2 - kappa * p * r)) + return null; + (alpha, beta) = (alpha.InnerSimplified, beta.InnerSimplified); + var alphaKappaD = (alpha * kappa * determinant).InnerSimplified; + // p_0 .. p_(j - 1), and p_j = 0, from the top down. + var coefficients = new Entity[j + 1]; + coefficients[j] = Number.Integer.Zero; + for (var m = j; m >= 1; m--) + { + var divisor = (a2 * (m + 2 * beta + 1).Expand()).InnerSimplified; + if (VanishesIdentically(divisor)) + return null; + Entity rest = m == j ? Number.Integer.One + : -((alphaKappaD + b2 * (m + beta + 1).Expand()) * coefficients[m] + c2 * (m + 1) * coefficients[m + 1]); + coefficients[m - 1] = Functions.PartialFractions.Bare((rest / divisor).Simplify()); + } + var second = j >= 2 ? coefficients[1] : Number.Integer.Zero; + if (!VanishesIdentically((alphaKappaD + b2 * (beta + 1)) * coefficients[0] + c2 * second)) + return null; + Entity polynomial = Number.Integer.Zero; + for (var m = 0; m < j; m++) + if (!VanishesIdentically(coefficients[m])) + polynomial += coefficients[m] * (m == 0 ? Number.Integer.One : m == 1 ? x : MathS.Pow(x, m)); + var answer = (constant * MathS.Pow(quotient, alpha) * MathS.Pow(quadratic, (beta + 1).Expand().InnerSimplified) * polynomial).InnerSimplified; + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; + + // T^alpha as it is written, or as e^(alpha ln(T)) where it is not folded yet. + static (Entity Quotient, Entity Exponent)? AsAPowerOfAQuotient(Entity factor, Entity.Variable x) + { + if (factor is Powf(Divf(var above, var below) quotient, var exponent) && !exponent.ContainsNode(x) + && above.ContainsNode(x) && below.ContainsNode(x)) + return (quotient, exponent); + if (factor is not Powf(var e, var power) || e != MathS.e || !power.ContainsNode(x)) + return null; + var logarithms = power.Nodes.OfType().Where(node => node.Base == MathS.e && node.ContainsNode(x)).Distinct().ToList(); + if (logarithms.Count != 1 || logarithms[0].Antilogarithm is not Divf(var top, var bottom) || !top.ContainsNode(x) || !bottom.ContainsNode(x)) + return null; + var placeholder = Variable.CreateUnique(power, "l_quot"); + if (!TreeAnalyzer.TryGetPolyLinear(power.Replace(node => node == logarithms[0] ? placeholder : node), placeholder, out var multiple, out var free) + || multiple.ContainsNode(x) || !TreeAnalyzer.IsZero(free)) + return null; + return (logarithms[0].Antilogarithm, multiple); + } + } + /// /// A function of the hyperbolic tangent of a linear, under u = tanh(a + b x): the parser /// writes tanh(L) as (W - 1)/(W + 1) and coth(L) as (W + 1)/(W - 1) @@ -37776,15 +37894,18 @@ static bool HoldsAPythagoreanPair(Entity sum) /// negative power goes through u: w - c is a constant times g(u/2 + beta)^2 /// for a sine or a cosine g, so x/g is constant on each side of 0, and the /// powers -1 and -2 are sine and cosine integrals of (a + b u)/(2 b) - /// (hyperbolic ones for the hyperbolic cosine), the rest the recurrence run backwards. Exact on - /// the real domain; off it x/g can change sign across a cut. + /// (hyperbolic ones for the hyperbolic cosine and sine), the rest the recurrence run backwards. + /// Exact on the real domain; off it x/g can change sign across a cut. For the inverse + /// hyperbolic sine g is cosh(u/2 + c pi/4): cosh(u/2) + c sinh(u/2) is + /// sqrt(2) times it, and its square is 1 + c sinh(u) = c d x^2. /// /// /// Each u has u'' = e w u'^3, and with c^2 = 1 (or -1) the radicand of /// u' is d x^2 (2 c + d x^2) up to sign, so H u' = 2 and /// (x (a + b u)^n + 2 e n b H (a + b u)^(n - 1))' is (a + b u)^n plus a multiple of /// (a + b u)^(n - 2). Rubi's 5.1.5, 5.2.5 and 7.2.5 (a + b acos(1 + d x^2))^4 - /// and the like were declined. + /// and the like were declined, and 7.1.5 1/(a + b asinh(i + d x^2))^n, which had no + /// row in the table below. /// https://github.com/asc-community/AngouriMath/issues/718 /// internal static Entity? SolveAPowerOfAnInverseFunctionOfAQuadraticByTheRecurrence(Entity expr, Entity.Variable x) @@ -37920,6 +38041,9 @@ static bool HoldsAPythagoreanPair(Entity sum) (Arcsinf, _, _, false) => (false, -MathS.pi / 4, false, false, -1), (Logf, -1, 1, true) => (true, Number.Integer.Zero, true, true, 1), (Logf, -1, 1, false) => (false, Number.Integer.Zero, true, false, 1), + // asinh with c = +-i: cosh(u/2) + c sinh(u/2) is sqrt(2) cosh(u/2 + c pi/4), whose + // square is 1 + c sinh(u) = c d x^2. + (Logf, -1, -1, _) => (false, c * MathS.pi / 4, true, false, 1), _ => null, }; if (table is not var (gIsSine, beta, hyperbolic, cosine, hSign)) @@ -37950,6 +38074,10 @@ static bool HoldsAPythagoreanPair(Entity sum) integrals[k - 2] = (x * MathS.Pow(a_bu, k) + 2 * e * k * b * h * MathS.Pow(a_bu, k - 1) - integrals[k]) / (4 * e * MathS.Sqr(b) * k * (k - 1)); var answer = integrals[n].InnerSimplified; + // With c = +-i the integrand is real nowhere, so it is compared where it is complex: + // the argument of the logarithm stays off its cut there, its imaginary part the sign of c. + if (squareAtZero == -1) + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; return Functions.PartialFractions.DerivativeDiffersWhereReal(answer, expr, x) ? null : answer; } } diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index a35a41750..309a1ac1a 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -687,6 +687,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A power of a sine, cosine or hyperbolic one over a binomial in a power of it, over the // binomial's roots, before the rules that spent the budget on `sinh(y)^4/(a + b sinh(y)^3)`. if ((answer = IndefiniteIntegralSolver.SolveAPowerOverABinomialInAPowerOfASineOverItsRoots(expr, x)) is { }) return answer; + // x^j T^alpha Q^beta for a quotient of linears T and a multiple Q of their product, where the + // antiderivative is T^alpha Q^(beta + 1) times a polynomial of degree j - 1, before the rewrite + // below hands `e^(6 artanh(a x)) x^2/(c - a^2 c x^2)^19` to partial fractions of degree 38. + if ((answer = IndefiniteIntegralSolver.SolveAPowerOfXBesidePowersOfAQuotientOfLinearsAndTheirProductAlgebraically(expr, x)) is { }) return answer; // An exponential of a whole multiple of an inverse hyperbolic tangent, written algebraically, // before the rules that timed out on `e^(artanh(a x)) x^2/(c - a^2 c x^2)^2`. if ((answer = IndefiniteIntegralSolver.SolveAnExponentialOfAnInverseHyperbolicTangentAlgebraically(expr, x)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs index 30284e6a1..0befaf528 100644 --- a/Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/APowerOfAnInverseFunctionOfAQuadraticIntegralTest.cs @@ -76,6 +76,19 @@ public void ANegativePowerIsAnsweredAndDifferentiatesBack(string integrand, int DifferentiatesBack(integrand, ("a", -0.75), ("b", -0.7), ("d", 0.15 * sign)); } + // A negative power of the inverse hyperbolic sine of c + d x^2 with c = +-i, complex at every + // x, in hyperbolic sine and cosine integrals: cosh(u/2) + c sinh(u/2) is sqrt(2) cosh(u/2 + c pi/4). + [Theory] + [InlineData("1/(a+b*asinh(i+d*x^2))", 1)] + [InlineData("1/(a+b*asinh(i+d*x^2))^2", -1)] + [InlineData("1/(a+b*asinh(-i+d*x^2))^3", 1)] + [InlineData("1/(a+b*asinh(-i+d*x^2))^4", -1)] + public void ANegativePowerOfAnInverseHyperbolicSineIsAnsweredAndDifferentiatesBack(string integrand, int sign) + { + DifferentiatesBack(integrand, ("a", 0.6), ("b", 1.25), ("d", 0.3 * sign)); + DifferentiatesBack(integrand, ("a", -0.75), ("b", -0.7), ("d", 0.15 * sign)); + } + // Half an odd power of the inverse hyperbolic cosine, in the error functions of // sqrt(a + b u)/sqrt(2 b). [Theory] diff --git a/Sources/Tests/UnitTests/Calculus/APowerOfXBesidePowersOfAQuotientOfLinearsAndTheirProductIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfXBesidePowersOfAQuotientOfLinearsAndTheirProductIntegralTest.cs new file mode 100644 index 000000000..ec6dfaff9 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerOfXBesidePowersOfAQuotientOfLinearsAndTheirProductIntegralTest.cs @@ -0,0 +1,56 @@ +// +// 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 T^alpha Q^beta, a power of a quotient of linears beside a power of a multiple of + /// their product, where the antiderivative is T^alpha Q^(beta + 1) times a polynomial of + /// degree j - 1: e^(n artanh(a x)) x^2 (c - a^2 c x^2)^(-1 - n^2/2) was left + /// unevaluated, and the rest ran past the budget. Rubi's 7.3.6. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APowerOfXBesidePowersOfAQuotientOfLinearsAndTheirProductIntegralTest + { + [Theory] + [InlineData("e^(n*atanh(a*x))*x^2*(c - a^2*c*x^2)^(-1 - n^2/2)")] + [InlineData("e^(atanh(a*x)/2)*x^2/(c - a^2*c*x^2)^(9/8)")] + [InlineData("e^(6*atanh(a*x))*x^2/(c - a^2*c*x^2)^19")] + [InlineData("e^(5*atanh(a*x))*x^2/(c - a^2*c*x^2)^(27/2)")] + [InlineData("x^2/(e^(5*atanh(a*x))*(c - a^2*c*x^2)^(27/2))")] + [InlineData("x*((x + 2)/(3*x - 1))^(1/3)*(3*x^2 + 5*x - 2)^(-8/15)")] + public void AsTheQuotientAndTheQuadraticTimesAPolynomial(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 2000, $"{text.Length} characters of answer for {integrand}"); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("c", 0.7).Substitute("n", 0.9); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + // On both sides of the roots, -0.77 and 0.77 for a = 1.3, and of 1/3 and -2. + foreach (var at in new[] { -2.9, -1.4, -0.4, 0.2, 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}"); + } + } +}