diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 757f04982..f964211c6 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -618,6 +618,19 @@ interval without a case for each. Only where nothing else answers the integrand | `"acoth(x)^2/(1-x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 1,304 characters, with `(ln((x + 1)/(x - 1)) - ln((1 + x)/(1 - x)))/2` for `K` | | `"acoth(c*x)^2/(1-c^2*x^2)^2".ToEntity().Integrate("x")` | `integral(...)` | 1,681 characters | +### Half-odd powers of a hyperbolic sine beside an exponential are read as a derivative + +**An answer where there was none.** Rubi's 6.7.1 `e^(a + b x) sinh(c + d x)^(3/2) - 3/4 e^(a + b x) d^2/((b^2 - 9/4 d^2) sqrt(sinh(c + d x)))` +was declined. Each of its terms alone is not elementary, and the sum is the derivative of +`e^(a + b x) (4 b sinh^(3/2) - 6 d cosh sinh^(1/2))/(4 b^2 - 9 d^2)`. Such a sum is now read as the derivative of +`e^(a + b x)` times a sum of `sinh^p cosh^q`, over the integrand's monomials and the ones their derivative reaches, with +`cosh^2 = 1 + sinh^2`. Its coefficients are a linear system +([#718](https://github.com/asc-community/AngouriMath/issues/718)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"e^(a+b*x)*sinh(c+d*x)^(3/2)-3/4*e^(a+b*x)*d^2/((b^2-9/4*d^2)*sqrt(sinh(c+d*x)))".ToEntity().Integrate("x")` | `integral(...)` | 249 characters | + ### A function of the hyperbolic tangent alone is integrated under `u = tanh` **Answers where there were none.** `tanh(x)/sqrt(a + b tanh(x)^2 + c tanh(x)^4)` was declined, and so were 35 more of diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index dc44a289e..2ad189101 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -29768,6 +29768,148 @@ node is Powf(var @base, var exponent) && @base == MathS.e && powers.TryGetValue( return Folded(written.Substitute(constant, named - s * slope * x)).Substitute(s, sign); } + /// + /// An exponential of a linear times half-odd powers of a hyperbolic sine of another linear, and the + /// cosine to the first power beside them, as the derivative of the same kind of sum: with + /// E = e^(a + b x), S = sinh(u), C = cosh(u) and u = c + d x, + /// (E S^p C^q)' = E (b S^p C^q + d p S^(p - 1) C^(q + 1) + d q S^(p + 1) C^(q - 1)), with + /// C^2 = 1 + S^2, so F is a sum of E S^p C^q over the integrand's monomials and their + /// partners, and its coefficients a linear system. Rubi's 6.7.1 + /// e^(a + b x) sinh(c + d x)^(3/2) - 3/4 e^(a + b x) d^2/((b^2 - 9/4 d^2) sqrt(sinh(c + d x))) + /// is (e^(a + b x) (b sinh^(3/2) - 3/2 d cosh sinh^(1/2))/(b^2 - 9/4 d^2))', where each term + /// alone is not elementary and the split declined both. Powers of one sine multiply on every + /// branch, and no power of a power is taken. Checked by its derivative. + /// https://github.com/asc-community/AngouriMath/issues/718 + /// + internal static Entity? SolveHalfOddPowersOfAHyperbolicSineBesideAnExponentialAsADerivative(Entity expr, Entity.Variable x) + { + if (!Integration.AnsweringTheQuestionAsked) + return null; + var terms = Sumf.LinearChildren(expr); + if (terms.Count < 2) + return null; + // u: the argument of the first e^u - e^(-u) or e^u + e^(-u) in the integrand. + Entity? u = null; + foreach (var node in expr.Nodes) + if (node is Minusf or Sumf && node.DirectChildren.ToList() is [Powf(var e1, var p1), Powf(var e2, var p2)] + && e1 == MathS.e && e2 == MathS.e && TreeAnalyzer.TryGetPolyLinear(p1, x, out var s1, out var c1) && !TreeAnalyzer.IsZero(s1) + && TreeAnalyzer.TryGetPolyLinear(p2, x, out var s2, out var c2) + && Functions.PartialFractions.InLowestTermsOverTheSymbols(s1 + s2) == Number.Integer.Zero + && Functions.PartialFractions.InLowestTermsOverTheSymbols(c1 + c2) == Number.Integer.Zero) + { + u = p1; + break; + } + if (u is null || !TreeAnalyzer.TryGetPolyLinear(u, x, out var d, out _)) + return null; + bool IsPlusMinus(Entity first, Entity second) => first is Powf(var e1, var p1) && e1 == MathS.e && p1 == u + && second is Powf(var e2, var p2) && e2 == MathS.e && Functions.PartialFractions.InLowestTermsOverTheSymbols(p2 + u) == Number.Integer.Zero; + Entity? exponentOfE = null; + // (constant, power of S, power of C), the exponential's exponent kept apart. + (Entity K, ERational P, int Q)? Read(Entity node) + { + if (!node.ContainsNode(x)) + return (node, ERational.Zero, 0); + switch (node) + { + case Minusf(var a, var b) when IsPlusMinus(a, b): + return (Number.Integer.Create(2), ERational.One, 0); + case Sumf(var a, var b) when IsPlusMinus(a, b): + return (Number.Integer.Create(2), ERational.Zero, 1); + case Powf(var e0, var exponent) when e0 == MathS.e && TreeAnalyzer.TryGetPolyLinear(exponent, x, out _, out _): + if (exponentOfE is null) + exponentOfE = exponent; + return exponentOfE == exponent ? (Number.Integer.One, ERational.Zero, 0) : null; + case Mulf(var a, var b): + return Read(a) is var (ka, pa, qa) && Read(b) is var (kb, pb, qb) ? (ka * kb, pa.Add(pb), qa + qb) : null; + case Divf(var a, var b): + return Read(a) is var (ka2, pa2, qa2) && Read(b) is var (kb2, pb2, qb2) ? (ka2 / kb2, pa2.Subtract(pb2), qa2 - qb2) : null; + // A power of the sine alone, or a whole power of anything read. + case Powf(var @base, Number.Integer { EInteger: var whole }) when whole.CanFitInInt32() && Read(@base) is var (k, pw, qw): + return (MathS.Pow(k, Number.Integer.Create(whole)), pw.Multiply(ERational.FromEInteger(whole)), qw * whole.ToInt32Unchecked()); + case Powf(var @base, Number.Rational { ERational: var r }) when Read(@base) is var (k2, ps, qs) && qs == 0 && ps.Equals(ERational.One): + return (MathS.Pow(k2, Number.Rational.Create(r)), r, 0); + default: + return null; + } + } + // The integrand's monomials, C^2 written 1 + S^2. + var target = new Dictionary<(ERational, int), Entity>(); + void Add(Dictionary<(ERational, int), Entity> into, ERational p, int q, Entity coefficient) + { + while (q >= 2) + { + Add(into, p.Add(ERational.FromInt32(2)), q - 2, coefficient); + q -= 2; + } + if (q < 0) + return; + into[(p, q)] = into.TryGetValue((p, q), out var already) ? already + coefficient : coefficient; + } + var halfOdd = false; + foreach (var term in terms) + { + if (Read(term) is not var (k, p, q) || q < 0 || q > 3) + return null; + halfOdd |= !p.Denominator.Equals(EInteger.One); + Add(target, p, q, k); + } + if (!halfOdd || exponentOfE is null || !TreeAnalyzer.TryGetPolyLinear(exponentOfE, x, out var b, out _)) + return null; + // The ansatz: every monomial and its partner, which the derivative reaches. + var basis = new List<(ERational P, int Q)>(); + foreach (var (p, q) in target.Keys) + foreach (var candidate in q == 0 ? new[] { (p, 0), (p.Subtract(ERational.One), 1) } : new[] { (p, 1), (p.Add(ERational.One), 0) }) + if (!basis.Contains(candidate)) + basis.Add(candidate); + var derivatives = basis.Select(m => + { + var into = new Dictionary<(ERational, int), Entity>(); + Add(into, m.P, m.Q, b); + if (!m.P.IsZero) + Add(into, m.P.Subtract(ERational.One), m.Q + 1, d * Number.Rational.Create(m.P)); + if (m.Q > 0) + Add(into, m.P.Add(ERational.One), m.Q - 1, d * m.Q); + return into; + }).ToList(); + // Square, one row for each monomial of the ansatz, which holds every one of the integrand's: the + // rows past it, the lowest power the derivative reaches, are zero for a solution, which an + // elimination with symbols does not decide; the derivative is checked at the end instead. + var keys = basis.Select(m => (m.P, m.Q)).ToList(); + var matrix = keys.Select(key => derivatives.Select(map => map.TryGetValue(key, out var v) ? v.InnerSimplified : Number.Integer.Zero).ToArray()).ToArray(); + // Each compound coefficient of the integrand a name while the system is solved: the elimination + // with symbols reduces a row to `0 = c` with c zero only as a value, `-3/4 d^2/(b^2 - 9/4 d^2)` + // less `d/2` times the second unknown, and declines. + var named = new List<(Variable Name, Entity Value)>(); + Entity Named(Entity value) + { + if (value is Number or Variable) + return value; + var name = Variable.CreateUnique(expr + named.Aggregate((Entity)Number.Integer.Zero, (all, pair) => all + pair.Name), "k_es"); + named.Add((name, value)); + return name; + } + var rhs = keys.Select(key => target.TryGetValue(key, out var v) ? Named(v.InnerSimplified) : Number.Integer.Zero).ToArray(); + if (!Functions.PartialFractions.TrySolveLinearWithSymbols(matrix, rhs, out var values)) + return null; + var sine = (MathS.Pow(MathS.e, u) - MathS.Pow(MathS.e, -u)) / 2; + var cosine = (MathS.Pow(MathS.e, u) + MathS.Pow(MathS.e, -u)) / 2; + Entity sum = Number.Integer.Zero; + for (var j = 0; j < basis.Count; j++) + { + var value = Functions.PartialFractions.InLowestTermsOverTheSymbols(named.Aggregate(values[j], (all, pair) => all.Substitute(pair.Name, pair.Value))); + if (value == Number.Integer.Zero) + continue; + var (p, q) = basis[j]; + Entity monomial = p.IsZero ? Number.Integer.One : MathS.Pow(sine, Number.Rational.Create(p)); + if (q == 1) + monomial *= cosine; + sum += value * monomial; + } + var answer = MathS.Pow(MathS.e, exponentOfE) * sum; + return Functions.PartialFractions.DerivativeHoldsAtSampledPoints(answer, expr, x) ? answer : null; + } + /// /// The sign of a real-valued function of among the factors of the /// integrand is constant on every interval between its zeros, and goes in front of the diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 92a33c32d..6d7dc76e0 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -1796,6 +1796,9 @@ internal static bool ASumOverRootsWouldAnswer // a sign on. Last, so that every one answered as written keeps that answer: parts answers a // polynomial beside it. if ((answer = IndefiniteIntegralSolver.SolveByNamingAnInverseSineOfASineLinearOnEachInterval(expr, x, integrateByParts)) is { }) return answer; + // And half-odd powers of a hyperbolic sine beside an exponential as the derivative of the same kind, + // last: each term alone is not elementary, and the split declined them. + if ((answer = IndefiniteIntegralSolver.SolveHalfOddPowersOfAHyperbolicSineBesideAnExponentialAsADerivative(expr, x)) is { }) return answer; return null; } /// diff --git a/Sources/Tests/UnitTests/Calculus/HalfOddPowersOfAHyperbolicSineBesideAnExponentialIntegralTest.cs b/Sources/Tests/UnitTests/Calculus/HalfOddPowersOfAHyperbolicSineBesideAnExponentialIntegralTest.cs new file mode 100644 index 000000000..ef6a173d5 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/HalfOddPowersOfAHyperbolicSineBesideAnExponentialIntegralTest.cs @@ -0,0 +1,47 @@ +// +// 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 +{ + /// + /// Half-odd powers of a hyperbolic sine beside an exponential of another linear, read as the + /// derivative of the same kind of sum: Rubi's 6.7.1 + /// e^(a + b x) sinh(c + d x)^(3/2) - 3/4 e^(a + b x) d^2/((b^2 - 9/4 d^2) sqrt(sinh(c + d x))) + /// is (e^(a + b x) (4 b sinh^(3/2) - 6 d cosh sinh^(1/2))/(4 b^2 - 9 d^2))', and each term alone + /// is not elementary, so the split declined both. + /// #718 + /// + [Trait("Area", "Calculus")] + public sealed class HalfOddPowersOfAHyperbolicSineBesideAnExponentialIntegralTest + { + [Theory] + // Rubi 6.7.1:1141. + [InlineData("e^(a + b*x)*sinh(c + d*x)^(3/2) - 3/4*e^(a + b*x)*d^2/((b^2 - 9/4*d^2)*sqrt(sinh(c + d*x)))")] + public void IsAnsweredAsADerivative(string integrand) + { + var integral = integrand.ToEntity().Integrate("x"); + var text = integral.Stringize(); + Assert.DoesNotContain("integral(", text); + Entity Pinned(Entity e) => e.Substitute("a", 0.4).Substitute("b", 1.3).Substitute("c", 0.2).Substitute("d", 0.7); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in new[] { 0.3, 0.9, 1.4, 2.1, 3.3 }) + { + 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}"); + } + } + } +}