diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 0f440ae58..070dcb02e 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -620,6 +620,22 @@ are right on both sides of zero ([#718](https://github.com/asc-community/Angouri | `"x^2/(1 + x^4)^(3/4)".ToEntity().Integrate("x")` | `integral(...)` | logarithms and an arctangent of `(1 + x^4)^(1/4)/x` | | `"x^6*(3 + 4*x^4)^(1/4)".ToEntity().Integrate("x")` | `integral(...)` | the same in `(3 + 4 x^4)^(1/4)/x` | +### A power of a cosine and a sine plus their amplitude is integrated + +**Answers where there were none.** `sqrt(5 + 4 cos(x) + 3 sin(x))` was declined, with the rest of +Rubi's 4.7.7 powers of `a + b cos(y) + c sin(y)` with `a^2 = b^2 + c^2`, several after searches past +the budget. `b cos(y) + c sin(y)` is `R cos(y - phi)` for `R = sqrt(b^2 + c^2)`, so the base is +`a (1 ± cos(y - phi))`, a square of the half angle; a half-odd power of it, or a negative whole one, +is integrated in closed form now, through `T = b sin(y) - c cos(y)`, with no `phi` in the answer +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sqrt(5 + 4*cos(x) + 3*sin(x))".ToEntity().Integrate("x")` | `integral(...)` | `2 (4 sin(x) - 3 cos(x))/sqrt(5 + 4 cos(x) + 3 sin(x))` | +| `"1/sqrt(5 + 4*cos(x) + 3*sin(x))".ToEntity().Integrate("x")` | `integral(...)` | `ln((1 + q)/(1 - q))/sqrt(10)` for `q = (4 sin(x) - 3 cos(x))/(sqrt(10) sqrt(5 + 4 cos(x) + 3 sin(x)))` | +| `"(-5 + 4*cos(x) + 3*sin(x))^(3/2)".ToEntity().Integrate("x")` | `integral(...)` | `T sqrt(S)/(3/2) - (40/3) T/sqrt(S)` for `T = 4 sin(x) - 3 cos(x)` and `S` the base, imaginary as the integrand is | +| `"1/(b*cos(g + f*x) + c*sin(g + f*x) - sqrt(b^2 + c^2))^3".ToEntity().Integrate("x")` | `integral(...)` | `b sin(g + f x) - c cos(g + f x)` times powers of the base | + ### A rational function of the sine and cosine with symbols in it is integrated by the half angle `sin(x)^2/(a + b cos(x))` was left unevaluated while `sin(x)^2/(2 + 3 cos(x))` was answered. Under diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index c88422e78..120c9f71b 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -3550,6 +3550,148 @@ Entity Step(ERational power) return Step(n).InnerSimplified; } + /// + /// (a + b cos(y) + c sin(y))^n for a^2 = b^2 + c^2, y linear in + /// x and n half-odd or a negative whole number, in closed form. + /// + /// + /// + /// b cos(y) + c sin(y) is R cos(y - phi) for R = sqrt(b^2 + c^2), so + /// the base is a (1 ± cos(y - phi)), a square of the half angle as + /// a + a sin(y) is, and 's identity + /// holds for it. Rubi's sqrt(5 + 4 cos(x) + 3 sin(x)) and the rest of its 4.7.7 + /// with a^2 = b^2 + c^2 were declined or past the budget. The answer needs no + /// phi: with S the base and T = b sin(y) - c cos(y), + /// T' = S - a, S' = -T and T^2 = S (2a - S), the last of which is + /// a^2 = b^2 + c^2, and from them + /// + /// + /// d/dy (T S^(n-1)) = n S^n - a (2n - 1) S^(n-1) + /// + /// + /// which steps n down to 1/2, where the integral is 2T/sqrt(S), or + /// up to -1/2 or -1, where it is + /// (2/sqrt(2a)) atanh(T/(sqrt(2a) sqrt(S))) and T/(a S). Each base is + /// checked by differentiating it with those three identities alone. + /// + /// + /// For a real a of either sign: the derivations use only + /// S^(3/2) = S sqrt(S) and (sqrt(2a) sqrt(S))^2 = 2a S, which hold for + /// the principal powers, and the argument of the hyperbolic arctangent stays between + /// -1 and 1, one less its square being S/(2a), which is not negative + /// since S has the sign of a. An antiderivative on every interval between + /// the zeros of S, where T changes sign. a^2 = b^2 + c^2 is decided, + /// not assumed, and the cosine and the sine are both there: one alone is + /// 's. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + /// + internal static Entity? SolveAPowerOfACosineAndASinePlusTheirAmplitude(Entity expr, Entity.Variable x) + { + // A constant times one power of the base, the power above the bar or below it. + if (!TryReadAsQuotient(expr, out var above, out var below)) + (above, below) = (expr, Number.Integer.One); + Entity constant = Number.Integer.One; + Entity? @base = null; + ERational? exponent = null; + foreach (var (side, sign) in new[] { (above, 1), (below, -1) }) + foreach (var factor in Mulf.LinearChildren(side)) + { + if (!factor.ContainsNode(x)) + { + constant = sign > 0 ? constant * factor : constant / factor; + continue; + } + if (@base is not null) + return null; + (@base, exponent) = factor is Powf(var inner, var power) && power.Evaled is Number.Rational rational + ? (inner, rational.ERational) + : (factor, ERational.One); + if (sign < 0) + exponent = exponent.Negate(); + } + if (@base is null || exponent is null) + return null; + var n = exponent; + // Half-odd, or a negative whole number, and of a modest size: a positive whole power + // is a polynomial in the sine and cosine. + var isHalfOdd = n.Denominator.Equals(EInteger.FromInt32(2)); + if (!isHalfOdd && !(n.Denominator.Equals(EInteger.One) && n.Sign < 0) + || n.Abs().CompareTo(ERational.FromInt32(MaximumAmplitudePower)) > 0) + return null; + // a + b cos(y) + c sin(y), read off the sum. + Entity a = Number.Integer.Zero; + Entity? b = null, c = null, argument = null; + foreach (var term in Sumf.LinearChildren(@base)) + { + if (!term.ContainsNode(x)) + { + a = a == Number.Integer.Zero ? term : a + 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 null; + var inner = function?.DirectChildren.First(); + if (inner is null || argument is not null && inner != argument) + return null; + argument = inner; + if (function is Cosf) + { + if (b is not null) return null; + b = coefficient; + } + else + { + if (c is not null) return null; + c = coefficient; + } + } + if (b is null || c is null || argument is null || a == 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; + // a^2 = b^2 + c^2, decided rather than assumed. + var excess = (MathS.Sqr(a) - MathS.Sqr(b) - MathS.Sqr(c)).InnerSimplified; + if (excess.Evaled is not Number.Complex { IsZero: true } + && !(excess.Vars.Any() && Functions.PartialFractions.Bare(excess.Simplify()).Evaled is Number.Complex { IsZero: true })) + return null; + if (a.Evaled is Number.Complex { IsZero: true }) + return null; + + var s = @base; + var t = b * MathS.Sin(argument) - c * MathS.Cos(argument); + var half = ERational.Create(1, 2); + Entity Integral(ERational power) + { + if (power.CompareTo(half) == 0) + return 2 * t / MathS.Sqrt(s); + if (power.CompareTo(half.Negate()) == 0) + return 2 / MathS.Sqrt(2 * a) * MathS.Hyperbolic.Artanh(t / (MathS.Sqrt(2 * a) * MathS.Sqrt(s))); + if (power.CompareTo(ERational.FromInt32(-1)) == 0) + return t / (a * s); + var current = Number.Rational.Create(power); + if (power.Sign > 0) + return t * MathS.Pow(s, Number.Rational.Create(power.Subtract(ERational.One))) / current + + a * (2 * current - 1) / current * Integral(power.Subtract(ERational.One)); + return ((current + 1) * Integral(power.Add(ERational.One)) - t * MathS.Pow(s, current)) / (a * (2 * current + 1)); + } + return (constant * Integral(n) / rate).InnerSimplified; + } + + /// + /// The largest power, either way, that + /// steps through: each step is a term of the answer. + /// + private const int MaximumAmplitudePower = 8; + /// /// Fractional powers of a ± a sin(y), or of a ± a cos(y), beside anything /// rational in the sine and cosine of y, by the half angle at which they are diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 94ed6ccd9..9f1119d6b 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -858,6 +858,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // cosine, by the half angle at which they are squares: `1 + sin(y)` is `2 sin(u)^2`. Before // the substitution search, which spent twenty seconds on the radicals of the sine. if ((answer = IndefiniteIntegralSolver.SolveByTheHalfAngleWhereOnePlusASineIsASquare(expr, x, integrateByParts)) is { }) return answer; + // And a power of `a + b cos(y) + c sin(y)` with `a^2 = b^2 + c^2`, the same square + // turned by a phase, in closed form: the substitution search spent the budget on it. + if ((answer = IndefiniteIntegralSolver.SolveAPowerOfACosineAndASinePlusTheirAmplitude(expr, x)) is { }) return answer; // And symbolic powers of `a ± a sin(y)` beside a power of `g cos(y)`, by u = sin(y), where // `(1 + u)(1 - u)` is the cosine's square: the half angle wants numeric powers, and the // substitution search spent the budget on these. diff --git a/Sources/Tests/UnitTests/Calculus/CosineAndSinePlusTheirAmplitudeIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/CosineAndSinePlusTheirAmplitudeIntegralTest.cs new file mode 100644 index 000000000..04f0c47c7 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/CosineAndSinePlusTheirAmplitudeIntegralTest.cs @@ -0,0 +1,54 @@ +// +// 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 power of a + b cos(y) + c sin(y) with a^2 = b^2 + c^2, which is + /// a (1 ± cos(y - phi)), a square of the half angle: half-odd powers either way and + /// negative whole ones, in closed form through T = b sin(y) - c cos(y). Rubi's 4.7.7. + /// #718 + /// + /// + /// Compared as complex numbers on both sides of zero and away from the zeros of the base: + /// with a negative the base is not positive anywhere, and the integrand is imaginary + /// on the whole line. + /// + [Trait("Area", "Calculus")] + public sealed class CosineAndSinePlusTheirAmplitudeIntegralTest + { + [Theory] + [InlineData("sqrt(5 + 4*cos(x) + 3*sin(x))")] + [InlineData("(5 + 4*cos(x) + 3*sin(x))^(5/2)")] + [InlineData("1/sqrt(5 + 4*cos(x) + 3*sin(x))")] + [InlineData("1/(5 + 4*cos(x) + 3*sin(x))^(3/2)")] + [InlineData("(-5 + 4*cos(x) + 3*sin(x))^(3/2)")] + [InlineData("1/(b*cos(g + f*x) + c*sin(g + f*x) + sqrt(b^2 + c^2))^(5/2)")] + [InlineData("(b*cos(g + f*x) + c*sin(g + f*x) - sqrt(b^2 + c^2))^(3/2)")] + [InlineData("1/(b*cos(g + f*x) + c*sin(g + f*x) - sqrt(b^2 + c^2))^3")] + public void IsWrittenThroughTheDerivativeOfItsBase(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pinned(Entity e) => e.Substitute("b", 0.9).Substitute("c", 0.7).Substitute("g", 0.4).Substitute("f", 1.3); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { -1.9, -1.1, -0.4, 0.4, 1.1, 1.9 }) + { + 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}"); + } + } + } +}