# Mathematica 13.0 simplifies trigonometric integral wrong

I calculate the following general expression:

FullSimplify[Integrate[2/L*Sin[(2*\[Pi]*q*x)/L+\[Phi]1]*Sin[(2*\[Pi]*q*x)/L+\[Phi]2]*Sin[(2*\[Pi]*qt*x)/L+\[Phi]3],{x,0,L}],Element[q|qt,Integers]]


and I get:

0


This is wrong in general because of the following special case. For qt = 2q:

FullSimplify[Integrate[2/L*Sin[(2*\[Pi]*q*x)/L+\[Phi]1]*Sin[(2*\[Pi]*q*x)/L+\[Phi]2]*Sin[(2*\[Pi]*qt*x)/L+\[Phi]3]/.{qt->2*q},{x,0,L}],Element[q|qt,Integers]]


I get the sine of some phases:

1/2 Sin[\[Phi]1+\[Phi]2-\[Phi]3]


I have confirmed the special case with pen and paper. Can somebody please tell me why Mathematica screws the general case up like this? It is a plain old mistake by the most sophisticated symbolic maths engine on the planet. Or...?

It's not too uncommon for Integrate to return something which has a special case that can't be evaluated properly without taking a limit. In your case there is a denominator equal to 0 when qt == 2q, but you get the correct formula from Limit[yourIntegral, qt -> 2 q].

When FullSimplify sees the denominator it concludes qt != 2q which makes the expression 0.

