From fdf6629e0480f33910d7a7df82bdf4e462d741f5 Mon Sep 17 00:00:00 2001 From: Rafael Vuijk Date: Mon, 28 Sep 2026 04:42:59 +0000 Subject: [PATCH] The sine and cosine integrals are integrated to A polynomial in sines and cosines of one linear, beside a polynomial, over powers of linears is integrated to Si and Ci, by their definitions: int sin(q u)/u du = Si(q u) and int cos(q u)/u du = Ci(q u). Each monomial sin^s cos^c is a sum of sines and cosines of multiples of the argument, and under u = the linear every term is u^m sin(q u) or u^m cos(q u): elementary for m >= 0, Si or Ci for m = -1, and by parts toward them below. Several linears are split into partial fractions over them first, after dividing out the polynomial part, and an argument a + b x^r beside a power of x is the same question under u = x^r. Rubi's 4.1.10, 4.1.11, 4.1.12 and 4.7.3. Part of #1501. Co-Authored-By: Claude Opus 5.5 --- BREAKING-CHANGES.md | 17 + Sources/.editorconfig | 3 + .../Integration/IndefiniteIntegralSolver.cs | 377 ++++++++++++++++++ .../Integration/Integration.Definition.cs | 3 + .../TrigonometricIntegralIntegrationTest.cs | 130 ++++++ .../Calculus/TrigonometricTowerAnsatzTest.cs | 2 +- 6 files changed, 531 insertions(+), 1 deletion(-) create mode 100644 Sources/Tests/UnitTests/Calculus/TrigonometricIntegralIntegrationTest.cs diff --git a/BREAKING-CHANGES.md b/BREAKING-CHANGES.md index 8d7aca937..88079fc99 100644 --- a/BREAKING-CHANGES.md +++ b/BREAKING-CHANGES.md @@ -1654,6 +1654,23 @@ was refused, so each was read as a product with an undeclared variable | `"Chi(1)".Simplify()` | `Chi` | `Chi(1)`, and `0.8378669409802082…` evaluated | | `"Si(x)".Differentiate("x")` | `Si` | `sin(x) / x` | +### The sine and cosine integrals are integrated to + +`sin(x)/x` was left unintegrated. A polynomial in sines and cosines of a linear, beside a +polynomial, over powers of linears is now integrated to `Si` and `Ci`: each monomial of the sines +and cosines is a sum of sines and cosines of multiples of the argument, and under `u` = the linear +each term is `u^m sin(q u)` or `u^m cos(q u)`, which is `Si(q u)` or `Ci(q u)` for `m = -1` and by +parts toward them below. Several linears are split into partial fractions over them first, and an +argument `a + b x^r` beside a power of `x` is the same question under `u = x^r` +([#1501](https://github.com/asc-community/AngouriMath/issues/1501)). + +| Input | Was (2.5.0) | Now | +|---|---|---| +| `"sin(x)/x".Integrate("x")` | `integral(sin(x) / x, x)` | `Si(x) + C` | +| `"sin(x)^2/x".Integrate("x")` | `integral(sin(x) ^ 2 / x, x)` | `1/2 * ln(x) + -1/2 * Ci(2 * x) + C` | +| `"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` | + ### 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 acc9c6cdb..238627be5 100644 --- a/Sources/.editorconfig +++ b/Sources/.editorconfig @@ -666,3 +666,6 @@ file_header_template=\nCopyright (c) 2019-2026 Angouri.\nAngouriMath is licensed [Tests/UnitTests/Convenience/TrigonometricIntegralTest.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/TrigonometricIntegralIntegrationTest.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 15f1476b7..80bce58cc 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/IndefiniteIntegralSolver.cs @@ -3888,6 +3888,383 @@ void Add(int k, bool isSine, Entity coefficient) .ToList(); } + /// + /// A polynomial in sines and cosines of one linear, times a polynomial, over a whole power + /// of a linear: P(x) T(sin(a + b x), cos(a + b x))/(e + f x)^n, onto the sine and + /// cosine integrals. Each monomial sin^s cos^c of T is a sum of sines and + /// cosines of multiples of the argument, and under u = e + f x, with + /// alpha = a - b e/f and k = b/f, sin(j (a + b x)) is + /// sin(j alpha) cos(j k u) + cos(j alpha) sin(j k u) and the cosine likewise, so every + /// term is u^m sin(q u) or u^m cos(q u): elementary for m >= 0, + /// Si(q u) and Ci(q u) for m = -1, which is their definition, and below + /// that by parts toward them. The constant term an even power leaves is a power of u + /// alone. Several linears below the bar are split into partial fractions over them, each + /// term the question for one; and an argument a + b x^r beside a power of x is + /// the same question under u = x^r, where x^p dx is u^((p + 1)/r - 1) du/r. + /// Rubi's 4.1.10, (c + d x)^m (a + b sin(e + f x))^n, 4.1.11, + /// sin(c + d x)/(x^m (a + b x)^n), and 4.1.12, (e x)^m (a + b sin(c + d x^n))^p, + /// with m negative. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + /// + /// Only where a sine or cosine integral is left: otherwise the product is elementary, and + /// the rules that write it shortest answer it. + /// + internal static Entity? SolveATrigonometricOfALinearOverAPowerOfALinear(Entity expr, Entity.Variable x) + { + // A sine or a cosine of the variable, before anything is read. + if (expr is not (Divf or Mulf or Powf or Sinf or Cosf) || ASineOrCosineArgument(expr, x) is not { } argument) + return null; + if (TreeAnalyzer.TryGetPolyLinear(argument, x, out var rate, out _) && !TreeAnalyzer.IsZero(rate)) + return OverAPowerOfALinear(expr, x, argument); + return UnderAPowerOfTheVariable(expr, x, argument); + } + + /// + /// where the argument is + /// , linear in . Two or more linears below + /// the bar are split into partial fractions over them first, each term the one-linear + /// question. + /// + private static Entity? OverAPowerOfALinear(Entity expr, Entity.Variable x, Entity argument) + { + var (sine, cosine) = (MathS.Sin(argument), MathS.Cos(argument)); + var s = Variable.CreateUnique(expr, "s"); + var c = Variable.CreateUnique(expr + s, "c"); + Entity constant = Number.Integer.One; + var linears = new List<(Entity Linear, int Power)>(); + Entity? polynomial = null; + Entity? trigonometric = null; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + { + constant = underneath ? constant / factor : constant * factor; + continue; + } + // A polynomial in the sine and the cosine, above the bar. + if (factor.ContainsNode(sine) || factor.ContainsNode(cosine)) + { + var read = factor.Replace(node => node == sine ? s : node == cosine ? c : node); + if (underneath || read.ContainsNode(x)) + return null; + trigonometric = trigonometric is null ? read : trigonometric * read; + 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; + // A linear, or a whole power of one, below the bar or as a negative power. + if (exponent < 0) + { + if (!TreeAnalyzer.TryGetPolyLinear(@base, x, out var slopeOfIt, out _) || TreeAnalyzer.IsZero(slopeOfIt)) + return null; + var at = linears.FindIndex(pair => pair.Linear == @base); + if (at < 0) + linears.Add((@base, -exponent)); + else + linears[at] = (@base, linears[at].Power - exponent); + 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 (trigonometric is null || linears.Count == 0 || linears.Count > 4 || linears.Sum(pair => pair.Power) > 12 + || MultipleAnglesOf(trigonometric, s, c) is not { } angles) + return null; + if (linears.Count == 1) + return AgainstAPowerOfALinear(polynomial, linears[0].Linear, linears[0].Power, angles, argument, x, needsAnIntegral: true) is var (one, _) && one is not null + ? Functions.PartialFractions.Bare((constant * one).InnerSimplified) + : null; + // Over each linear by partial fractions, and each term the one-linear question; the + // polynomial part first, since the split is of a proper fraction. + var denominator = linears.Aggregate((Entity)Number.Integer.One, (product, pair) => product * MathS.Pow(pair.Linear, pair.Power)); + Entity numerator = polynomial ?? Number.Integer.One; + var terms = new List(); + if (numerator.ContainsNode(x) && TreeAnalyzer.PolynomialLongDivision(numerator, denominator, genericCase: true, inTermsOf: x) is var (quotient, remainder)) + { + if (!TreeAnalyzer.IsZero(quotient)) + terms.Add(quotient); + // The remainder comes over the divisor already. + numerator = remainder is Divf(var left, _) ? left : (numerator - quotient * denominator).Expand().InnerSimplified; + } + if (!TreeAnalyzer.IsZero(numerator)) + { + if (!Functions.PartialFractions.TrySplitOverWrittenFactors(numerator, denominator, x, out var decomposition)) + return null; + // The sum of the fractions, over a constant it may come divided by. + var fractions = decomposition.InnerSimplified; + if (fractions is Divf(var over, var by) && !by.ContainsNode(x)) + fractions = Sumf.LinearChildren(over).Aggregate((Entity)Number.Integer.Zero, (all, one) => all + one / by); + terms.AddRange(Sumf.LinearChildren(fractions)); + } + Entity sum = Number.Integer.Zero; + var special = false; + foreach (var term in terms) + { + var (above, linear, power) = OverOneLinear(term, x); + var (integral, integrals) = AgainstAPowerOfALinear(above, linear, power, angles, argument, x, needsAnIntegral: false); + if (integral is null) + return null; + sum += integral; + special |= integrals; + } + return special ? Functions.PartialFractions.Bare((constant * sum).InnerSimplified) : null; + } + + /// + /// A term of a partial-fraction decomposition as what is above the bar, the linear below + /// it and its power, or the term itself above and no linear. + /// + private static (Entity Above, Entity? Linear, int Power) OverOneLinear(Entity term, Entity.Variable x) + { + Entity above = Number.Integer.One; + Entity? linear = null; + var power = 0; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(term)) + { + 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 && @base.ContainsNode(x) && (linear is null || linear == @base)) + { + linear = @base; + power -= exponent; + } + else + above = underneath ? above / factor : above * factor; + } + return (above, linear, power); + } + + /// + /// T, a polynomial in and standing for + /// the sine and the cosine of one argument, as sines and cosines of multiples of it: each + /// entry a multiple, whether a sine, and its coefficient, the multiple 0 standing for 1. + /// + private static Dictionary<(int Multiple, bool IsSine), Entity>? MultipleAnglesOf(Entity trigonometric, Entity.Variable s, Entity.Variable c) + { + if (!TreeAnalyzer.TryGetPolynomial(trigonometric, s, out var inTheSine)) + return null; + var angles = new Dictionary<(int Multiple, bool IsSine), Entity>(); + foreach (var ofOnePower in inTheSine) + { + var (sines, ofTheSines) = (ofOnePower.Key, ofOnePower.Value); + if (sines.Sign < 0 || !sines.CanFitInInt32() || !TreeAnalyzer.TryGetPolynomial(ofTheSines, c, out var inTheCosine)) + return null; + foreach (var term in inTheCosine) + { + var (cosines, coefficient) = (term.Key, term.Value); + if (cosines.Sign < 0 || !cosines.CanFitInInt32() || coefficient.ContainsNode(s) || coefficient.ContainsNode(c) + || sines.ToInt32Checked() + cosines.ToInt32Checked() > MaximumTrigonometricPower) + return null; + var (sinePower, cosinePower) = (sines.ToInt32Checked(), cosines.ToInt32Checked()); + var terms = (sinePower + cosinePower) switch + { + 0 => new List<(int Multiple, bool IsSine, Entity Coefficient)> { (0, false, Number.Integer.One) }, + 1 => new List<(int Multiple, bool IsSine, Entity Coefficient)> { (1, sinePower == 1, Number.Integer.One) }, + _ => MultipleAngles(sinePower, cosinePower), + }; + foreach (var (multiple, isSine, weight) in terms) + angles[(multiple, isSine)] = angles.TryGetValue((multiple, isSine), out var already) + ? already + coefficient * weight : coefficient * weight; + } + } + return angles; + } + + /// + /// int P(x) T(x)/L^n dx, with T given by its multiple angles and L a + /// linear or, where is , x and + /// n = 0; and whether a sine or cosine integral is left in it. Where + /// and none is, the answer is left to other rules. + /// + private static (Entity? Integral, bool Special) AgainstAPowerOfALinear(Entity? polynomial, Entity? linear, int power, + Dictionary<(int Multiple, bool IsSine), Entity> angles, Entity argument, Entity.Variable x, bool needsAnIntegral) + { + linear ??= x; + if (!TreeAnalyzer.TryGetPolyLinear(linear, x, out var f, out var e) + || !TreeAnalyzer.TryGetPolyLinear(argument, x, out var b, out var a)) + return (null, false); + var u = Variable.CreateUnique(linear + argument + (polynomial ?? Number.Integer.One), "u"); + var inU = polynomial is null ? Number.Integer.One : polynomial.Substitute(x, (u - e) / f); + if (!TreeAnalyzer.TryGetPolynomial(inU, u, out var powers) || powers.Keys.Any(degree => !degree.CanFitInInt32()) + || needsAnIntegral && powers.Keys.All(degree => degree.CompareTo(EInteger.FromInt32(power)) >= 0)) + return (null, false); + var alpha = (a - b * e / f).InnerSimplified; + var rateInU = (b / f).InnerSimplified; + Entity sum = Number.Integer.Zero; + var special = false; + foreach (var angle in angles) + { + var ((multiple, isSine), coefficient) = (angle.Key, angle.Value); + if (multiple == 0) + { + // What even powers leave: a Laurent polynomial in u. + foreach (var monomial in powers) + { + var m = monomial.Key.ToInt32Checked() - power; + sum += coefficient * monomial.Value * (m == -1 ? MathS.Ln(u) : MathS.Pow(u, m + 1) / (m + 1)); + } + continue; + } + // sin(j alpha + q u) = sin(j alpha) cos(q u) + cos(j alpha) sin(q u), and the + // cosine's; a negative q turned positive, since sin(-y) = -sin(y) and + // cos(-y) = cos(y), so that Ci is asked a positive argument where u is. + var phase = (multiple * alpha).InnerSimplified; + var (sineOfThePhase, cosineOfThePhase) = phase.Evaled is Number.Real { IsNegative: true } + ? (-MathS.Sin((-phase).InnerSimplified), MathS.Cos((-phase).InnerSimplified)) + : (MathS.Sin(phase), MathS.Cos(phase)); + var (onTheSine, onTheCosine) = isSine + ? (cosineOfThePhase, sineOfThePhase) + : (-sineOfThePhase, cosineOfThePhase); + var frequency = (multiple * rateInU).InnerSimplified; + if (frequency.Evaled is Number.Real { IsNegative: true }) + (frequency, onTheSine) = ((-frequency).InnerSimplified, -onTheSine); + Entity elementary = Number.Integer.Zero, si = Number.Integer.Zero, ci = Number.Integer.Zero; + foreach (var monomial in powers) + { + var m = monomial.Key.ToInt32Checked() - power; + var ofTheSine = APowerTimesASineOrACosine(m, sine: true, frequency, u); + var ofTheCosine = APowerTimesASineOrACosine(m, sine: false, frequency, u); + elementary += monomial.Value * (onTheSine * ofTheSine.Elementary + onTheCosine * ofTheCosine.Elementary); + si += monomial.Value * (onTheSine * ofTheSine.Si + onTheCosine * ofTheCosine.Si); + ci += monomial.Value * (onTheSine * ofTheSine.Ci + onTheCosine * ofTheCosine.Ci); + } + (si, ci) = ((coefficient * si).InnerSimplified, (coefficient * ci).InnerSimplified); + special |= !TreeAnalyzer.IsZero(si) || !TreeAnalyzer.IsZero(ci); + sum += coefficient * elementary + si * MathS.Si(frequency * u) + ci * MathS.Ci(frequency * u); + } + if (needsAnIntegral && !special) + return (null, false); + return ((sum / f).Substitute(u, linear).InnerSimplified, special); + } + + /// + /// where the argument is + /// a + b x^r, beside a power of : under u = x^r, + /// x^p T(sin(a + b x^r), cos(a + b x^r)) dx is u^((p + 1)/r - 1) T(sin(a + b u), cos(a + b u)) du/r. + /// + private static Entity? UnderAPowerOfTheVariable(Entity expr, Entity.Variable x, Entity argument) + { + // r, from x^r, b/x or b/x^k in the argument. + Entity? r = null; + foreach (var node in argument.Nodes) + { + Entity? candidate = node switch + { + Powf(var @base, var exponent) when @base == x && !exponent.ContainsNode(x) => exponent, + Divf(var above, var below) when below == x && !above.ContainsNode(x) => Number.Integer.MinusOne, + Divf(var above, Powf(var @base, var exponent)) when @base == x && !above.ContainsNode(x) && !exponent.ContainsNode(x) => -exponent, + _ => null, + }; + if (candidate is null) + continue; + if (r is not null && r != candidate) + return null; + r = candidate; + } + if (r is null || TreeAnalyzer.IsZero(r)) + return null; + var u = Variable.CreateUnique(expr, "u"); + // x^r written as u where it was read, rather than x as u^(1/r): (u^(1/n))^n is not + // folded for a symbolic n. + var argumentInU = argument.Replace(node => node switch + { + Powf(var @base, var exponent) when @base == x && exponent == r => u, + Divf(var above, var below) when below == x && r == Number.Integer.MinusOne => above * u, + Divf(var above, Powf(var @base, var exponent)) when @base == x && (-exponent) == r => above * u, + _ => node, + }).InnerSimplified; + if (argumentInU.ContainsNode(x) || !TreeAnalyzer.TryGetPolyLinear(argumentInU, u, out var slope, out _) || TreeAnalyzer.IsZero(slope)) + return null; + var (sine, cosine) = (MathS.Sin(argument), MathS.Cos(argument)); + Entity constant = Number.Integer.One; + Entity p = Number.Integer.Zero; + Entity trigonometric = Number.Integer.One; + var sawOne = false; + foreach (var (factor, underneath) in FactorsOfTheIntegrand(expr)) + { + if (!factor.ContainsNode(x)) + constant = underneath ? constant / factor : constant * factor; + else if (factor == x) + p += underneath ? -1 : 1; + else if (factor is Powf(var @base, var exponent) && @base == x && !exponent.ContainsNode(x)) + p += underneath ? -exponent : exponent; + else if (!underneath && (factor.ContainsNode(sine) || factor.ContainsNode(cosine))) + { + var inU = factor.Replace(node => node == sine ? MathS.Sin(argumentInU) : node == cosine ? MathS.Cos(argumentInU) : node); + if (inU.ContainsNode(x)) + return null; + trigonometric *= inU; + sawOne = true; + } + else + return null; + } + // Only a whole power of u leaves the shape the linear case reads. + // Simplified, since a symbolic r cancels only there: x^(-1 - r) leaves -r/r - 1, and + // bare, since the cancellation holds where r is not 0, which it is not here. + if (!sawOne || Functions.PartialFractions.Bare(((p + 1) / r - 1).Simplify()) is not Number.Integer { EInteger: var whole } + || !whole.CanFitInInt32() || whole.Sign >= 0) + return null; + var integrand = constant / r * MathS.Pow(u, whole.ToInt32Checked()) * trigonometric; + return OverAPowerOfALinear(integrand, u, argumentInU)?.Substitute(u, MathS.Pow(x, r)).InnerSimplified; + } + + /// + /// int u^m sin(q u) du, or the cosine's, as its elementary part and the coefficients + /// of Si(q u) and Ci(q u) in it: Si(q u) and Ci(q u) themselves for + /// m = -1; below that + /// int u^m sin(q u) = u^(m + 1) sin(q u)/(m + 1) - q/(m + 1) int u^(m + 1) cos(q u) and + /// int u^m cos(q u) = u^(m + 1) cos(q u)/(m + 1) + q/(m + 1) int u^(m + 1) sin(q u); and + /// above it -u^m cos(q u)/q + m/q int u^(m - 1) cos(q u) and + /// u^m sin(q u)/q - m/q int u^(m - 1) sin(q u), all by parts. + /// + private static (Entity Elementary, Entity Si, Entity Ci) APowerTimesASineOrACosine(int m, bool sine, Entity q, Entity.Variable u) + { + if (m == -1) + return sine + ? (Number.Integer.Zero, Number.Integer.One, Number.Integer.Zero) + : (Number.Integer.Zero, Number.Integer.Zero, Number.Integer.One); + if (m < -1) + { + var (elementary, si, ci) = APowerTimesASineOrACosine(m + 1, !sine, q, u); + var by = (sine ? -q / (m + 1) : q / (m + 1)).InnerSimplified; + // u^(m + 1)/(m + 1), written as -1/(|m + 1| u^|m + 1|) + var written = -(sine ? MathS.Sin(q * u) : MathS.Cos(q * u)) / (-(m + 1) * MathS.Pow(u, -(m + 1))); + return ((written + by * elementary).InnerSimplified, (by * si).InnerSimplified, (by * ci).InnerSimplified); + } + var first = sine ? -MathS.Pow(u, m) * MathS.Cos(q * u) / q : MathS.Pow(u, m) * MathS.Sin(q * u) / q; + if (m == 0) + return (first.InnerSimplified, Number.Integer.Zero, Number.Integer.Zero); + var (below, _, _) = APowerTimesASineOrACosine(m - 1, !sine, q, u); + return ((first + (sine ? m / q : -m / q) * below).InnerSimplified, Number.Integer.Zero, Number.Integer.Zero); + } + + /// + /// The argument of a sine or a cosine of among the factors of + /// , above the bar and read through whole powers and sums, or + /// . + /// + private static Entity? ASineOrCosineArgument(Entity expr, Entity.Variable x) => expr switch + { + Sinf(var argument) when argument.ContainsNode(x) => argument, + Cosf(var argument) when argument.ContainsNode(x) => argument, + Powf(var @base, Number.Integer) => ASineOrCosineArgument(@base, x), + Mulf(var left, var right) => ASineOrCosineArgument(left, x) ?? ASineOrCosineArgument(right, x), + Divf(var numerator, _) => ASineOrCosineArgument(numerator, x), + Sumf(var left, var right) => ASineOrCosineArgument(left, x) ?? ASineOrCosineArgument(right, x), + Minusf(var left, var right) => ASineOrCosineArgument(left, x) ?? ASineOrCosineArgument(right, x), + _ => null, + }; + /// /// int P(x) e^(a x) (c cos(b x) + d sin(b x)) dx, closed by as many rounds of parts /// as P has degree. diff --git a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs index 8c07bd052..9fd928433 100644 --- a/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs +++ b/Sources/AngouriMath/Functions/Continuous/Integration/Integration.Definition.cs @@ -594,6 +594,9 @@ 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; + // 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; // 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. diff --git a/Sources/Tests/UnitTests/Calculus/TrigonometricIntegralIntegrationTest.cs b/Sources/Tests/UnitTests/Calculus/TrigonometricIntegralIntegrationTest.cs new file mode 100644 index 000000000..ae1db9b47 --- /dev/null +++ b/Sources/Tests/UnitTests/Calculus/TrigonometricIntegralIntegrationTest.cs @@ -0,0 +1,130 @@ +// +// 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 sine and cosine integrals: sines and cosines of a linear + /// over a power of a linear. + /// https://github.com/asc-community/AngouriMath/issues/1501 + /// + [Trait("Area", "Calculus")] + public sealed class TrigonometricIntegralIntegrationTest + { + 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"), ("k", "2") }; + + /// + /// Under u = c + d x, each term is u^m sin(q u) or u^m cos(q u), and + /// m = -1 is Si(q u) or Ci(q u), by parts below that. Rubi's 4.1.10, + /// (c + d x)^m sin(a + b x) with m negative. + /// + [Theory] + [InlineData("sin(x)/x")] + [InlineData("cos(x)/x")] + [InlineData("sin(2*x + 1)/(3*x + 4)")] + [InlineData("cos(3*x - 1)/(1 - x)")] + [InlineData("sin(x)/x^2")] + [InlineData("cos(x)/x^3")] + [InlineData("x*sin(x)/(x + 1)^2")] + [InlineData("(1 + x^2)*cos(x)/x^2")] + [InlineData("x^(-2)*sin(2*x)")] + [InlineData("sin(-2*x)/x")] + [InlineData("sin(a + b*x)/(c + d*x)")] + [InlineData("cos(a + b*x)/(c + k*x)^2")] + public void ASineOrACosineOfALinearOverAPowerOfALinear(string integrand) + => DifferentiatesBack(integrand, AroundZero, Pins); + + /// + /// A power of a sine or a cosine, or a product of both, is a sum of sines and cosines of + /// multiples of the argument, and an even power leaves a constant term, whose integral is + /// a logarithm: int sin(x)^2/x = ln(x)/2 - Ci(2 x)/2. + /// + [Theory] + [InlineData("sin(x)^2/x")] + [InlineData("cos(x)^2/x")] + [InlineData("sin(x)*cos(x)/x")] + [InlineData("sin(x)^3/x^2")] + [InlineData("sin(2*x + 1)^2*cos(2*x + 1)/(x + 1)")] + public void APowerOfASineOrACosineOverAPowerOfALinear(string integrand) + => DifferentiatesBack(integrand, AroundZero, Pins); + + /// + /// Two or more linears below the bar are split into partial fractions over them, and + /// each term is the one-linear question. Rubi's 4.1.11, sin(c + d x)/(x^m (a + b x)^n). + /// + [Theory] + [InlineData("sin(x)/(x*(x + 1))")] + [InlineData("cos(2*x + 1)/((x - 1)*(x + 2)^2)")] + [InlineData("sin(c + d*x)/(x^2*(a + b*x))")] + [InlineData("x^3*sin(x)^2/((x + 1)*(x + 3))")] + public void OverSeveralLinears(string integrand) + => DifferentiatesBack(integrand, AroundZero, Pins); + + /// + /// An argument a + b x^r beside a power of x is the same question under + /// u = x^r, for a whole, a fractional, a negative and a symbolic r. Rubi's + /// 4.1.12, (e x)^m (a + b sin(c + d x^n))^p. + /// + [Theory] + [InlineData("sin(x^2)/x")] + [InlineData("(a + b*sin(c + d*x^2))^2/x^3")] + [InlineData("sin(sqrt(x))/x")] + [InlineData("sin(a + b/x)")] + [InlineData("x^2*sin(a + b/x)^2")] + [InlineData("x^(-1 - k)*sin(a + b*x^k)")] + public void UnderAPowerOfTheVariable(string integrand) + => DifferentiatesBack(integrand, new[] { 0.35, 0.9, 1.45, 2.3 }, Pins); + + [Theory] + [InlineData("sin(x)/x", "Si(x)")] + [InlineData("cos(x)/x", "Ci(x)")] + [InlineData("sin(2*x)/x", "Si(2 * x)")] + [InlineData("sin(x)*cos(x)/x", "1/2 * Si(2 * x)")] + public void TheDefinitionsAreTheAnswers(string integrand, string expected) + => Assert.Equal(expected.ToEntity() + "C".ToEntity(), integrand.ToEntity().Integrate("x")); + + /// + /// Where the cosine integrals of two terms cancel, the antiderivative is elementary, and + /// simplified it is written without them: (x cos(x) - sin(x))/x^2 is the derivative + /// of sin(x)/x. + /// + [Fact] + public void AnElementaryAntiderivativeKeepsNoSineOrCosineIntegral() + { + var integral = DifferentiatesBack("(x*cos(x) - sin(x))/x^2", AroundZero).Simplify().Stringize(); + Assert.DoesNotContain("Si(", integral); + Assert.DoesNotContain("Ci(", integral); + } + } +} diff --git a/Sources/Tests/UnitTests/Calculus/TrigonometricTowerAnsatzTest.cs b/Sources/Tests/UnitTests/Calculus/TrigonometricTowerAnsatzTest.cs index 7ac4fcb34..a8571bc71 100644 --- a/Sources/Tests/UnitTests/Calculus/TrigonometricTowerAnsatzTest.cs +++ b/Sources/Tests/UnitTests/Calculus/TrigonometricTowerAnsatzTest.cs @@ -93,7 +93,7 @@ public void TheAnswerIsInTheSineAndCosine() /// no P for any Q, and nothing else reads the shape. /// [Theory] - [InlineData("sin(x)/x")] + [InlineData("tan(x)/x")] [InlineData("x/sin(x)")] [InlineData("x*tan(x)")] [InlineData("x*sin(x)/(1 + x^2)")]