diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index cee0dcaa1..698516842 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -707,6 +707,21 @@ interval between the poles of the tangent, as the substitution's are | `"(A + B*tan(x))/(a + i*a*tan(x))^2".ToEntity().Integrate("x")` | `integral(...)` | the same | | `"1/(a + i*a*tan(c + d*x))^3".ToEntity().Integrate("x")` | `integral(...)` | the same in `tan(c + d x)` | +### A cosine and a sine over a power of another such sum are integrated through its derivative + +**Answers where there were none.** `sin(x)/(a + b cos(x) + c sin(x))` was declined, with the rest of +Rubi's 4.7.7 quotients of `A + B cos(x) + C sin(x)` by a whole power of `a + b cos(x) + c sin(x)`, +the squares and cubes after searches past the budget. The numerator is a multiple of the +denominator, one of its derivative and a constant, and each power of the denominator is brought +down to its reciprocal by an identity of derivatives; the reciprocal is answered as before +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sin(x)/(a + b*cos(x) + c*sin(x))".ToEntity().Integrate("x")` | `integral(...)` | `c x/(b^2 + c^2) - b ln(a + b cos(x) + c sin(x))/(b^2 + c^2)` and a multiple of the reciprocal's antiderivative | +| `"1/(a + b*cos(x) + c*sin(x))^2".ToEntity().Integrate("x")` | `integral(...)` | `(b sin(x) - c cos(x))/(a + b cos(x) + c sin(x))` less `a` times the reciprocal's antiderivative, over `b^2 + c^2 - a^2` | +| `"(k + p*cos(x) + q*sin(x))/(a + b*cos(x) + c*sin(x))^3".ToEntity().Integrate("x")` | `integral(...)` | powers of the denominator and the reciprocal's antiderivative | + ### An integrand with the imaginary unit in it is not simplified by the rule that scales the variable **Answers where there were none.** The rule that scales the variable by the integrand's one symbol diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index db5f3fcf6..f5c4c2674 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -8885,6 +8885,161 @@ private static (Entity Coefficient, Entity Degree)? TheMonomial(Entity expr, Ent : null; } + /// + /// (A + B cos(y) + C sin(y))/(a + b cos(y) + c sin(y))^n, for a whole n, by + /// the denominator, its derivative and a constant: the numerator is + /// alpha S + beta S' + gamma for S the base, and each power of S + /// comes down to 1/S, which the half-angle tangent answers. + /// + /// + /// + /// With S' = c cos(y) - b sin(y), matching the cosine, the sine and the constant + /// gives alpha = (B b + C c)/(b^2 + c^2), beta = (B c - C b)/(b^2 + c^2) and + /// gamma = A - a alpha; beta S'/S^n integrates to a logarithm or a power of + /// S. For the powers, with T = b sin(y) - c cos(y) = -S', + /// T^2 = b^2 + c^2 - (S - a)^2, from which + /// + /// + /// d/dy (T S^m) = (1 + m) S^(m+1) - a (1 + 2m) S^m + m (a^2 - b^2 - c^2) S^(m-1) + /// + /// + /// which for m = 1 - k writes int S^(-k) through int S^(1-k) and + /// int S^(2-k), down to int 1/S and y. Rubi's + /// (A + B cos(x) + C sin(x))/(a + b cos(x) + c sin(x)) and its kin were declined: + /// under the half-angle tangent each is a rational function over a quadratic with every + /// coefficient a symbol, beside 1 + t^2, and the square of the denominator was + /// answered with coefficients past the fiftieth degree in them. + /// + /// + /// Exact, as the answer for 1/S it is built on: a linear combination and an + /// identity of derivatives. a^2 = b^2 + c^2, where the base is a square of the + /// half angle and the recurrence divides by zero, is declined here, as is + /// b^2 + c^2 = 0, where the base is a + b e^(±i y); both decided. A constant + /// numerator over the first power is 1/S itself, left to the rules that answer it. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + /// + internal static Entity? SolveACosineAndASineOverAPowerOfAnother(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!TryReadAsQuotient(expr, out var numerator, out var denominator)) + return null; + var (@base, power) = denominator is Powf(var inner, var exponent) && exponent.Evaled is Number.Integer { EInteger: var whole } + ? (inner, whole) + : (denominator, EInteger.One); + if (power.Sign <= 0 || power.CompareTo(EInteger.FromInt32(MaximumPowerOfACosineAndASine)) > 0) + return null; + Entity? argument = null; + bool TryRead(Entity sum, out Entity constant, out Entity cosine, out Entity sine) + { + constant = cosine = sine = Number.Integer.Zero; + foreach (var term in Sumf.LinearChildren(sum)) + { + if (!term.ContainsNode(x)) + { + constant = constant == Number.Integer.Zero ? term : constant + term; + continue; + } + Entity coefficient = Number.Integer.One; + Entity? function = null; + foreach (var factor in Mulf.LinearChildren(term)) + if (!factor.ContainsNode(x)) + coefficient = coefficient == Number.Integer.One ? factor : coefficient * factor; + else if (function is null && factor is Sinf or Cosf) + function = factor; + else + return false; + var inner = function!.DirectChildren.First(); + if (argument is not null && inner != argument) + return false; + argument = inner; + if (function is Cosf) + { + if (cosine != Number.Integer.Zero) return false; + cosine = coefficient; + } + else + { + if (sine != Number.Integer.Zero) return false; + sine = coefficient; + } + } + return true; + } + // The base with a constant, a cosine and a sine: with one of the two functions alone the + // half angle reads it, and without the constant it is one cosine turned by a phase, + // which the rotation answers. + if (!TryRead(@base, out var a, out var b, out var c) || a == Number.Integer.Zero || b == Number.Integer.Zero + || c == Number.Integer.Zero || !TryRead(numerator, out var bigA, out var bigB, out var bigC) || argument is null) + return null; + if (power.Equals(EInteger.One) && bigB == Number.Integer.Zero && bigC == Number.Integer.Zero) + return null; + if (!TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out _) || rate.ContainsNode(x) + || rate.Evaled is Number.Complex { IsZero: true }) + return null; + static bool IsZero(Entity value) + { + var simplified = value.InnerSimplified; + return simplified.Evaled is Number.Complex { IsZero: true } + || simplified.Vars.Any() && Functions.PartialFractions.Bare(simplified.Simplify()).Evaled is Number.Complex { IsZero: true }; + } + var squares = MathS.Sqr(b) + MathS.Sqr(c); + var discriminant = MathS.Sqr(a) - squares; + if (IsZero(squares) || IsZero(discriminant)) + return null; + + var s = @base; + var t = b * MathS.Sin(argument) - c * MathS.Cos(argument); + var alpha = (bigB * b + bigC * c) / squares; + var beta = (bigB * c - bigC * b) / squares; + var gamma = bigA - a * alpha; + // int S^(-k) dx, from int 1/S as the integrator answers it and x. + Entity? reciprocal = null; + var computed = new Dictionary(); + Entity? Reciprocal(int k) + { + if (k == 0) + return x; + if (computed.TryGetValue(k, out var known)) + return known; + Entity? result; + if (k == 1) + result = reciprocal ??= Integration.ComputeIndefiniteIntegral(1 / s, x, integrateByParts); + else + { + if (Reciprocal(k - 1) is not { } previous || Reciprocal(k - 2) is not { } beforeThat) + return null; + result = (t * MathS.Pow(s, 1 - k) / rate - (2 - k) * beforeThat + a * (3 - 2 * k) * previous) + / ((1 - k) * discriminant); + } + if (result is not null) + computed[k] = result; + return result; + } + var n = power.ToInt32Checked(); + Entity answer = Number.Integer.Zero; + if (alpha != Number.Integer.Zero && !IsZero(alpha)) + { + if (Reciprocal(n - 1) is not { } lower) + return null; + answer += alpha * lower; + } + if (beta != Number.Integer.Zero && !IsZero(beta)) + answer += beta / rate * (n == 1 ? MathS.Ln(s) : MathS.Pow(s, 1 - n) / (1 - n)); + if (gamma != Number.Integer.Zero && !IsZero(gamma)) + { + if (Reciprocal(n) is not { } same) + return null; + answer += gamma * same; + } + return answer == Number.Integer.Zero ? null : answer; + } + + /// + /// The largest power of the denominator + /// steps through: each step is a term of the answer. + /// + private const int MaximumPowerOfACosineAndASine = 6; + /// /// A sum of a cosine and a sine of one argument turned into one cosine: /// a cos(y) + b sin(y) is R cos(u) for R = sqrt(a^2 + b^2) and diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 494173d9e..d956c6a3d 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -877,6 +877,10 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // Several trigonometric arguments that are multiples of one linear with an offset or a // symbolic slope, written in that linear: before the substitution search, which reads // each function on its own. + // A cosine and a sine over a power of another such sum, through the denominator, its + // derivative and a constant, down to the reciprocal of the base: before the substitution + // search, which reads the quotient term by term. + if ((answer = IndefiniteIntegralSolver.SolveACosineAndASineOverAPowerOfAnother(expr, x, integrateByParts)) is { }) return answer; if ((answer = IndefiniteIntegralSolver.SolveByWritingMultiplesOfOneLinearArgument(expr, x, integrateByParts)) is { }) return answer; // A constant out of a fractional power of a trigonometric factor: `sqrt(b sec(x))` is // `sqrt(b) sqrt(sec(x))` for a positive `b`, which meets the other powers of the diff --git a/Sources/Tests/UnitTests/Calculus/CosineAndSineOverAPowerOfAnotherIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/CosineAndSineOverAPowerOfAnotherIntegralTest.cs new file mode 100644 index 000000000..adafc2b9e --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/CosineAndSineOverAPowerOfAnotherIntegralTest.cs @@ -0,0 +1,51 @@ +// +// 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 + B cos(x) + C sin(x))/(a + b cos(x) + c sin(x))^n through the denominator, its + /// derivative and a constant, each power of the denominator brought down to its reciprocal. + /// Rubi's 4.7.7. + /// #718 + /// + /// + /// Checked by differentiating back with the symbols pinned, inside (-pi, pi), where the + /// half-angle tangent the reciprocal is answered in is continuous. + /// + [Trait("Area", "Calculus")] + public sealed class CosineAndSineOverAPowerOfAnotherIntegralTest + { + [Theory] + [InlineData("sin(x)/(a + b*cos(x) + c*sin(x))")] + [InlineData("(k + q*cos(x))/(a + b*cos(x) + c*sin(x))")] + [InlineData("(k + p*cos(x) + q*sin(x))/(a + b*cos(x) + c*sin(x))^2")] + [InlineData("(p*cos(x) + q*sin(x))/(a + b*cos(x) + c*sin(x))^3")] + [InlineData("1/(a + b*cos(x) + c*sin(x))^2")] + public void IsWrittenThroughItsDenominator(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pinned(Entity e) => e.Substitute("a", 2.3).Substitute("b", 0.9).Substitute("c", 0.7) + .Substitute("k", 1.1).Substitute("p", 0.6).Substitute("q", 0.8); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { -2.6, -1.6, -0.4, 0.4, 1.4, 2.0 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + 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) + Math.Abs((double)want.ImaginaryPart)), + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + } +}