diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index ffaca613a..9dbee9310 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1703,6 +1703,21 @@ argument `a + b x^r` beside a power of `x` is the same question under `u = x^r` | `"sin(x)/(x*(1+x))".Integrate("x")` | `integral(sin(x) / (x * (1 + x)), x)` | `Si(x) + -cos(1) * Si(1 + x) + sin(1) * Ci(1 + x) + C` | | `"sin(a+b/x)/x".Integrate("x")` | `integral(sin(a + b / x) / x, x)` | `-(cos(a) * Si(b * 1 / x) + sin(a) * Ci(b * 1 / x)) + C` | +### The hyperbolic sine and cosine integrals are integrated to + +`sinh(x)/x` was left unintegrated. `sinh` and `cosh` arrive as exponentials, and a polynomial in +exponentials of a linear, beside a polynomial, over a whole power of a linear is now integrated to +`Shi` and `Chi`. Under `u` = the linear each term is `u^m e^(k u)`, the exponential integral's, and +two terms of opposite rates pair: `A Ei(k u) + B Ei(-k u) = (A + B) Chi(k u) + (A - B) Shi(k u)`, up +to a constant ([#1501](https://github.com/asc-community/AngouriMath/issues/1501)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sinh(x)/x".Integrate("x")` | `integral((e ^ x - e ^ (-x)) / 2 / x, x)` | `Shi(x) + C` | +| `"cosh(x)/x".Integrate("x")` | `integral((e ^ x + e ^ (-x)) / 2 / x, x)` | `Chi(x) + C` | +| `"sinh(x)^2/x".Integrate("x")` | `integral(((e ^ x - e ^ (-x)) / 2) ^ 2 / x, x)` | `-1/2 * ln(x) + 1/2 * Chi(2 * x) + C` | +| `"cosh(a+b*x)/(c+d*x)".Integrate("x")` | `integral((e ^ (a + b * x) + e ^ (-(a + b * x))) / 2 / (c + d * x), x)` | `Chi` and `Shi` of `b (c + d x)/d`, with `cosh(a - b c/d)/d` and `sinh(a - b c/d)/d` in front, written in exponentials | + ### A definite integral is the limit at a bound where its antiderivative is undefined `integral(x * e^(-x), x, 0, +oo)` was `NaN`, this library's way of saying a value does not diff --git a/Sources/.editorconfig b/Sources/.editorconfig index 435b4248b..80e2bd225 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -672,3 +672,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Core/DowncastingOffComparisonTest.cs] file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n + +[Tests/UnitTests/Calculus/HyperbolicIntegralIntegrationTest.cs] +file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed under MIT.\nDetails: https://github.com/asc-community/AngouriMath/blob/master/LICENSE.md.\nWebsite: https://am.angouri.org.\n diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs index 992904f4d..999bf3bb6 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -3350,6 +3350,209 @@ private static (Entity Elementary, Entity EiCoefficient) APowerTimesAnExponentia return (exponential * sum, Number.Integer.Zero); } + /// + /// A sum of exponentials of the variable, times a polynomial, over a whole power of a + /// linear: P(x) sum_j c_j e^(r_j x)/(e + f x)^n, which is how a polynomial in + /// sinh and cosh of a linear arrives, onto the hyperbolic sine and cosine + /// integrals. Under u = e + f x each term is u^m e^(k u), the exponential + /// integral's (), and two of opposite rates pair: + /// A Ei(k u) + B Ei(-k u) = (A + B) Chi(k u) + (A - B) Shi(k u), since + /// Ei(y) + Ei(-y) = 2 Chi(y) and Ei(y) - Ei(-y) = 2 Shi(y) for y > 0, and + /// for y < 0 up to a constant. A rate of 0, which an even power leaves, is a power + /// of u alone. Rubi's 6.1.1, (c + d x)^m (a + b sinh(e + f x))^n, with + /// m negative. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + /// + /// After , which answers one + /// exponential alone; and only where an integral is left. + /// + internal static Entity? SolveAHyperbolicOfALinearOverAPowerOfALinear(Entity expr, Entity.Variable x) + { + if (expr is not (Divf or Mulf) || !HasASumOfExponentials(expr, x)) + return null; + Entity constant = Number.Integer.One; + Entity? linear = null; + var power = 0; + Entity? polynomial = null; + Entity? exponentials = null; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + var (@base, exponent) = factor is Powf(var raised, Number.Integer whole) && whole.EInteger.CanFitInInt32() && !whole.EInteger.IsZero + ? (raised, whole.EInteger.ToInt32Checked()) : (factor, 1); + if (underneath) + exponent = -exponent; + if (exponent < 0 && TreeAnalyzer.TryGetPolyLinear(@base, x, out var slopeOfIt, out _) && !TreeAnalyzer.IsZero(slopeOfIt)) + { + if (linear is not null && linear != @base) + return null; + linear = @base; + power -= exponent; + continue; + } + if (underneath) + return null; + if (factor.Nodes.Any(node => node is Powf(var b, var e) && !b.ContainsNode(x) && e.ContainsNode(x))) + { + exponentials = exponentials is null ? factor : exponentials * factor; + continue; + } + if (TreeAnalyzer.TryGetPolynomial(factor, x, out var monomials) && monomials.Keys.All(degree => degree.Sign >= 0 && degree.CompareTo(EInteger.FromInt32(12)) <= 0)) + { + polynomial = polynomial is null ? factor : polynomial * factor; + continue; + } + return null; + } + if (exponentials is null || linear is null || power > 12 || !TreeAnalyzer.TryGetPolyLinear(linear, x, out var f, out var e0)) + return null; + // The exponentials as a sum of c e^(r x), gathered by the rate r. + var rates = new List<(Entity Rate, Entity Coefficient)>(); + foreach (var term in Sumf.LinearChildren(exponentials.Expand())) + { + if (ReadAnExponentialTerm(term, x) is not { } read) + return null; + // The same rate may be written two ways, `2 (-b)` and `-2 b`. + var at = rates.FindIndex(pair => pair.Rate == read.Rate + || pair.Rate.Vars.Any() && TreeAnalyzer.IsZero(Functions.PartialFractions.Bare((pair.Rate - read.Rate).Simplify()))); + if (at < 0) + rates.Add(read); + else + rates[at] = (read.Rate, rates[at].Coefficient + read.Coefficient); + } + var u = Variable.CreateUnique(expr, "u"); + var inU = polynomial is null ? Number.Integer.One : polynomial.Substitute(x, (u - e0) / f); + if (!TreeAnalyzer.TryGetPolynomial(inU, u, out var powers) || powers.Keys.All(degree => degree.CompareTo(EInteger.FromInt32(power)) >= 0) + || powers.Keys.Any(degree => !degree.CanFitInInt32())) + return null; + Entity sum = Number.Integer.Zero; + // The constant in front goes into each coefficient, so that the numbers fold. + var overall = (constant / f).InnerSimplified; + // The coefficient of Ei(k u) for each rate k in u, to be paired. + var integrals = new List<(Entity Rate, Entity Coefficient)>(); + foreach (var (rate, coefficientOfTheRate) in rates) + { + var coefficient = (overall * coefficientOfTheRate).InnerSimplified; + if (TreeAnalyzer.IsZero(rate)) + { + // The terms of rate 0 come from several products and may cancel, as in + // `cosh(y) sinh(y)^3`, where `e^(2 a) e^(-2 a)` has to be seen to be 1. + var constantPart = Functions.PartialFractions.Bare(coefficient.Simplify()); + if (TreeAnalyzer.IsZero(constantPart)) + continue; + foreach (var monomial in powers) + { + var m = monomial.Key.ToInt32Checked() - power; + sum += constantPart * monomial.Value * (m == -1 ? MathS.Ln(u) : MathS.Pow(u, m + 1) / (m + 1)); + } + continue; + } + // e^(r x) = e^(-r e/f) e^(k u) with k = r/f. + var k = (rate / f).InnerSimplified; + var shifted = (coefficient * MathS.Pow(MathS.e, (-rate * e0 / f).InnerSimplified)).InnerSimplified; + var exponentialOfU = MathS.Pow(MathS.e, k * u); + Entity elementary = Number.Integer.Zero, ei = Number.Integer.Zero; + foreach (var monomial in powers) + { + var (below, eiBelow) = APowerTimesAnExponential(monomial.Key.ToInt32Checked() - power, k, exponentialOfU, u); + elementary += monomial.Value * below; + ei += monomial.Value * eiBelow; + } + sum += shifted * elementary; + ei = (shifted * ei).InnerSimplified; + if (!TreeAnalyzer.IsZero(ei)) + integrals.Add((k, ei)); + } + if (integrals.Count == 0) + return null; + // Pairs of opposite rates as Chi and Shi, the positive rate written. + while (integrals.Count > 0) + { + var (k, a) = integrals[0]; + integrals.RemoveAt(0); + // Simplified, since `d/f` and `-d/f` do not cancel inner-simplified. + var opposite = integrals.FindIndex(pair => TreeAnalyzer.IsZero(Functions.PartialFractions.Bare((pair.Rate + k).Simplify()))); + if (opposite < 0) + { + sum += Functions.PartialFractions.Bare(a.Simplify()) * MathS.Ei(k * u); + continue; + } + var b = integrals[opposite].Coefficient; + integrals.RemoveAt(opposite); + if (k.Evaled is Number.Real { IsNegative: true }) + (k, a, b) = ((-k).InnerSimplified, b, a); + // The two coefficients simplified, which is where `e e^(-4/3)` becomes `e^(-1/3)`. + sum += Functions.PartialFractions.Bare((a + b).Simplify()) * MathS.Chi(k * u) + + Functions.PartialFractions.Bare((a - b).Simplify()) * MathS.Shi(k * u); + } + return Functions.PartialFractions.Bare(sum.Substitute(u, linear).InnerSimplified); + } + + /// + /// A term c e^(r x) of an expanded sum of exponentials, as its rate and coefficient: + /// a product of constants and of constants to linears in , or + /// . + /// + private static (Entity Rate, Entity Coefficient)? ReadAnExponentialTerm(Entity term, Entity.Variable x) + { + Entity rate = Number.Integer.Zero; + Entity coefficient = Number.Integer.One; + foreach (var (written, underneath) in FactorsOfTheIntegrand(term)) + if (!Read(written, underneath ? -1 : 1)) + return null; + // The rate is a sum of the slopes, `b + b - b - b` for four exponentials of `a + b x`, + // and inner simplification leaves that as `2 b + 2 (-b)`, a zero nothing reads: a + // rate with a symbol in it is simplified. + var folded = rate.InnerSimplified; + if (folded.Vars.Any()) + folded = Functions.PartialFractions.Bare(folded.Simplify()); + return (folded, coefficient.InnerSimplified); + + // A factor to a whole power: a constant, an exponential of a linear, or a product or + // a quotient of those -- `(e^(2 x)/4)^2`, as the expansion of `sinh(x)^2 cosh(x)^2` + // leaves one, is `e^(4 x)/16`. + bool Read(Entity factor, int times) + { + while (factor is Powf(var inner, Number.Integer whole) && whole.EInteger.CanFitInInt32()) + (factor, times) = (inner, times * whole.EInteger.ToInt32Checked()); + if (!factor.ContainsNode(x)) + { + coefficient *= times == 1 ? factor : MathS.Pow(factor, times); + return true; + } + switch (factor) + { + case Mulf(var left, var right): + return Read(left, times) && Read(right, times); + case Divf(var numerator, var denominator): + return Read(numerator, times) && Read(denominator, -times); + case Powf(var @base, var exponent) when !@base.ContainsNode(x) + && TreeAnalyzer.TryGetPolyLinear(exponent, x, out var slope, out var intercept): + var logarithm = @base == MathS.e ? (Entity)Number.Integer.One : MathS.Ln(@base); + rate += times * slope * logarithm; + coefficient *= MathS.Pow(@base, times * intercept); + return true; + default: + return false; + } + } + } + + /// Whether a factor of is a sum with an exponential of the variable in it. + private static bool HasASumOfExponentials(Entity expr, Entity.Variable x) => expr switch + { + Sumf or Minusf => expr.Nodes.Any(node => node is Powf(var b, var e) && !b.ContainsNode(x) && e.ContainsNode(x)), + Powf(var @base, Number.Integer) => HasASumOfExponentials(@base, x), + Mulf(var left, var right) => HasASumOfExponentials(left, x) || HasASumOfExponentials(right, x), + Divf(var numerator, _) => HasASumOfExponentials(numerator, x), + _ => false, + }; + /// /// A factor of , read through products and quotients, that is a /// constant to a linear in , above the bar, or . diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 9fd928433..171b8ff6f 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -594,6 +594,8 @@ private static Entity Normalized(Entity expr, Entity.Variable x) => if ((answer = IndefiniteIntegralSolver.SolveAProductOfExponentialsOfQuadratics(expr, x)) is { }) return answer; // An exponential of a linear over a power of a linear, onto the exponential integral. if ((answer = IndefiniteIntegralSolver.SolveAnExponentialOfALinearOverAPowerOfALinear(expr, x)) is { }) return answer; + // And a sum of exponentials, which is how sinh and cosh arrive, onto Shi and Chi. + if ((answer = IndefiniteIntegralSolver.SolveAHyperbolicOfALinearOverAPowerOfALinear(expr, x)) is { }) return answer; // Sines and cosines of a linear over a power of a linear, onto Si and Ci under u = the // linear, where the search would take the quotient by parts without end. if ((answer = IndefiniteIntegralSolver.SolveATrigonometricOfALinearOverAPowerOfALinear(expr, x)) is { }) return answer; diff --git a/Sources/Tests/UnitTests/Calculus/HyperbolicIntegralIntegrationTest.cs b/Sources/Tests/UnitTests/Calculus/HyperbolicIntegralIntegrationTest.cs new file mode 100644 index 000000000..c95628726 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/HyperbolicIntegralIntegrationTest.cs @@ -0,0 +1,77 @@ +// +// 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 System.Linq; +using AngouriMath.Extensions; +using Xunit; + +namespace AngouriMath.Tests.Calculus +{ + /// + /// The integrals that answer in the hyperbolic sine and cosine integrals: sums of + /// exponentials of a linear, which is how sinh and cosh arrive, over a power of + /// a linear. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + [Trait("Area", "Calculus")] + public sealed class HyperbolicIntegralIntegrationTest + { + private static Entity DifferentiatesBack(string integrand, double[] points, params (string Name, string Value)[] pins) + { + var integral = integrand.ToEntity().Integrate("x"); + Assert.DoesNotContain("integral(", integral.Stringize()); + Entity Pinned(Entity e) => pins.Aggregate(e, (current, pin) => current.Substitute(pin.Name, pin.Value.ToEntity())); + var derivative = Pinned(integral.Substitute("C", 0)).Differentiate("x"); + var original = Pinned(integrand.ToEntity()); + foreach (var at in points) + { + var got = derivative.Substitute("x", at).EvalNumerical(); + var want = original.Substitute("x", at).EvalNumerical(); + 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)); + Assert.True(difference / scale < 1e-9, $"d/dx of {integral} is {got} at x = {at}, where {integrand} is {want}"); + } + return integral; + } + + /// Either side of 0 and of every pole below. + private static readonly double[] AroundZero = { -1.7, 0.35, 0.9, 1.45 }; + + private static readonly (string, string)[] Pins = + { ("a", "1/3"), ("b", "2/3"), ("c", "5/4"), ("d", "1/2") }; + + /// + /// Under u = c + d x, each exponential is u^m e^(k u), and the two of + /// opposite rates pair: A Ei(k u) + B Ei(-k u) = (A + B) Chi(k u) + (A - B) Shi(k u). + /// Rubi's 6.1.1 and 6.2.1, (c + d x)^m sinh(a + b x)^n with m negative, and + /// 6.7.1's products of powers, whose expansion writes one rate several ways -- and a zero + /// rate as 2 b + 2 (-b). + /// + [Theory] + [InlineData("sinh(x)/x")] + [InlineData("cosh(x)/x")] + [InlineData("sinh(2*x + 1)/(3*x + 4)")] + [InlineData("cosh(x)/x^2")] + [InlineData("x*sinh(x)/(x + 1)^2")] + [InlineData("sinh(x)^2/x")] + [InlineData("sinh(x)*cosh(x)/x")] + [InlineData("cosh(a + b*x)/(c + d*x)")] + [InlineData("cosh(a + b*x)*sinh(a + b*x)^3/x")] + [InlineData("sinh(a + b*x)^2*cosh(a + b*x)^2/(c + d*x)")] + [InlineData("(e^x + 2*e^(3*x))/x")] + public void ASumOfExponentialsOverAPowerOfALinear(string integrand) + => DifferentiatesBack(integrand, AroundZero, Pins); + + [Theory] + [InlineData("sinh(x)/x", "Shi(x)")] + [InlineData("cosh(x)/x", "Chi(x)")] + [InlineData("sinh(2*x)/x", "Shi(2 * x)")] + public void TheDefinitionsAreTheAnswers(string integrand, string expected) + => Assert.Equal(expected.ToEntity() + "C".ToEntity(), integrand.ToEntity().Integrate("x")); + } +}