diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index db6f97ad8..89c0b89fc 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -6394,6 +6394,29 @@ no integrand without one changes: Rubi's independent suites and a sample of its | `"e^(b*x)*Ei(b*x)/x".Integrate("x")` | `Ei * e ^ (b * x) + C`, with `Ei` a variable | `Ei(b * x) ^ 2 / 2 + C` | | `"Si(b*x)*sin(b*x)/x".Integrate("x")` | `Si * -cos(b * x) + C`, with `Si` a variable | `Si(b * x) ^ 2 / 2 + C` | +### A sum of powers of an inverse trigonometric or error function is integrated by parts down the powers + +**Answers where there were none.** `x^3/((1 + a^2 x^2) atan(a x)^3) - 3 x^2/(2 a atan(a x)^2)`, +`x asin(x)^(3/2)/(1 - x^2)^2 - 3 x/(8 (1 - x^2) sqrt(asin(x)))` and `erf(b x) e^(-b^2 x^2) (1/x^3 + b^2/x)` with its +`erfc` and `erfi` twins, Rubi's 5.1.4, 5.3.4 and 8.1, were declined: no term is elementary alone, and the split and +the parts each read one term at a time. A sum of terms each a power of one function `K`, whose derivative is free of +it, times a cofactor free of it is now integrated by parts from the top power, +`integral A K^p = F K^p - integral p F K^(p - 1) K'` with `F` the antiderivative of `A`, and the last integral joins +the terms a power below. Where every power cancels so, the +parts are the answer, `(-x^3/(2 a atan(a x)^2))'` for the first; where what is left is free of `K` it is asked. The +powers a whole number apart and at most five apart, a cancelling group decided at sampled points, and the answer +differentiated back. Three sums of the same kind the unreleased master answered come back shorter +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x^3/((1+a^2*x^2)*arctan(a*x)^3)-3/2*x^2/(a*arctan(a*x)^2)".ToEntity().Integrate("x")` | `integral(...)` | `-x^3/(2 a arctan(a x)^2)`, 51 characters | +| `"x*arcsin(x)^(3/2)/(1-x^2)^2-3/8*x/((1-x^2)*sqrt(arcsin(x)))".ToEntity().Integrate("x")` | `integral(...)` | 99 characters | +| `"erf(b*x)/(e^(b^2*x^2)*x^3)+b^2*erf(b*x)/(e^(b^2*x^2)*x)".ToEntity().Integrate("x")` | `UnrecognizedFunctionParseException`: there is no function `erf` | 219 characters | +| `"erfi(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfi(b*x)/(e^(b^2*x^2)*x)".ToEntity().Integrate("x")` | `integral(...)`, reading `erfi` as a symbol | 131 characters | +| `"erf(x)^2+4*x*erf(x)*e^(-x^2)/sqrt(pi)".ToEntity().Integrate("x")` | `UnrecognizedFunctionParseException`; 352 characters on the unreleased master | `x erf(x)^2`, 18 characters | +| `"2*x/arccos(x)+x^2/(sqrt(1-x^2)*arccos(x)^2)".ToEntity().Integrate("x")` | `integral(...)`; 121 characters on the unreleased master | `x^2/arccos(x)`, 33 characters | + ### A square of a special function, and one beside a power of `x` and its elementary derivative, are integrated by parts `x Si(b x) sin(b x)` was left unevaluated where `Si(b x) sin(b x)` was not: by parts it is diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 3fd08cfbb..94b9d998a 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -15991,6 +15991,137 @@ static bool More(Entity difference) return changed ? rebuilt / bottom : null; } + /// + /// A sum of terms each a power of one function K whose derivative is free of it -- an + /// inverse trigonometric function or an error function -- times a cofactor free of it, + /// integrated by parts down the powers: ∫ A K^p = F K^p - ∫ p F K^(p - 1) K' for + /// F = ∫ A, the cofactor of the top power, and the last integral joins the terms a + /// power below. Where every power cancels so, the parts are the antiderivative: + /// x^3/((1 + a^2 x^2) atan(a x)^3) - 3 x^2/(2 a atan(a x)^2) is + /// (-x^3/(2 a atan(a x)^2))', in one step, and + /// x asin(x)^(3/2)/(1 - x^2)^2 - 3 x/(8 (1 - x^2) sqrt(asin(x))) in two. Where what is + /// left is free of K, it is asked: erf(b x) e^(-b^2 x^2) (1/x^3 + b^2/x) leaves + /// b e^(-2 b^2 x^2)/(sqrt(pi) x^2). + /// + /// + /// + /// The terms are not elementary apart -- ∫ x/((1 - x^2) sqrt(asin(x))) is + /// ∫ tan(t)/sqrt(t) dt -- so the split declined each, and by parts reads one term + /// at a time. Rubi's 5.1.4, 5.3.4 and 8.1 write such sums to test that the parts meet. + /// + /// + /// One such function in the sum, every term a rational power of it times a cofactor free + /// of it, the powers a whole number apart, at most five apart, and each step's + /// F free of K. A group that cancels is decided at sampled points, and the + /// answer is checked by differentiating it back. A remainder free of K is asked + /// only for a top power of two terms or more with a cofactor among them that is not + /// algebraic, which is where the cofactors integrate together and not apart: a sum of + /// algebraic cofactors, a + b asin(c x), is what the split answers term by term, + /// and this rule would only respell it. At the question asked only, and not inside its + /// own questions, which are free of K anyway. + /// + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + internal static Entity? SolveASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowers(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked || expr is not (Sumf or Minusf) || descendingThePowers) + return null; + Entity? kernel = null; + foreach (var node in expr.Nodes) + if (node is Arcsinf or Arccosf or Arctanf or Arccotanf or Erff or Erfcf or Erfif && node.ContainsNode(x)) + { + if (kernel is null) + kernel = node; + else if (node != kernel) + return null; + } + if (kernel is null) + return null; + var slope = kernel.Differentiate(x).InnerSimplified; + if (slope.ContainsNode(kernel)) + return null; + Number.Rational? PowerOfTheKernel(Entity factor) + => factor == kernel ? Number.Integer.One + : factor is Powf(var @base, Number.Rational power) && PowerOfTheKernel(@base) is { } inner ? inner * power + : null; + var terms = new List<(Number.Rational Power, Entity Cofactor)>(); + foreach (var term in Sumf.LinearChildren(expr)) + { + Number.Rational power = Number.Integer.Zero; + Entity cofactor = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(term)) + if (!factor.ContainsNode(kernel)) + cofactor *= factor; + else if (PowerOfTheKernel(factor) is { } read) + power += read; + else + return null; + terms.Add((power, cofactor)); + } + var top = terms.Select(term => term.Power).Aggregate((l, r) => Above(r, l) ? r : l); + var lowest = terms.Select(term => term.Power).Aggregate((l, r) => Above(l, r) ? r : l); + if (terms.Any(term => !(top - term.Power).ERational.Denominator.Equals(EInteger.One)) + || Above(top - lowest, Number.Integer.Create(5))) + return null; + var groups = new List<(Number.Rational Power, Entity Cofactor)>(); + foreach (var (power, cofactor) in terms) + Join(groups, power, cofactor); + var atTheTop = terms.Where(term => Same(term.Power, top)).ToList(); + var mayAskWhatIsLeft = atTheTop.Count >= 2 && atTheTop.Any(term => !IsAlgebraicIn(term.Cofactor, x)); + if (groups.Count < 2 && !mayAskWhatIsLeft) + return null; + descendingThePowers = true; + try + { + Entity parts = Number.Integer.Zero; + for (var step = 0; step < 8; step++) + { + groups.RemoveAll(group => Functions.PartialFractions.HoldsAtSampledPoints(group.Cofactor, Number.Integer.Zero, x)); + if (groups.Count == 0) + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(parts, expr, x) ? parts : null; + var (power, cofactor) = groups.Aggregate((l, r) => Above(r.Power, l.Power) ? r : l); + if (power.ERational.IsZero) + { + // Free of K: asked where the top power's cofactors integrate together only. + if (groups.Count > 1 || !mayAskWhatIsLeft + || Integration.ComputeAsAQuestionOfItsOwn(cofactor, x, integrateByParts) is not { } rest) + return null; + var answer = parts + rest; + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; + } + // Below zero the powers never reach a remainder free of K, so below the lowest + // written one nothing is left to cancel. + if (power.ERational.Sign < 0 && Above(lowest, power) + || Integration.ComputeAsAQuestionOfItsOwn(cofactor, x, integrateByParts) is not { } integrated + || integrated.ContainsNode(kernel)) + return null; + parts += integrated * (Same(power, Number.Integer.One) ? kernel : MathS.Pow(kernel, power)); + groups.RemoveAll(group => Same(group.Power, power)); + Join(groups, power - Number.Integer.One, -power * integrated * slope); + for (var i = 0; i < groups.Count; i++) + groups[i] = (groups[i].Power, groups[i].Cofactor.InnerSimplified); + } + return null; + } + finally + { + descendingThePowers = false; + } + static bool Same(Number.Rational a, Number.Rational b) => a.ERational.CompareTo(b.ERational) == 0; + static bool Above(Number.Rational a, Number.Rational b) => a.ERational.CompareTo(b.ERational) > 0; + static void Join(List<(Number.Rational Power, Entity Cofactor)> groups, Number.Rational power, Entity cofactor) + { + var at = groups.FindIndex(group => Same(group.Power, power)); + if (at < 0) + groups.Add((power, cofactor)); + else + groups[at] = (power, groups[at].Cofactor + cofactor); + } + } + + /// Set while asks its questions. + [System.ThreadStatic] private static bool descendingThePowers; + /// /// A product of powers times a sum that is the derivative of the product with some of /// the powers raised by one: e^x x^2 ln(x)^2 (3 + (3 + x) ln(x)) is diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 9223eb21e..eb38403d9 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -932,6 +932,11 @@ private static Entity WithItsPiecewiseSumsCombined(Entity expr) // recurrence that raises the whole power to -1: closed, where the substitution for the // root makes a rational function over a power of a symbolic quadratic. if ((answer = IndefiniteIntegralSolver.SolveANegativePowerOfALinearBesideAPowerOfAnotherByTheRecurrence(expr, x, integrateByParts)) is { }) return answer; + // A sum of powers of one inverse trigonometric or error function beside cofactors free + // of it, by parts down the powers from the top: `x^3/((1 + a^2 x^2) atan(a x)^3) - + // 3 x^2/(2 a atan(a x)^2)` is the derivative of one product, and neither term is + // elementary on its own. Before the split, which loses it. + if ((answer = IndefiniteIntegralSolver.SolveASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowers(expr, x, integrateByParts)) is { }) return answer; // A product of powers times a sum that is the derivative of the product with some // of the powers raised by one -- `e^x x^2 ln(x)^2 (3 + (3 + x) ln(x))` -- read off // the sum, for symbolic exponents too. Before the split, which loses it. diff --git a/Sources/Tests/UnitTests/Calculus/ASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowersIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowersIntegralTest.cs new file mode 100644 index 000000000..0cba42696 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/ASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowersIntegralTest.cs @@ -0,0 +1,60 @@ +// +// 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 +{ + /// + /// A sum of powers of one inverse trigonometric or error function beside cofactors free of + /// it, by parts down the powers from the top. Its terms are not elementary apart, so the + /// split declined each, and the sum was declined or ran out of time. Rubi's 5.1.4, 5.3.4 + /// and 8.1. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class ASumOfPowersOfAnInverseOrErrorFunctionByPartsDownThePowersIntegralTest + { + [Theory] + // The parts cancel every power below the top: one step, and two. + [InlineData("x^3/((1+a^2*x^2)*arctan(a*x)^3)-3/2*x^2/(a*arctan(a*x)^2)")] + [InlineData("x*arcsin(x)^(3/2)/(1-x^2)^2-3/8*x/((1-x^2)*sqrt(arcsin(x)))")] + [InlineData("1/((1-x^2)*arcsin(x)^2)-x/((1-x^2)^(3/2)*arcsin(x))")] + [InlineData("2*x/arccos(x)+x^2/(sqrt(1-x^2)*arccos(x)^2)")] + [InlineData("erf(x)^2+4*x*erf(x)*e^(-x^2)/sqrt(pi)")] + // What is left is free of the function, and the cofactors integrate together only. + [InlineData("erf(b*x)/(e^(b^2*x^2)*x^3)+b^2*erf(b*x)/(e^(b^2*x^2)*x)")] + [InlineData("erfc(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfc(b*x)/(e^(b^2*x^2)*x)")] + [InlineData("erfi(b*x)/(e^(b^2*x^2)*x^3)+b^2*erfi(b*x)/(e^(b^2*x^2)*x)")] + public void ByPartsDownThePowers(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("b", 0.7); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + var compared = 0; + // Inside (0, 1), where the arcsine and the arccosine are real and positive. + foreach (var at in new[] { 0.15, 0.3, 0.45, 0.6, 0.75, 0.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}"); + } + } +}