diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index a224d3b20..7d193e04f 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -6593,6 +6593,25 @@ leaves are read as quotients of the sine and the cosine | `"1/(a+b*asin(c*x))".Integrate("x")` | `integral(1 / (a + b * arcsin(c * x)), x)` | `Si` and `Ci` of `(a + b arcsin(c x))/b` | | `"x/((c+a^2*c*x^2)^2*atan(a*x))".Integrate("x")` | `integral(x / ((c + a ^ 2 * c * x ^ 2) ^ 2 * arctan(a * x)), x)` | `1/2 * Si(2 * arctan(a * x)) / (a ^ 2 * c ^ 2) + C` | +### A power of x over a power of a quadratic times an inverse trigonometric function is integrated by parts against one term + +**Answers where there were none.** `x^3 (a + b F(c x))/(d + e x^2)^3` for `F` the inverse sine, tangent, secant and +cosecant, Rubi's 5.1.4, 5.3.4, 5.5.1 and 5.6.1, was declined. By parts against the integrator's own antiderivative of +`x^3/(d + e x^2)^3`, a piecewise on `e = 0`, spent the budget on every arm. Where `m + 2p + 3 = 0`, +`x^m (d + e x^2)^p` has the one-term antiderivative `v = x^(m + 1) (d + e x^2)^(p + 1)/(d (m + 1))`, Rubi's, which +vanishes at 0, and the integral is now `v (a + b F(c x)) - b integral v F(c x)'` with `F(c x)'` algebraic. +`v' = x^m (d + e x^2)^p` is an identity of rational functions, so the parts hold wherever the integrand is defined, +on any branch of `F`. Only +for a whole `p` of `-3` or less, and not where `d + e x^2` is a multiple of the radicand `F` holds, which the +substitution under `F` answers ([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x^3*(a + b*arcsin(c*x))/(d + k*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 6,013 characters | +| `"x^3*(a + b*arctan(c*x))/(d + k*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 4,066 characters | +| `"x^3*(a + b*arcsec(c*x))/(d + k*x^2)^3".ToEntity().Integrate("x")` | `integral(...)` | 1,772 characters | +| `"x^5*(a + b*arctan(c*x))/(d + k*x^2)^4".ToEntity().Integrate("x")` | `integral(...)` | 6,176 characters | + ### The arcsecant beside a root of `x^2 - 1` is integrated right below `-1` too **Answers where there were none.** Under `x = sec(u)`, the substitution that undoes the arcsecant, diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 1c01f2891..1aae3fb8c 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -32123,6 +32123,134 @@ Piecewise piecewise when piecewise.Cases.Any() => piecewise.Cases.First().Expres return (alpha / (p0 * slope) * product).InnerSimplified; } + /// + /// x^m (d + e x^2)^p (a + b F(c x)), F an inverse trigonometric function and + /// m + 2p + 3 = 0, by parts against the antiderivative of x^m (d + e x^2)^p that is + /// one term, v = x^(m + 1) (d + e x^2)^(p + 1)/(d (m + 1)): + /// integral = v (a + b F(c x)) - b integral v F(c x)', and F(c x)' is algebraic. + /// + /// + /// + /// v' = x^m (d + e x^2)^p ((m + 1)(d + e x^2) + 2 (p + 1) e x^2)/(d (m + 1)), and the bracket + /// is (m + 1) d + (m + 2p + 3) e x^2, which is (m + 1) d for every e: an + /// identity of rational functions, so the parts are exact wherever the integrand is defined, on + /// any branch of F. It is Rubi's choice of v, and the one that vanishes at + /// x = 0. + /// + /// + /// The integrator's own antiderivative of x^3/(d + e x^2)^3 is a piecewise on e = 0, + /// and by parts against it spent the budget on every arm: Rubi's 5.1.4, 5.3.4, 5.5.1 and 5.6.1 + /// x^3 (a + b F(c x))/(d + e x^2)^3, for F the inverse sine, tangent, secant and + /// cosecant, ran out of time, where the remainder against this v is answered in under a + /// second. + /// + /// + /// Narrow on purpose. p whole and at most -3, so that v is rational and + /// m = -2p - 3 at least 3: m = 1, x/(d + e x^2)^2, is answered by parts as it + /// was, and a half-odd p is answered as it was beside the inverse sine and tangent and is + /// elliptic beside the secant. And not for e = +-c^2 d, the quadratic a multiple of + /// F's own radicand, which the substitution under F answers. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveAPowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByParts(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked || integratingAPowerOverAQuadraticBesideAnInverse) + return null; + // The one inverse trigonometric function, of c x. + Entity? function = null; + foreach (var node in expr.Nodes) + if ((node is Arcsinf or Arccosf or Arctanf or Arccotanf or Arcsecantf or Arccosecantf) && node.ContainsNode(x)) + { + if (function is not null && function != node) + return null; + function = node; + } + var argument = function switch + { + Arcsinf(var u) => u, Arccosf(var u) => u, Arctanf(var u) => u, + Arccotanf(var u) => u, Arcsecantf(var u) => u, Arccosecantf(var u) => u, + _ => null + }; + if (function is null || argument is null || !TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out var offset) + || rate.ContainsNode(x) || !Functions.PartialFractions.IsZeroAsAValue(offset) || Functions.PartialFractions.IsZeroAsAValue(rate)) + return null; + // Above the bar one factor linear in it, a + b F, and the rest; below, nothing with it. + var (numerator, denominator) = Functions.SingleQuotient.Of(expr); + if (denominator.ContainsNode(function)) + return null; + Entity? linear = null; + Entity rest = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(numerator)) + { + if (!factor.ContainsNode(function)) + rest *= factor; + else if (linear is null) + linear = factor; + else + return null; + } + var t = Variable.CreateUnique(expr, "t_inv"); + if (linear is null || !TreeAnalyzer.TryGetPolyLinear(linear.Substitute(function, t), t, out var slope, out var intercept) + || slope.ContainsNode(x) || intercept.ContainsNode(x) || Functions.PartialFractions.IsZeroAsAValue(slope)) + return null; + // The rest over the denominator as K x^m Q^p, Q = d + e x^2. + Entity constant = Number.Integer.One; + Entity? quadratic = null; + var m = 0; + var p = 0; + foreach (var (side, sign) in new[] { (rest, 1), (denominator, -1) }) + foreach (var factor in Mulf.LinearChildren(side)) + { + if (!factor.ContainsNode(x)) + { + constant = sign > 0 ? constant * factor : constant / factor; + continue; + } + var (@base, power) = factor is Powf(var b0, Number.Integer { EInteger: var e0 }) && e0.CanFitInInt32() + ? (b0, e0.ToInt32Unchecked()) : (factor, 1); + if (@base == x) + m += sign * power; + else if (quadratic is null || quadratic == @base) + { + quadratic = @base; + p += sign * power; + } + else + return null; + } + if (quadratic is null || p > -3 || m != -2 * p - 3 + || !TreeAnalyzer.TryGetPolynomial(quadratic, x, out var read) + || read.Keys.Any(key => !key.IsZero && !key.Equals(EInteger.FromInt32(2))) + || !read.TryGetValue(EInteger.Zero, out var d) || !read.TryGetValue(EInteger.FromInt32(2), out var e) + || d.ContainsNode(x) || e.ContainsNode(x) + || Functions.PartialFractions.IsZeroAsAValue(d) || Functions.PartialFractions.IsZeroAsAValue(e) + || Functions.PartialFractions.IsZeroAsAValue((e - rate * rate * d).Expand()) + || Functions.PartialFractions.IsZeroAsAValue((e + rate * rate * d).Expand())) + return null; + var v = constant * MathS.Pow(x, Number.Integer.Create(m + 1)) * MathS.Pow(quadratic, Number.Integer.Create(p + 1)) + / (d * (m + 1)); + // The remainder holds no inverse function, so it is not this integrand again; the flag is for + // whatever its own search writes. + var remainder = (slope * v * function.Differentiate(x)).InnerSimplified; + Entity? integral; + integratingAPowerOverAQuadraticBesideAnInverse = true; + try + { + integral = Integration.ComputeAsAQuestionOfItsOwn(remainder, x, integrateByParts); + } + finally + { + integratingAPowerOverAQuadraticBesideAnInverse = false; + } + if (integral is null) + return null; + var answer = v * linear - integral; + return Functions.PartialFractions.DerivativeDiffersWhereReal(answer, expr, x) ? null : answer; + } + + [System.ThreadStatic] private static bool integratingAPowerOverAQuadraticBesideAnInverse; + /// /// x^m (a + b atan(c x))^p over a power of d (1 +- i c x), under x = tan(u)/c, where /// 1 + i tan(u) is e^(i u)/cos(u): the integrand in u is diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 5d138bd72..21a4510fe 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -707,6 +707,10 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // A rational function times the logarithm of a nested product of powers of linears, by parts, // its logarithm's derivative read off as written: `ln(e (f (a + b x)^p (c + d x)^q)^r)`. if ((answer = IndefiniteIntegralSolver.SolveARationalFunctionTimesALogarithmOfPowersOfLinears(expr, x, integrateByParts)) is { }) return answer; + // x^m (d + e x^2)^p (a + b F(c x)) with m + 2p + 3 = 0 and F an inverse trigonometric function, by + // parts against v = x^(m + 1) (d + e x^2)^(p + 1)/(d (m + 1)), before by parts against the + // integrator's piecewise v timed out on `x^3 (a + b asin(c x))/(d + e x^2)^3`. + if ((answer = IndefiniteIntegralSolver.SolveAPowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByParts(expr, x, integrateByParts)) is { }) return answer; // A function of an inverse of something not linear times its derivative up to a factor // constant though written with x: `atan(c x/sqrt(a - c^2 x^2))^2/sqrt(d - c^2 d x^2/a)`. Early, // since it fires only where that factor holds x, and the rules before it spent 28 seconds on diff --git a/Sources/Tests/UnitTests/Calculus/APowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByPartsIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/APowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByPartsIntegralTest.cs new file mode 100644 index 000000000..302a84283 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/APowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByPartsIntegralTest.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 System.Globalization; +using System.Linq; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// x^m (d + k x^2)^p (a + b F(c x)) with m + 2p + 3 = 0 and F an inverse + /// trigonometric function, by parts against v = x^(m + 1) (d + k x^2)^(p + 1)/(d (m + 1)). + /// By parts against the integrator's own antiderivative, a piecewise on k = 0, ran out of + /// time. Rubi's 5.1.4, 5.3.4, 5.5.1 and 5.6.1. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class APowerOfXOverAPowerOfAQuadraticTimesAnInverseTrigonometricFunctionByPartsIntegralTest + { + private const string Inside = "-1.3 -0.9 -0.4 0.3 0.8 1.2"; + private const string Outside = "-4.1 -2.9 -2.1 2.1 2.9 4.1"; + + [Theory] + [InlineData("x^3*(a + b*arcsin(c*x))/(d + k*x^2)^3", Inside)] + [InlineData("x^3*(a + b*arctan(c*x))/(d + k*x^2)^3", Inside)] + [InlineData("x^3*(a + b*arcsec(c*x))/(d + k*x^2)^3", Outside)] + [InlineData("x^3*(a + b*arccsc(c*x))/(d + k*x^2)^3", Outside)] + [InlineData("x^3*(a + b*arccos(c*x))/(d + k*x^2)^3", Inside)] + [InlineData("x^3*(a + b*arccot(c*x))/(d + k*x^2)^3", Inside)] + [InlineData("x^3*arcsin(c*x)/(d + k*x^2)^3", Inside)] + [InlineData("x^5*(a + b*arctan(c*x))/(d + k*x^2)^4", Inside)] + public void ByPartsAgainstTheOneTermAntiderivative(string integrand, string points) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.True(text.Length < 20000, $"{text.Length} characters of answer for {integrand}"); + foreach (var (d, k) in new[] { (2.1, 0.9), (-1.7, 1.3), (1.9, -0.8) }) + { + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.6) + .Substitute("d", d).Substitute("k", k); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + foreach (var at in points.Split(' ').Select(point => double.Parse(point, CultureInfo.InvariantCulture))) + { + 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}, d = {d}, k = {k}, where the integrand is {want}"); + } + Assert.True(compared >= 5, $"only {compared} points could be compared for {integrand} at d = {d}, k = {k}"); + } + } + } +}