• Why not return a conditional such as "0 if qt != 2q"? Then the user could check the special case themselves. That would make mathematica much more useful for mathematicians... (I am actually used to correct conditionals returned by Mathematica. I don't expect Mathematica to make mistakes like this.) May 22, 2022 at 15:58
• That is great, but ClearAll[q,qt,L,\[Phi]1 ,\[Phi]2 ,\[Phi]3 ];Integrate[ 2/L*Sin[(2*\[Pi]*q*x)/L + \[Phi]1]*Sin[(2*\[Pi]*q*x)/L + \[Phi]2]* Sin[(2*\[Pi]*qt*x)/L + \[Phi]3], {x, 0, L}, Assumptions -> {q, qt} \[Element] Integers && L > 0 && \[Phi]1 > 0 && \[Phi]3 > 0 && \[Phi]2 > 0, GenerateConditions -> True] on a fresh kernel returns a long expression. May 22, 2022 at 16:24
• @dimitsev FunctionSingularities[res, {q, qt}] returns the conditions you seek. The issue is, as Coolwater pointed out, is the simplification functions make transformations that are generically true. FWIW, this handles one of the singularities: res /. Sin[z_] :> z*Sinc[z] // Simplify. May 22, 2022 at 16:38

Here's a way to eliminate the singularities in the answer, provided you can be happy removing a singularity of the form Sin[u]/u with Sinc[u]. We end up with a very nice simplification.

There are three singularities, which may be found by inspection or by computing

FunctionSingularities[res, {q, qt, ϕ1, ϕ2, ϕ3}]
(*  qt == 0 || 4 q^2 - qt^2 == 0  *)


The singularities at u == 2 q - qt and u == 2 u + qt are more easily dealt with via substitutions:

qt ->  2 q + u
qt -> -2 q + u


The basic idea is to expand the result into a sum of terms that have 1/u and terms that do not. This is possible in the OP's result but obviously not in general. You have to check this. Since in the OP's case the singularity at 1/u is removable, we apply the sinc identity on the terms with 1/u. The case where u = qt has a particularly simple solution that does not require separating the terms. The other two use Cases[] to separate the terms.

res = FullSimplify@ Integrate[
2/L*Sin[(2*π*q*x)/L + ϕ1]*Sin[(2*π*q*x)/L + ϕ2]*
Sin[(2*π*qt*x)/L + ϕ3], {x, 0, L}];
singqt = res /. Sin[z : Pi qt] :> z*Sinc[z] // Simplify;
singminus = singqt /. qt -> 2 q + u // Simplify // TrigExpand;
singminus = Simplify[
Total@ Cases[singminus, t_ /; ! FreeQ[t, 1/u]] /.
{1/u -> Pi*Sinc[Pi u]/Sin[Pi u]}
] + Total@ Cases[singminus, t_ /; FreeQ[t, 1/u]] /. u -> qt - 2 q;
singplus = singminus /. qt -> -2 q + u // Simplify // TrigExpand;
singplus = Simplify[
Total@ Cases[singplus, t_ /; ! FreeQ[t, 1/u]] /.
{1/u -> Pi*Sinc[Pi u]/Sin[Pi u]}
] + Total@ Cases[singplus, t_ /; FreeQ[t, 1/u]] /. u -> qt + 2 q;
integral = singplus // FullSimplify
(*
1/2 (2 Cos[ϕ1 - ϕ2] Sin[π qt + ϕ3] Sinc[π qt] +
Sin[2 π q - π qt + ϕ1 + ϕ2 - ϕ3] Sinc[π (-2 q + qt)] -
Sin[2 π q + π qt + ϕ1 + ϕ2 + ϕ3] Sinc[π (2 q + qt)])
*)


Check:

res - integral // FunctionExpand // FullSimplify
(*  0  *)


Since integral is continuous, it makes sense that there are no conditions on the result res from an integration point of view. That res has removable singularities may seem unfortunate to some, but others sometimes wish Mathematica would stick to elementary functions when it can.

Update

I should have taken it one more step:

FullSimplify[integral, {q, qt} ∈ Integers]
(*
((-1)^qt *
(KroneckerDelta[qt] 2 Cos[ϕ1 - ϕ2] Sin[ϕ3] +
KroneckerDelta[2*q - qt] Sin[ϕ1 + ϕ2 - ϕ3] -
KroneckerDelta[2 q + qt] Sin[ϕ1 + ϕ2 + ϕ3]))/2
*)


Update 2

Here's a simplified way to get the simplified integral, which can then be simplified further to get KroneckerDelta[] as above:

integral = res // Apart[#, q] & //
Map[# /. Sin[Pi z_] :> Pi z Sinc[Pi z] &]


Under natural assumptions one obtains

Integrate[2/L*Sin[(2*\[Pi]*q*x)/L + \[Phi]1]*
Sin[(2*\[Pi]*q*x)/L + \[Phi]2]*
Sin[(2*\[Pi]*qt*x)/L + \[Phi]3], {x, 0, L},
Assumptions -> {q, qt} \[Element] Integers &&
L > 0 && \[Phi]1 > 0 && \[Phi]3 > 0 && \[Phi]2 > 0,
GenerateConditions -> True]


((-4 q^2 + qt^2) Cos[\[Phi]1 - \[Phi]2 - \[Phi]3] - qt (2 q + qt) Cos[\[Phi]1 + \[Phi]2 - \[Phi]3] + 2 q qt Cos[4 \[Pi] q - 2 \[Pi] qt + \[Phi]1 + \[Phi]2 - \[Phi]3] + qt^2 Cos[4 \[Pi] q - 2 \[Pi] qt + \[Phi]1 + \[Phi]2 - \[Phi]3] - 4 q^2 Cos[\[Phi]1 - \[Phi]2 + \[Phi]3] + qt^2 Cos[\[Phi]1 - \[Phi]2 + \[Phi]3] + 4 q^2 Cos[2 \[Pi] qt + \[Phi]1 - \[Phi]2 + \[Phi]3] - qt^2 Cos[2 \[Pi] qt + \[Phi]1 - \[Phi]2 + \[Phi]3] + 4 q^2 Cos[2 \[Pi] qt - \[Phi]1 + \[Phi]2 + \[Phi]3] - qt^2 Cos[2 \[Pi] qt - \[Phi]1 + \[Phi]2 + \[Phi]3] + 2 q qt Cos[\[Phi]1 + \[Phi]2 + \[Phi]3] - qt^2 Cos[\[Phi]1 + \[Phi]2 + \[Phi]3] - 2 q qt Cos[4 \[Pi] q + 2 \[Pi] qt + \[Phi]1 + \[Phi]2 + \[Phi]3] + qt^2 Cos[ 4 \[Pi] q + 2 \[Pi] qt + \[Phi]1 + \[Phi]2 + \[Phi]3])/(4 \[Pi] qt (-4 q^2 + qt^2))

Limit[%, qt -> 0]


((4 \[Pi] q Cos[\[Phi]1 - \[Phi]2] + Sin[\[Phi]1 + \[Phi]2] - Sin[4 \[Pi] q + \[Phi]1 + \[Phi]2]) Sin[\[Phi]3])/(4 \[Pi] q)

Limit[%%, qt -> 2 q]


(1/(16 \[Pi] q))(2 Cos[\[Phi]1 - \[Phi]2 - \[Phi]3] + 2 Cos[\[Phi]1 - \[Phi]2 + \[Phi]3] - 2 Cos[4 \[Pi] q + \[Phi]1 - \[Phi]2 + \[Phi]3] - 2 Cos[4 \[Pi] q - \[Phi]1 + \[Phi]2 + \[Phi]3] - Cos[\[Phi]1 + \[Phi]2 + \[Phi]3] + Cos[8 \[Pi] q + \[Phi]1 + \[Phi]2 + \[Phi]3] + 8 \[Pi] q Sin[\[Phi]1 + \[Phi]2 - \[Phi]3])

Limit[%%%, qt -> -2 q]


-(1/(16 \[Pi] q))(2 Cos[\[Phi]1 - \[Phi]2 - \[Phi]3] - 2 Cos[4 \[Pi] q + \[Phi]1 - \[Phi]2 - \[Phi]3] - 2 Cos[4 \[Pi] q - \[Phi]1 + \[Phi]2 - \[Phi]3] - Cos[\[Phi]1 + \[Phi]2 - \[Phi]3] + Cos[8 \[Pi] q + \[Phi]1 + \[Phi]2 - \[Phi]3] + 2 Cos[\[Phi]1 - \[Phi]2 + \[Phi]3] + 8 \[Pi] q Sin[\[Phi]1 + \[Phi]2 + \[Phi]3])

FullSimplify[-(1/(16 \[Pi] q)) (2 Cos[\[Phi]1 - \[Phi]2 - \[Phi]3] -
2 Cos[4 \[Pi] q + \[Phi]1 - \[Phi]2 - \[Phi]3] -
2 Cos[4 \[Pi] q - \[Phi]1 + \[Phi]2 - \[Phi]3] -
Cos[\[Phi]1 + \[Phi]2 - \[Phi]3] +
Cos[8 \[Pi] q + \[Phi]1 + \[Phi]2 - \[Phi]3] +
2 Cos[\[Phi]1 - \[Phi]2 + \[Phi]3] +
8 \[Pi] q Sin[\[Phi]1 + \[Phi]2 + \[Phi]3]),
Assumptions -> {q, qt} \[Element] Integers]


-(1/2) Sin[\[Phi]1 + \[Phi]2 + \[Phi]3]

• The same is obtained only with Integrate[ 2/L*Sin[(2*\[Pi]*q*x)/L + \[Phi]1]*Sin[(2*\[Pi]*q*x)/L + \[Phi]2]* Sin[(2*\[Pi]*qt*x)/L + \[Phi]3], {x, 0, L}, GenerateConditions -> True]. May 22, 2022 at 19:42
• It is clear without firing a cannon at sparrows the only possible singularities are at the zeroes of the denominator (4 \[Pi] qt (-4 q^2 + qt^2)). May 22, 2022 at 19:50