diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 6dc4bb0d0..a24ce827b 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1131,6 +1131,33 @@ quotients that nothing below reads. Rubi's 6.3.2 and 6.5.3 | `"sech(a+2*ln(c/x^(1/2)))^3".Integrate("x")` | left unevaluated | the antiderivative | | `"x*tanh(a+2*ln(x))^2".Integrate("x")` | left unevaluated | the antiderivative | +### An exponential times a trigonometric at the same frequency no longer divides by zero, and a phase is expanded + +`x e^(-i x) cos(x)` was answered `NaN` at every point of its domain. The closed form +`int e^(a x)(c cos(b x) + d sin(b x)) = e^(a x)(...)/(a^2 + b^2)` is guarded against the resonance +`a = ±i b`, where the product is a sum of two exponentials and the formula divides by zero -- but +the guard tested whether `a^2 + b^2` *evaluates* to a number, and with a symbolic frequency +`(-i f)^2 + f^2` is `-f^2 + f^2`, which `InnerSimplified` does not collect. So the guard passed, +the formula divided by a zero it could not see, and the answer was `NaN`: a wrong answer, not a +decline. Both places that test it -- the table pattern and the polynomial route -- see the +resonance symbolically now. + +With that fixed, two capabilities: a **phase** in the trigonometric's argument is expanded by the +angle-sum identity where an exponential stands beside it (`cos(f x + p)` is +`cos(p) cos(f x) - sin(p) sin(f x)`; the rule reads one frequency and no phase, an exponential's +own offset being a constant factor where a trigonometric's is not), and **`A + i A tan(z)` below +the bar** is written as `A e^(i z)/cos(z)`, with `A + i A cot(z)` as `i A e^(-i z)/sin(z)` -- +exactly, since `cos z + i sin z` is `e^(i z)`. Beside a polynomial that is the shape the closed +rule answers, where the imaginary unit in the coefficient is read by nothing else. Rubi's 4.3.10, +4.4.1 and 4.4.10 ([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x*e^(-i*x)*cos(x)".Integrate("x")` | left unevaluated | `x^2/4 + e^(-2 i x)(1 + 2 i x)/8` up to the form -- where `master` between the releases answered `NaN` | +| `"x*e^(2*x)*cos(3*x+1)".Integrate("x")` | left unevaluated | the antiderivative | +| `"x/(2+2*i*tan(x))".Integrate("x")` | left unevaluated | the antiderivative | +| `"(c+d*x)/(a+i*a*tan(pe+f*x))".Integrate("x")` | left unevaluated | the antiderivative | + ### `binomial(n, k)` is a function **Addition, not silent.** The binomial coefficient is a node, `Entity.Binomialf`, spelled diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 7c8cf4439..9dce20aa0 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -3359,7 +3359,10 @@ void Add(int k, bool isSine, Entity coefficient) var exponential = MathS.Pow(MathS.e, rate * x); var cosine = MathS.Cos(frequency * x); var sine = MathS.Sin(frequency * x); - if (scale == 0 || scale.Evaled is Number.Complex { IsZero: true }) + // Symbolically as well as numerically: with a symbolic frequency `f` the rate `-i f` + // leaves `-f^2 + f^2`, which `InnerSimplified` does not collect. + if (scale == 0 || scale.Evaled is Number.Complex { IsZero: true } + || scale.Vars.Any() && Functions.PartialFractions.Bare(scale.Simplify()).Evaled is Number.Complex { IsZero: true }) { // Resonant: by Euler, `c cos(bx) + d sin(bx)` is // `(c - i d) e^(i b x)/2 + (c + i d) e^(-i b x)/2`, and each term beside @@ -14479,6 +14482,134 @@ static bool IsAWholeNumberOfSteps(ERational difference) return (inFront * antiderivative - frontSlope * integralOfAntiderivative).InnerSimplified; } + /// + /// A sine or cosine of a linear form with a phase, beside an exponential, expanded by the + /// angle-sum identity so that both are of the same bare multiple of the variable: + /// cos(f x + p) is cos(p) cos(f x) - sin(p) sin(f x). The closed rule for a + /// polynomial times an exponential times a trigonometric reads one frequency and no + /// phase -- an exponential's own offset is a constant factor, a trigonometric's is not -- + /// and a phase is what every row of Rubi's 4.3.10 and 4.4.10 carries once + /// A + i A tan(pe + f x) is written as an exponential over a cosine. + /// + /// + /// Only beside an exponential of the variable, which is the shape that rule answers: the + /// expansion doubles the terms of every sine and cosine it touches, and the sum split + /// would then answer `sin(a + b x)` as two integrals where one closed rule answers it as + /// one. Once: the expanded form has no phase left. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveByExpandingATrigonometricPhaseBesideAnExponential(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + // An exponential of the variable among the factors, and a phase to expand. + if (!expr.Nodes.Any(node => node is Powf(var @base, var exponent) && !@base.ContainsNode(x) && exponent.ContainsNode(x))) + return null; + var expanded = expr.Replace(node => + { + Entity argument; + bool isSine; + switch (node) + { + case Sinf(var inner): (argument, isSine) = (inner, true); break; + case Cosf(var inner): (argument, isSine) = (inner, false); break; + default: return node; + } + if (!argument.ContainsNode(x) || !TreeAnalyzer.TryGetPolyLinear(argument, x, out var slope, out var offset) + || TreeAnalyzer.IsZero(slope) || offset.ContainsNode(x) || TreeAnalyzer.IsZero(offset) + || offset.Evaled is Number.Integer(0)) + return node; + var bare = (slope * x).InnerSimplified; + return isSine + ? MathS.Sin(offset) * MathS.Cos(bare) + MathS.Cos(offset) * MathS.Sin(bare) + : MathS.Cos(offset) * MathS.Cos(bare) - MathS.Sin(offset) * MathS.Sin(bare); + }); + if (expanded == expr) + return null; + return Integration.ComputeAsAQuestionOfItsOwn(expanded.Expand().InnerSimplified, x, integrateByParts); + } + + /// + /// A + i A tan(z) is A e^(i z)/cos(z), and A + i A cot(z) is + /// i A e^(-i z)/sin(z) -- exactly, wherever the tangent is defined, since + /// cos(z) + i sin(z) is e^(i z). A power of one of them below the bar is + /// then a power of the cosine times an exponential, and beside a polynomial that is a + /// shape the closed rules answer: (c + d x)^2/(a + i a tan(pe + f x))^2 is a + /// search past the budget as written and a second once rewritten. + /// + /// + /// The four signs are the four identities: A - i A tan(z) is + /// A e^(-i z)/cos(z) and A - i A cot(z) is -i A e^(i z)/sin(z). + /// Rubi's 4.3.10 and 4.4.10, where every such row carries the imaginary unit in the + /// coefficient and nothing else reads it. The same question, not one of its own: what is + /// handed on is the integrand in another spelling. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveByWritingAnImaginaryTangentAsAnExponential(Entity expr, Entity.Variable x, bool integrateByParts) + { + // Below the bar only: `(A + i A tan(z))^3` above it is a polynomial in the tangent, + // which the rules for those answer in the tangent and more shortly, where the + // reciprocal is what nothing reads. + var (above, below) = Functions.SingleQuotient.Of(Functions.SingleQuotient.Combine(expr)); + if (!below.ContainsNode(x)) + return null; + var rewritten = below.Replace(node => + { + if (node is not Sumf and not Minusf || !node.ContainsNode(x)) + return node; + Entity constant = Number.Integer.Zero; + Entity? coefficient = null; + Entity? argument = null; + var isTangent = false; + foreach (var term in Sumf.LinearChildren(node)) + { + if (!term.ContainsNode(x)) + { + constant += term; + continue; + } + if (coefficient is not null) + return node; + Entity factors = Number.Integer.One; + foreach (var factor in Mulf.LinearChildren(term)) + switch (factor) + { + case Tanf(var inner) when argument is null: + (argument, isTangent) = (inner, true); + break; + case Cotanf(var inner) when argument is null: + (argument, isTangent) = (inner, false); + break; + default: + if (factor.ContainsNode(x)) + return node; + factors *= factor; + break; + } + if (argument is null) + return node; + coefficient = factors; + } + if (coefficient is null || argument is null || constant == Number.Integer.Zero) + return node; + // The ratio is the imaginary unit one way or the other, and nothing else. + // Bare: `i a/a` simplifies to `i provided not a = 0`, and a condition is not a + // number to compare against. + var ratio = Functions.PartialFractions.Bare((coefficient / constant).InnerSimplified); + if (ratio.Evaled is not Number.Complex) + ratio = Functions.PartialFractions.Bare(ratio.Simplify()); + var plus = ratio.Evaled == MathS.i.Evaled; + if (!plus && ratio.Evaled != (-MathS.i).Evaled) + return node; + var exponential = MathS.Pow(MathS.e, ((plus == isTangent ? MathS.i : -MathS.i) * argument).InnerSimplified); + // A + i A tan(z) is A e^(iz)/cos(z); A + i A cot(z) is i A e^(-iz)/sin(z). + return isTangent + ? constant * exponential / MathS.Cos(argument) + : (plus ? MathS.i : -MathS.i) * constant * exponential / MathS.Sin(argument); + }); + return rewritten == below ? null : Integration.ComputeAsTheSameQuestion((above / rewritten).InnerSimplified, x, integrateByParts); + } + /// /// Bioche's first two rules for the hyperbolic functions: a rational function of /// sinh(y) and cosh(y) that is odd in the hyperbolic sine is a rational diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs index fc0de53c2..82052f481 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IntegralPatterns.cs @@ -420,8 +420,19 @@ private static bool IsZeroOnceSimplified(Entity expr) /// exponential times a sine or cosine divides by: zero for a rate of i times /// the frequency, where the product is a sum of two exponentials instead. /// + /// + /// Symbolically as well as numerically. With a symbolic frequency f, the rate + /// -i f -- which is what A + i A tan(f x) leaves once written as + /// A e^(i f x)/cos(f x) -- gives -f^2 + f^2, which + /// does not collect, so the evaluation is not a + /// number and the guard passed: the answer then divided by a zero it could not see and + /// was NaN at every point. + /// private static bool NotResonant(Entity rate, Entity frequency) - => (rate * rate + frequency * frequency).InnerSimplified.Evaled is not Entity.Number.Complex { IsZero: true }; + { + var scale = (rate * rate + frequency * frequency).InnerSimplified; + return scale.Evaled is not Entity.Number.Complex { IsZero: true } && !IsZeroOnceSimplified(scale); + } /// /// Whether is an exponential in , and diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 886a371e0..fafb4772a 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -501,6 +501,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // An exponential of a multiple of a logarithm is a power of the argument, which is // how every inverse hyperbolic function under an exponential arrives. if ((answer = IndefiniteIntegralSolver.SolveByFoldingAnExponentialOfALogarithm(expr, x, integrateByParts)) is { }) return answer; + // `A + i A tan(z)` is `A e^(i z)/cos(z)`, which beside a polynomial is a shape the + // closed rules answer, where the imaginary unit in the coefficient is read by none. + if ((answer = IndefiniteIntegralSolver.SolveByWritingAnImaginaryTangentAsAnExponential(expr, x, integrateByParts)) is { }) return answer; // A fractional power of a perfect square is the power of the modulus, sgn(P) P^(2r). if ((answer = IndefiniteIntegralSolver.SolveByTakingARootOfAPerfectSquare(expr, x, integrateByParts)) is { }) return answer; // x^(n - 1) g(x^n) with a symbolic n is g(u)/n under u = x^n. @@ -637,6 +640,9 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // rules, because none of them reads the exponential and all of them would have to // decline it. if ((answer = IndefiniteIntegralSolver.SolveAPolynomialTimesAnExponentialAndATrigonometric(expr, x)) is { }) return answer; + // And the same shape with a phase in the trigonometric's argument, expanded by the + // angle-sum identity so that the rule above reads one frequency and no phase. + if ((answer = IndefiniteIntegralSolver.SolveByExpandingATrigonometricPhaseBesideAnExponential(expr, x, integrateByParts)) is { }) return answer; // A half-integer power of `1 + sin` or `1 - cos` and their kin, closed by one // cancellation that `a^2 = b^2` allows. Before the reductions, which do not read a // fractional power at all. diff --git a/Sources/Tests/UnitTests/Calculus/PolynomialExponentialTrigTest.cs b/Sources/Tests/UnitTests/Calculus/PolynomialExponentialTrigTest.cs index 0cc4d0d43..e4cd459c3 100644 --- a/Sources/Tests/UnitTests/Calculus/PolynomialExponentialTrigTest.cs +++ b/Sources/Tests/UnitTests/Calculus/PolynomialExponentialTrigTest.cs @@ -205,5 +205,62 @@ public void TheSearchThisAvoidsOpening(string integrand) return; DifferentiatesBack(integrand); } + + /// + /// A phase in the trigonometric's argument, expanded by the angle-sum identity where an + /// exponential stands beside it: the closed rule reads one frequency and no phase, an + /// exponential's own offset being a constant factor where a trigonometric's is not. + /// + [Theory] + [InlineData("x*e^(2*x)*cos(3*x + 1)")] + [InlineData("x^2*e^x*sin(x + 2)")] + [InlineData("e^(3*x)*cos(2*x + 1/2)*sin(2*x + 1/2)")] + public void APhaseIsExpandedBesideAnExponential(string integrand) => DifferentiatesBack(integrand); + + /// + /// The resonant case, where the rate is the frequency times i and the closed + /// form's a^2 + b^2 is zero: the product is a sum of two exponentials instead. + /// With a symbolic frequency -f^2 + f^2 is not collected by + /// , so the guard read it as a nonzero number and the + /// answer divided by it -- NaN at every point, for every row of Rubi's 4.3.10 that + /// carries a + i a tan(pe + f x) below the bar. + /// #718 + /// + [Theory] + [InlineData("x*e^(-i*f*x)*cos(f*x)", "f")] + [InlineData("e^(i*f*x)*sin(f*x)", "f")] + [InlineData("x^2*e^(-i*x)*cos(x)", null)] + public void TheResonanceIsSeenSymbolically(string integrand, string? symbol) + { + var integral = integrand.ToEntity().Integrate("x").Substitute("C", 0); + Assert.DoesNotContain("NaN", integral.Stringize()); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pin(Entity e) => symbol is null ? e : e.Substitute(symbol, 1.3); + var derivative = Pin(integral).Differentiate("x"); + var original = Pin(integrand.ToEntity()); + foreach (var at in Points) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + Assert.False(got.IsNaN, $"the antiderivative of {integrand} differentiates to NaN at x = {at}"); + var difference = Math.Abs((double)(got - want).RealPart) + Math.Abs((double)(got - want).ImaginaryPart); + var scale = Math.Max(1.0, Math.Abs((double)want.RealPart) + Math.Abs((double)want.ImaginaryPart)); + Assert.True(difference / scale < 1e-9, + $"d/dx of the antiderivative of {integrand} is {got} at x = {at}, where the integrand is {want}"); + } + } + + /// + /// A + i A tan(z) below the bar is A e^(i z)/cos(z), and + /// A + i A cot(z) is i A e^(-i z)/sin(z): beside a polynomial that is the + /// shape above, where the imaginary unit in the coefficient is read by nothing else. + /// Rubi's 4.3.10 and 4.4.10. Above the bar the tangent's own rules answer it, and the + /// rewrite leaves it alone. + /// + [Theory] + [InlineData("x/(2 + 2*i*tan(x))")] + [InlineData("x^2/(3 - 3*i*cot(x))")] + [InlineData("(1 + 2*x)/(3 + 3*i*tan(1/2 + x))^2")] + public void AnImaginaryTangentBelowTheBarIsAnExponential(string integrand) => DifferentiatesBack(integrand); } }