diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 76e155998..4ca552def 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -541,6 +541,24 @@ tanh keep the answers the exponential substitution gives them ([#233](https://gi | `"sin(tanh(a+b*x))".ToEntity().Integrate("x")` | `integral(...)` | `Si` and `Ci` of `1 - tanh(a + b x)` and `1 + tanh(a + b x)` | | `"tanh(x)/sqrt(a+b*tanh(x)^2+c*tanh(x)^4)".ToEntity().Integrate("x")` | `integral(...)` | a logarithm of the root, piecewise on the sign of `a + b + c` | +### The derivative of a sum of exponentials over another is read off the integrand over the other's square + +**Answers where there were none.** Rubi's 6.7.1 `x (b - a sinh(x))/(a + b sinh(x))^2` and `x (b + a cosh(x))/(a + b cosh(x))^2` +were declined. The bracket over the square is the derivative of `-cosh(x)/(a + b sinh(x))`, and of +`sinh(x)/(a + b cosh(x))`. So by parts the integral is `x` times that quotient, plus `ln(a + b sinh(x))/b` or minus +`ln(a + b cosh(x))/b`. The integrals without x were answered with a piecewise on a discriminant, from the +substitution in `e^x`; they are now the quotient itself. For a sum `D` of exponentials, the rule reads the +integrand over `D^2` as `(N/D)'`, with `N` a constant times each exponential in `D` +([#233](https://github.com/asc-community/AngouriMath/issues/233)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"x*(b-a*sinh(x))/(a+b*sinh(x))^2".ToEntity().Integrate("x")` | `integral(...)` | 123 characters | +| `"x*(b+a*cosh(x))/(a+b*cosh(x))^2".ToEntity().Integrate("x")` | `integral(...)` | 119 characters | +| `"(b-a*sinh(x))/(a+b*sinh(x))^2".ToEntity().Integrate("x")` | `integral(...)` | 71 characters, where before this change master gave 7,861 | +| `"(b+a*cosh(x))/(a+b*cosh(x))^2".ToEntity().Integrate("x")` | `integral(...)` | 70 characters, where before this change master gave 7,631 | +| `"sech(x)^2*tanh(x)^2/(2+tanh(x)^3)^2".ToEntity().Integrate("x")` | no answer in thirty seconds | 60 characters, `-1/3/(2 + tanh(x)^3)`, in 55 ms; master gave 117 in 16 s and 9.7 GB | + ### A hyperbolic function of `a + 2 ln(x)` is integrated with `e^a` as a symbol too **Answers where there were none.** `tanh(a + 2 ln(x))^2/x^3` was declined, and so were nineteen more of Rubi's 6.3.2, diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 1ad367eb0..0e3f4d960 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -42689,6 +42689,162 @@ Entity Pinned(Entity e) } } + /// + /// The derivative of N/D, for D a sum of exponentials of x and constants and + /// N a constant plus a constant times each exponential in D, read off the + /// integrand over D^2: alone, or beside x by parts. Rubi's + /// x (b - a sinh(x))/(a + b sinh(x))^2 is x G' for G = -cosh(x)/(a + b sinh(x)), + /// so x G - integral G, and the integral of G is -ln(a + b sinh(x))/b. + /// + /// + /// + /// The exponentials' span is closed under the derivative, so (N/D)' = (N' D - N D')/D^2 + /// with N' D - N D' linear in N's constants, and matching it to the integrand's + /// numerator is a linear system: each exponential a whole power of the first one, w, + /// and the identity read at six values of w. The answer is differentiated back at sampled + /// points. + /// + /// + /// The split took the numerator's terms apart, and the exponential substitution answered each + /// over the square of a quadratic in e^x with a piecewise on its discriminant; beside x, + /// parts integrated that piecewise again and ran past the budget, where the trigonometric + /// twin is 's. + /// https://github.com/asc-community/AngouriMath/issues/233 + /// + /// + internal static Entity? SolveAsTheDerivativeOfAQuotientOverASumOfExponentials(Entity expr, Entity.Variable x, bool integrateByParts) + { + if (!Integration.AnsweringTheQuestionAskedOrOneBelow || expr is not (Mulf or Divf)) + return null; + // A constant, x at most once, a numerator and the square of one sum below the bar. + Entity constant = Number.Integer.One; + Entity numerator = Number.Integer.One; + Entity? sum = null; + var besideX = false; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + constant = underneath ? constant / factor : constant * factor; + else if (factor == x && !underneath && !besideX) + besideX = true; + else if (underneath && sum is null && factor is Powf(Sumf or Minusf, Number.Integer { EInteger: var two } ) squared && two.Equals(EInteger.FromInt32(2))) + sum = squared.Base; + else if (!underneath) + numerator = numerator == Number.Integer.One ? factor : numerator * factor; + else + return null; + } + bool OnlyThroughTheExponentials(Entity e) + => !e.Replace(node => node is Powf(var @base, var exponent) && @base == MathS.e && exponent.ContainsNode(x) ? Number.Integer.One : node).ContainsNode(x); + if (sum is null || !OnlyThroughTheExponentials(sum) || !OnlyThroughTheExponentials(numerator)) + return null; + // Each exponential a whole power w^n of the first one, w = e^(m + k x): `e^(-(c + d x))` is + // `w^(-1)` beside `e^(c + d x)`, with no `e^c` left to the system. + var exponentials = new List<(Entity Node, Entity Exponent, Entity Rate, int Power)>(); + foreach (var node in expr.Nodes) + if (node is Powf(var @base, var exponent) && @base == MathS.e && exponent.ContainsNode(x) + && !exponentials.Any(known => known.Node == node)) + { + if (!TreeAnalyzer.TryGetPolyLinear(exponent, x, out var rate, out var offset) || rate.ContainsNode(x) || offset.ContainsNode(x) + || Functions.PartialFractions.IsZeroAsAValue(rate)) + return null; + var multiple = exponentials.Count == 0 ? Number.Integer.One : Functions.PartialFractions.Bare((exponent / exponentials[0].Exponent).Simplify()); + if (multiple is not Number.Integer { EInteger: var power } || !power.CanFitInInt32()) + return null; + exponentials.Add((node, exponent, rate, power.ToInt32Unchecked())); + } + var inTheSum = exponentials.Where(exponential => sum.ContainsNode(exponential.Node)).ToList(); + if (inTheSum.Count is 0 or > 4) + return null; + // N = c_0 + sum c_i E_i: N' D - N D' is sum c_i (k_i E_i D - E_i D') - c_0 D'. N and N + c D + // have one derivative over D, so the system has a line of solutions; without c_0 first, + // which is the one of them with no constant where D has one: `-cosh(x)/(a + b sinh(x))`, + // whose integral is a logarithm, and not `-(a/b + e^x)/(a + b sinh(x))`, a constant away. + var derivativeOfTheSum = sum.Differentiate(x); + var columns = new List { -derivativeOfTheSum }; + foreach (var (node, _, rate, _) in inTheSum) + columns.Add(rate * node * sum - node * derivativeOfTheSum); + var target = constant * numerator; + if ((Above(withoutAConstant: true) ?? Above(withoutAConstant: false)) is not { } above) + return null; + Entity? At(Entity e, Number.Rational w) + { + var written = e.Replace(node => node is Powf(var @base, var exponent) && @base == MathS.e && exponent.ContainsNode(x) + && exponentials.FirstOrDefault(known => known.Node == node) is { Node: not null } known + ? MathS.Pow(w, Number.Integer.Create(known.Power)) + : node).InnerSimplified; + return written.ContainsNode(x) ? null : written; + } + var g = above / sum; + if (!Functions.PartialFractions.DerivativeHoldsAtSampledPoints(g, target / MathS.Pow(sum, 2), x)) + return null; + if (!besideX) + return g; + // Where N is a constant times D', the integral of G is that constant times ln(D): the + // constant read off the first exponential's coefficients, D being linear in them, and + // the answer checked. Asked as a question, `-cosh(x)/(a + b sinh(x))` spelled in + // exponentials went through the substitution in e^x and came back a piecewise. + var (firstNode, _, firstRate, _) = inTheSum[0]; + Entity WithTheExponentials(Entity e, Entity? one) + => e.Replace(node => node is Powf(var @base, var exponent) && @base == MathS.e && exponent.ContainsNode(x) + ? (node == one ? Number.Integer.One : Number.Integer.Zero) : node); + var inTheSumOfTheFirst = LowestOverTheSymbols((WithTheExponentials(sum, firstNode) - WithTheExponentials(sum, null)).InnerSimplified); + var inTheNumeratorOfTheFirst = LowestOverTheSymbols((WithTheExponentials(above, firstNode) - WithTheExponentials(above, null)).InnerSimplified); + Entity? answer = null; + if (!Functions.PartialFractions.IsZeroAsAValue(inTheSumOfTheFirst)) + { + var logarithmic = LowestOverTheSymbols(inTheNumeratorOfTheFirst / (firstRate * inTheSumOfTheFirst)); + var closed = (x * g - logarithmic * MathS.Ln(sum)).InnerSimplified; + if (Functions.PartialFractions.DerivativeHoldsAtSampledPoints(closed, expr, x)) + answer = closed; + } + if (answer is null) + { + if (Integration.ComputeAsAQuestionOfItsOwn(g, x, integrateByParts) is not { } integralOfG + || integralOfG.Nodes.Any(node => node == MathS.NaN || node is Integralf)) + return null; + var asked = (x * g - integralOfG).InnerSimplified; + if (Functions.PartialFractions.DerivativeHoldsAtSampledPoints(asked, expr, x)) + answer = asked; + } + return answer; + + // N, the system read at six values of w, or null where it has no solution. + Entity? Above(bool withoutAConstant) + { + var first = withoutAConstant ? 1 : 0; + var samples = new[] { Number.Rational.Create(2), Number.Rational.Create(3), Number.Rational.Create(5), + Number.Rational.Create(1, 2), Number.Rational.Create(1, 3), Number.Rational.Create(7) }; + var matrix = new Entity[samples.Length][]; + var rhs = new Entity[samples.Length]; + for (var j = 0; j < samples.Length; j++) + { + matrix[j] = new Entity[columns.Count - first]; + for (var k = first; k < columns.Count; k++) + if (At(columns[k], samples[j]) is { } entry) + matrix[j][k - first] = entry; + else + return null; + if (At(target, samples[j]) is not { } value) + return null; + rhs[j] = value; + } + if (!Functions.PartialFractions.TrySolveLinearWithSymbols(matrix, rhs, out var coefficients)) + return null; + Entity read = Number.Integer.Zero; + for (var k = 0; k < coefficients.Length; k++) + { + var coefficient = LowestOverTheSymbols(coefficients[k]); + if (Functions.PartialFractions.IsZeroAsAValue(coefficient)) + continue; + var column = k + first; + var term = column == 0 ? coefficient : coefficient == Number.Integer.One ? inTheSum[column - 1].Node : coefficient * inTheSum[column - 1].Node; + read = read == Number.Integer.Zero ? term : read + term; + } + return read == Number.Integer.Zero ? null : read; + } + } + /// /// x^m (a + b atan(c x))^p over a power of d (1 +- i c x), under x = tan(u)/c, where /// 1 + i tan(u) is e^(i u)/cos(u): the integrand in u is diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 1544c7c7b..54eb3cf8e 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -720,6 +720,11 @@ internal static bool ASumOverRootsWouldAnswer // `(a + a cos(u))^(3/2)/cos(u)^(5/4)`, and parts integrated `x (b + a sin(x))/(a + b sin(x))^2`'s // derivative into a function of tan(x/2) and ran out of time. if ((answer = IndefiniteIntegralSolver.SolveAsTheDerivativeOfASineRaisedFromNothingBesidePowersOfTheCosine(expr, x, integrateByParts)) is { }) return answer; + // And its hyperbolic twin, where the functions are exponentials: the derivative of a sum of + // them over another, read off the integrand over the other's square, alone or beside x by + // parts. The split and the exponential substitution answered `(b - a sinh(x))/(a + b sinh(x))^2` + // with a piecewise on a discriminant, and parts beside x ran past the budget integrating it. + if ((answer = IndefiniteIntegralSolver.SolveAsTheDerivativeOfAQuotientOverASumOfExponentials(expr, x, integrateByParts)) is { }) return answer; // A cosine beside a function of the sine with a power off the rationals in it, in the sine: // the substitution search declined `cos(x)/(sin(x) + sin(x)^sqrt(2))`. if ((answer = IndefiniteIntegralSolver.SolveACosineBesideAFunctionOfTheSineWithAnIrrationalPowerInIt(expr, x, integrateByParts)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/TheDerivativeOfAQuotientOverASumOfExponentialsIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/TheDerivativeOfAQuotientOverASumOfExponentialsIntegralTest.cs new file mode 100644 index 000000000..8e8408dbb --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/TheDerivativeOfAQuotientOverASumOfExponentialsIntegralTest.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 +{ + /// + /// The derivative of a sum of exponentials over another, read off the integrand over the other's + /// square: (b - a sinh(x))/(a + b sinh(x))^2 is the derivative of -cosh(x)/(a + b sinh(x)), + /// and beside x, by parts, Rubi's 6.7.1 x (b - a sinh(x))/(a + b sinh(x))^2 is + /// -x cosh(x)/(a + b sinh(x)) + ln(a + b sinh(x))/b. The split and the substitution in + /// e^x answered the first with a piecewise on a discriminant, and the second ran past the budget. + /// #233 + /// + [Trait("Area", "Calculus")] + public sealed class TheDerivativeOfAQuotientOverASumOfExponentialsIntegralTest + { + [Theory] + [InlineData("x*(b - a*sinh(x))/(a + b*sinh(x))^2")] + [InlineData("x*(b + a*cosh(x))/(a + b*cosh(x))^2")] + [InlineData("x*(b - a*sinh(c + d*x))/(a + b*sinh(c + d*x))^2")] + [InlineData("(b - a*sinh(x))/(a + b*sinh(x))^2")] + [InlineData("(b + a*cosh(c + d*x))/(a + b*cosh(c + d*x))^2")] + public void IsAnsweredWithoutCases(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Assert.DoesNotContain("piecewise", text); + Entity Pinned(Entity e) => e.Substitute("a", 1.3).Substitute("b", 0.7).Substitute("c", 0.4).Substitute("d", 0.9); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { -1.7, -0.6, 0.4, 1.3, 2.2 }) + { + var want = original.Substitute("x", at).EvalNumerical(); + var got = derivative.Substitute("x", at).EvalNumerical(); + Assert.False(got.IsNaN, $"d/dx of the antiderivative of {integrand} has no value at x = {at}, where the integrand is {want}"); + 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}"); + } + } + } +}