diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index a89f8e5df..e5e1db970 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1256,6 +1256,31 @@ here that writes one does: left inside, a rule below differentiated it, and | `"sqrt(a^2+2*a*b*x+b^2*x^2)*sqrt(c+pe*x+d*x^2)".Integrate("x")` | left unevaluated | the antiderivative | | `"(a^2+2*a*b*x+b^2*x^2)^(5/2)".Integrate("x")` | the antiderivative, through the substitution search | `sgn(a + b x) (a + b x)^6 |b|^5/(6 b)` up to the form | +### A power of the variable times a sine or cosine of a logarithm, and a power of a monomial + +`x^2 sin(a + b ln(c x^n))` and `(c x^n)^b` were both left as written. Two rules, each exact: + +**The trigonometric of a logarithm** is a closed form rather than a search. Two rounds of parts +close on the integrand -- the sine's remainder is the cosine's integral and the cosine's is the +sine's -- so the pair is *solved*: `int x^m sin(L) dx` is +`x^(m + 1)((m + 1) sin(L) - B cos(L))/((m + 1)^2 + B^2)` wherever `L' = B/x` for a constant `B`, +which `a + b ln(c x^n)` is with `B = b n`; the cosine's is the same with +`(m + 1) cos(L) + B sin(L)`. The hyperbolic twin needs no rule -- the library writes `sinh` as +exponentials, which fold against the logarithm -- and the sine and cosine are nodes that fold +against nothing. + +**A power of a monomial** is distributed: `(c x^n)^b` is `c^b x^(n b)`, for a positive `c`, and +only where `n` is not whole. With a whole `n` the integrand is real at a negative `x` too, and +there `(c x^n)^b` is a power of `|x|`, not of `x`: `(2 u^3)^(3/2)` is `2^(3/2) sgn(u) u^(9/2)`, +which the rules for a root of an even power already answer with its sign +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x^2*sin(a+b*ln(c*x^n))".Integrate("x")` | left unevaluated | `x^3(3 sin(L) - b n cos(L))/(9 + b^2 n^2)` | +| `"cos(a+b*ln(c*x^n))".Integrate("x")` | left unevaluated | `x(cos(L) + b n sin(L))/(1 + b^2 n^2)` | +| `"(c*x^n)^b".Integrate("x")` | left unevaluated | `c^b x^(n b + 1)/(n b + 1) provided c > 0` | + ### `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 ec8ebcba7..8f08b74cd 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -14583,6 +14583,154 @@ private static (Entity Factor, Entity Inside)? TryTakeACommonLinearFactorOfASum( return back.Nodes.Any(node => node == MathS.NaN) ? null : back; } + /// + /// A power of the variable times a sine or cosine of a logarithm, in closed form: + /// int x^m sin(L) dx is x^(m + 1)((m + 1) sin(L) - B cos(L))/((m + 1)^2 + B^2) + /// whenever L' = B/x for a constant B, which a + b ln(c x^n) is with + /// B = b n; the cosine's is the same with (m + 1) cos(L) + B sin(L). + /// + /// + /// Two rounds of parts close on the integrand -- the sine's remainder is the cosine's + /// integral and the cosine's is the sine's -- and the pair is solved rather than + /// iterated. Differentiating the answer is the whole proof: the cross terms cancel and + /// what is left is x^m sin(L)((m + 1)^2 + B^2). The hyperbolic twin needs no rule, + /// the library writing sinh as exponentials that fold against the logarithm; + /// the sine and cosine are nodes and fold against nothing. Rubi's 4.7.5. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveAPowerTimesATrigonometricOfALogarithm(Entity expr, Entity.Variable x, bool integrateByParts) + { + Entity coefficient = Number.Integer.One; + Entity? degree = null; + Entity? argument = null; + var isSine = false; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + coefficient = underneath ? coefficient / factor : coefficient * factor; + continue; + } + switch (factor) + { + case Sinf(var inner) when argument is null && !underneath: + (argument, isSine) = (inner, true); + continue; + case Cosf(var inner) when argument is null && !underneath: + (argument, isSine) = (inner, false); + continue; + case Variable bare when bare == x && degree is null: + degree = underneath ? Number.Integer.Create(-1) : Number.Integer.One; + continue; + case Powf(Variable bare, var power) when bare == x && degree is null && !power.ContainsNode(x): + degree = underneath ? (-power).InnerSimplified : power; + continue; + default: + return null; + } + } + if (argument is null || !argument.ContainsNode(x)) + return null; + // L' x is the constant B, which is what makes the pair close. + var b = (argument.Differentiate(x) * x).InnerSimplified; + b = Functions.PartialFractions.Bare(b.Simplify()); + if (b.ContainsNode(x) || b.Evaled is Number.Complex { IsZero: true }) + return null; + var m = degree ?? Number.Integer.Zero; + var mPlusOne = (m + Number.Integer.One).InnerSimplified; + // `(m + 1)^2 + B^2`, which the answer divides by. Symbolically as well as + // numerically: for `B = n sqrt(-9/n^2)` and `m = 2` it is `9 + (n sqrt(-9/n^2))^2`, + // which `InnerSimplified` leaves standing and `Simplify` folds to zero -- an + // imaginary `B` of the right size makes the pair of parts circular rather than + // closing, and the answer divided by nothing at all. + var scale = (MathS.Sqr(mPlusOne) + MathS.Sqr(b)).InnerSimplified; + if (scale.Evaled is Number.Complex { IsZero: true } + || scale.Vars.Any() && Functions.PartialFractions.Bare(scale.Simplify()).Evaled is Number.Complex { IsZero: true }) + return null; + var sine = MathS.Sin(argument); + var cosine = MathS.Cos(argument); + var bracket = isSine ? mPlusOne * sine - b * cosine : mPlusOne * cosine + b * sine; + var answer = coefficient * MathS.Pow(x, mPlusOne) * bracket / scale; + return answer.InnerSimplified; + } + + /// + /// A power of a monomial with the power distributed: (c x^n)^b is + /// c^b x^(n b), which the power rule answers where the monomial as written is read + /// by nothing -- x^2 sin(a + b ln(c x^n)) folds to a power of c x^n and + /// stops there. + /// + /// + /// Exact on the domain the integrand has. With a symbolic n, x^n is real + /// only for a positive x -- at a negative one it is a complex number or nothing -- + /// so the integrand is defined there and nowhere else, and on x > 0 the + /// identity holds for a positive c. That condition travels with the answer as + /// provided c > 0, unless c is a positive number already; the + /// (e x)^(n - 1) of Rubi's 6.5.2 carries the same one. A whole exponent is + /// 's, and needs no condition. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveByDistributingAPowerOfAMonomial(Entity expr, Entity.Variable x, bool integrateByParts) + { + Entity assumed = Entity.Boolean.True; + var changed = false; + var written = expr.Replace(node => + { + if (node is not Powf(var @base, var exponent) || exponent.ContainsNode(x) || exponent is Number.Integer + || !@base.ContainsNode(x)) + return node; + // The base as `c x^n`, with `c` and `n` free of the variable and `n` not whole -- + // a whole one is the distributing rule's, which owes no condition. + Entity constant = Number.Integer.One; + Entity? degree = null; + foreach (var factor in Mulf.LinearChildren(@base)) + { + if (!factor.ContainsNode(x)) + { + constant *= factor; + continue; + } + if (degree is not null) + return node; + switch (factor) + { + case Variable bare when bare == x: + degree = Number.Integer.One; + break; + case Powf(Variable bare, var power) when bare == x && !power.ContainsNode(x): + degree = power; + break; + default: + return node; + } + } + // A whole `n` only where the identity needs nothing: for a negative `x`, + // `x^n` is a real number when `n` is whole, and `(c x^n)^b` is then `|x|`'s + // power, not `x`'s -- `(2 u^3)^(3/2)` is `2^(3/2) sgn(u) u^(9/2)`, and + // distributing it without the sign is a wrong answer where `u < 0`, which the + // rules for a root of an even power answer correctly and this must leave to + // them. With a fractional or symbolic `n` the integrand is real only for a + // positive `x`, and there the distribution is exact. + if (degree is null || degree is Number.Integer || degree.Evaled is Number.Integer) + return node; + if (constant != Number.Integer.One && constant.Evaled is not Number.Real { IsPositive: true }) + { + if (constant.Vars.Any() is false) + return node; + var positive = new Greaterf(constant, Number.Integer.Zero); + assumed = assumed == Entity.Boolean.True ? positive : assumed & positive; + } + changed = true; + var distributed = MathS.Pow(x, (degree * exponent).InnerSimplified); + return constant == Number.Integer.One ? distributed : MathS.Pow(constant, exponent) * distributed; + }); + if (!changed || written == expr) + return null; + if (Integration.ComputeAsTheSameQuestion(written.InnerSimplified, x, integrateByParts) is not { } answer) + return null; + return assumed == Entity.Boolean.True ? answer : answer.Provided(assumed); + } + /// /// A pair of hyperbolic powers two apart whose coefficients kill the reduction's /// residual: cosh(y)^p - (p - 1)/p cosh(y)^(p - 2) is diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 2f7c62191..6df2872e3 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -511,6 +511,12 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => // A whole power of a product of a constant and the variable, as the product of // the powers, which is how the inverse hyperbolic secant and cosecant arrive. if ((answer = IndefiniteIntegralSolver.SolveByDistributingWholePowersOfProducts(expr, x, integrateByParts)) is { }) return answer; + // And a fractional or symbolic power of a monomial: `(c x^n)^b` is `c^b x^(n b)` + // for a positive `c`, on the `x > 0` where a symbolic `n` leaves the integrand real. + if ((answer = IndefiniteIntegralSolver.SolveByDistributingAPowerOfAMonomial(expr, x, integrateByParts)) is { }) return answer; + // A power of the variable times a sine or cosine of a logarithm, where two rounds of + // parts close on the integrand: `int x^m sin(a + b ln(c x^n))` in closed form. + if ((answer = IndefiniteIntegralSolver.SolveAPowerTimesATrigonometricOfALogarithm(expr, x, integrateByParts)) is { }) return answer; // A logarithm of a quotient that cancels with the functions in it as // indeterminates, which is how an inverse hyperbolic function of a hyperbolic // one arrives: atanh(tanh(u)) is 1/2 ln(e^(2u)) once its quotient is cancelled. diff --git a/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs index f2a3f0490..804a4638a 100644 --- a/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs +++ b/Sources/Tests/UnitTests/Calculus/ExponentialOfALogarithmIntegralTest.cs @@ -90,5 +90,70 @@ private static void DifferentiatesBack(string integrand) [InlineData("e^(1 + ln(x + 2)/2)")] [InlineData("tanh(ln(x))")] public void AHyperbolicFunctionOfALogarithm(string integrand) => DifferentiatesBack(integrand); + + /// + /// A power of the variable times a sine or cosine of a logarithm, in closed form: two + /// rounds of parts close on the integrand, since the sine's remainder is the cosine's + /// integral and the cosine's is the sine's, and the pair is solved rather than iterated. + /// int x^m sin(L) is x^(m+1)((m+1) sin(L) - B cos(L))/((m+1)^2 + B^2) + /// wherever L' = B/x. The hyperbolic twin needs no rule, `sinh` being written as + /// exponentials that fold against the logarithm; the sine and cosine are nodes. + /// Rubi's 4.7.5. + /// + [Theory] + [InlineData("x^2*sin(a + b*ln(c*x^n))", "a=0.4,b=1.3,c=1.7,n=2.1")] + [InlineData("cos(a + b*ln(c*x^n))", "a=0.4,b=1.3,c=1.7,n=2.1")] + [InlineData("x^m*cos(a + b*ln(c*x^n))", "a=0.4,b=1.3,c=1.7,n=2.1,m=1.4")] + [InlineData("sin(2 + 3*ln(x))", "")] + [InlineData("x^2*sin(a + b*ln(x))/x", "a=0.4,b=1.3")] + public void APowerTimesATrigonometricOfALogarithm(string integrand, string pins) + { + var integral = integrand.ToEntity().Integrate("x").Substitute("C", 0); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.DoesNotContain("NaN", integral.Stringize()); + Entity Pin(Entity e) + { + foreach (var pin in pins.Split(',', StringSplitOptions.RemoveEmptyEntries)) + { + var parts = pin.Split('='); + e = e.Substitute(parts[0], double.Parse(parts[1], System.Globalization.CultureInfo.InvariantCulture)); + } + return e; + } + var derivative = Pin(integral).Differentiate("x"); + var original = Pin(integrand.ToEntity()); + foreach (var at in new[] { 0.3, 0.7, 1.1, 1.9, 2.6 }) + { + 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 fractional or symbolic power of a monomial, distributed: (c x^n)^b is + /// c^b x^(n b) on the x > 0 where a symbolic n leaves the + /// integrand real, for a positive c -- which the answer says. + /// + [Fact] + public void APowerOfAMonomialIsDistributed() + { + var integral = "(c*x^n)^b".ToEntity().Integrate("x").Substitute("C", 0); + Assert.DoesNotContain("integral(", integral.Stringize()); + Assert.Contains(integral.Nodes, node => node is Entity.Providedf(_, var predicate) && predicate == "c > 0".ToEntity()); + var pinned = integral.Substitute("c", 1.7).Substitute("n", 2.1).Substitute("b", 0.6); + var derivative = pinned.Differentiate("x"); + var original = "(c*x^n)^b".ToEntity().Substitute("c", 1.7).Substitute("n", 2.1).Substitute("b", 0.6); + foreach (var at in new[] { 0.3, 0.7, 1.1, 1.9, 2.6 }) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + Assert.True(Math.Abs((double)(got - want).RealPart) < 1e-9, $"at x = {at}: {got} for {want}"); + } + } } }