# Symbolic Integration gives different results by changing the constants. Is it a bug?

Problem existing in Version 10.4.1 or earlier and persisting through Version 11.3.0

I have the following integral (I know it could be simplified, but this particular form is important for the question) $$\int_0^\pi \frac{ \Delta p^2 (1-F) \sin^2(\theta -\theta_r+ \Omega\tau )}{\pi \left(1-\Delta p^2 \cos ^2(\theta -\theta_r+ \Omega\tau )\right)} \mathrm{d} \theta.$$ I have Mathematica 11.0.1.0 on Ubuntu 16.04 (but I tried on a Mac as well with the same version). If I start a fresh kernel and without assuming anything I compute the integral in the following way:

Integrate[1/π (Δp^2 Sin[Ω τ + θ- θr]^2 (1-F))/(1 - Δp^2 Cos[Ω τ + θ - θr]^2),{θ,0,π}]


I get:

(* 1 - F *)


While, if I use $\Omega \tau = a$ in the following way

Integrate[1/π (Δp^2 Sin[a + θ- θr]^2 (1-F))/(1 - Δp^2 Cos[a + θ - θr]^2),{θ,0,π}]


I get

(* ConditionalExpression[(-1 + F) (-1 + Sqrt[1/(1 - Δp^2)] - Δp^2 Sqrt[1/
(1 - Δp^2)]),Im[θr] == Im[a] && -(π/2) + Re[a] <= Re[θr] <= π/2 + Re[a]] *)


Apart from the result, it is problematic that just by giving a different name to a constant I get different results. Is this a bug?

Now, regarding the result, this is integral can also be done by hand... anyway one can easily see that the first one is wrong because the result obviously has to depend on $\Delta p$, since for $\Delta p = 0$ everything is zero...

Does anybody have a clue?

• I tried on a friend's mathematica 10 and the result is 1-F in both cases... it is the wrong result, but at least it's consistent! – justmyfault Mar 17 '17 at 16:37
• I can reproduce the strange behavior, it looks like a bug to me. – kiara Mar 17 '17 at 17:16
• Sorry i edited the code but forget to put the Ω τ to a in the cos. And i cannot edit this again. – kiara Mar 17 '17 at 17:41
• Don't worry, I fixed it... I'll try to write to Wolfram Support I guess... – justmyfault Mar 17 '17 at 19:27
• I agree. I tried various options to obtain the correct result but without success. It certainly looks like a bug. In fact, some of the conditions that DSolve returns when using a instead of Ω τ even seem unnecessary. Do you plan to report all this to Wolfram, Inc? – bbgodfrey Mar 18 '17 at 11:12

Both values seem to be correct, depending on the parameters; however, the dependence on the parameters is not correctly expressed in the ConditionalExpression result. It may be that the OP has in mind a real integral, but Mathematica, by default, treats the parameters as complex numbers. By a standard application of the Residue Theorem, an integral of the form $(1/2\pi)\int_0^{2\pi} R(\cos t,\sin t)\,dt$, where $R(x,y)$ is a rational function, is equal to the sum of the residues of $${1\over z}\cdot R\left({z+z^-1 \over 2},{z-z^-1 \over 2i}\right)$$ over the poles in $|z|<1$. The OP's integral is of this form ($2\theta=t$): $$\int_0^\pi \frac{ \Delta p^2 (1-F) \sin^2(\theta -\theta_r+ \Omega\tau )}{\pi \left(1-\Delta p^2 \cos ^2(\theta -\theta_r+ \Omega\tau )\right)} \; d\theta = \int_0^{2\pi} \frac{ \Delta p^2 (1-F) \sin^2(t/2-t_0)}{2\pi \left(1-\Delta p^2 \cos ^2(t/2-t_0)\right)} \;dt$$ $$\qquad= {1\over2\pi}\int_0^{2\pi} \frac{\Delta p^2 (1-F) (\sin t \sin 2 {t_0}+\cos t \cos 2 {t_0}-1)}{\left({\Delta p}^2+{\Delta p}^2 \sin t \sin 2 {t_0}+{\Delta p}^2 \cos t \cos 2 {t_0}-2\right)} \;dt$$

Here is are plots of the numerical integration for an imaginary value of $t_0$ (code below): Fig. 1. Plots of the integral computed with NIntegrate[] for $\Delta p=x+i\,y$, $F=0$, $t_0=0.1i$; the middle plot shows that the integral equals $1-F = 1$ outside a figure-eight region (where all the poles are inside the unit disk).

Here is an implementation of the above idea:

op = 1/π (Δp^2 Sin[Ω τ + θ - θr]^2 (1 - F))/(1 - Δp^2 Cos[Ω τ + θ - θr]^2) /.
{θr -> t0 + Ω τ, θ -> t/2};
integrand = op/2 // TrigReduce // TrigExpand //
Simplify[#,                 (* Complexity function unwanted trig forms *)
ComplexityFunction -> (LeafCount[#] +
1000 Count[#,
Sin[1/2 (t - 2 t0)] | Sin[t/2 - t0] | Cos[t - 2 t0],
Infinity] &)] &;

intz = Block[{F = 0},
integrand
] /. {Sin[t] -> (z - 1/z)/(2 I), Cos[t] -> (z + 1/z)/2} // Simplify;
poles = Solve[z*Denominator[intz] == 0, z] // Simplify;
residues = Residue[1/z*intz, {z, z /. #}] & /@ poles // FullSimplify;

(* the value of the integral *)
val = 2 Pi*residues.UnitStep[1 - Abs[z /. poles]] // PiecewiseExpand // FullSimplify The third condition cannot happen, but Mathematica does not figure that out. (The two distinct absolute values are reciprocals of each other, and the greater one has a minimum value of 1. This implies the first two pieces cover all cases. But that's difficult to compute.) Fig. 2. Plots of the symbolic value val of the integral for $\Delta p=x+i\,y$, $F=0$, $t_0=0.1i$; the middle plot show that the integral equals $1-F = 1$ outside a figure-eight region (where all the poles are inside the unit disk).

Remark on the computational issue. One reason for the different results that the OP experiences is probably the increase in complexity of using a product of constants Ω τ over a single constant a. I think this is worth emphasizing: Mathematica will try to reduce the conditions and perhaps tries to break down the conditions into terms involving the parts of Ω and τ separately. This might lead to different results than with a simpler a. (It does not, however, explain the incorrectness in the result.) Sometimes, with symbolics it pays to simplify the parameters, such replacing Ω τ by a, or indeed both terms Ω τ - θr by t0. They can then be substituted back at the end of the computation.

Code for figure 1:

Block[{Δp = x + I y, t0 = 0 + t1 I, t1 = 0.1, F = 0},
iplotAbs = Plot3D[
Abs@NIntegrate[
1/π (Δp^2 Sin[t - t0]^2 (1 - F))/(1 - Δp^2 Cos[t - t0]^2),
{t, 0, π}, Method -> "Trapezoidal", PrecisionGoal -> 4,
AccuracyGoal -> 4, MaxRecursion -> 20],
{x, -1.04/t1, 1.04/t1}, {y, -1.04/t1, 1.04/t1},
MaxRecursion -> 3, PlotLabel -> "Modulus of integral",
MeshFunctions -> {#3 &}, AxesLabel -> {x, y}
];
] // AbsoluteTiming
(*  {141.737, Null}  *)

Block[{Δp = x + I y, t0 = 0 + t1 I, t1 = 0.1, F = 0},
iplotArg = Plot3D[
Arg@NIntegrate[
1/π (Δp^2 Sin[t - t0]^2 (1 - F))/(1 - Δp^2 Cos[t - t0]^2),
{t, 0, π}, Method -> "Trapezoidal", PrecisionGoal -> 4,
AccuracyGoal -> 4, MaxRecursion -> 20],
{x, -1.04/t1, 1.04/t1}, {y, -1.04/t1, 1.04/t1},
MaxRecursion -> 3, AxesLabel -> {x, y}, MeshFunctions -> {#3 &},
PlotRange -> All, PlotLabel -> "Argument of integral"
];
] // AbsoluteTiming
(*  {145.414, Null}  *)

GraphicsRow[{iplotAbs,
Show[iplotAbs, PlotRange -> {0.8, 1.2}, ViewPoint -> {0, 3.5, 0.7}],
iplotArg}, ImageSize -> 600]


Code for Figure 2:

Block[{Δp = x + I y, t0 = 0 + t1 I, t1 = 0.1, F = 0},
GraphicsRow[
Plot3D[
#[val],
{x, -1.04/t1, 1.04/t1}, {y, -1.04/t1, 1.04/t1},
MaxRecursion -> 3, MeshFunctions -> {#3 &}, ImageSize -> 300,
AxesLabel -> {x, y}, PlotLabel -> #["val"],
Exclusions -> {Automatic, {x == 0, # === Arg}}
] & /@
{Abs, Arg}, ImageSize -> 600
]
]

• +1 just for the nice looking plots, as the math is too hard for me to follow :) – Nasser Mar 19 '17 at 4:31
• Thank you for the interesting and deep answer... So as far I understand I should fill a bug report to wolfram for the incorrect integral and not for the inconsistency, because it might be expected? – justmyfault Mar 19 '17 at 11:31
• Yes, please file a bug report. From the comments to this older, similar problem, it seems that getting a general-purpose algorithm to do this right is harder to get right than to solve a single instance. But each example probably helps them find improvements. – Michael E2 Mar 19 '17 at 12:57

With suitable Assumptions, the two formulations give equivalent results:

Integrate[1/π (Δp^2 Sin[Ω τ + θ - θr]^2 (1 - F))/(1 - Δp^2 Cos[Ω τ + θ - θr]^2),
{θ, 0, π}, Assumptions -> -1 < Δp < 1 && 0 < θr < π && 0 < Ω τ < π]
(* ConditionalExpression[(-1 + F) (-1 + 2 Sqrt[1 - Δp^2]), π + 2 θr < 2 τ Ω] *)

Integrate[1/π (Δp^2 Sin[a + θ - θr]^2 (1 - F))/(1 - Δp^2 Cos[a + θ - θr] 2),
{θ, 0, π}, Assumptions -> -1 < Δp < 1 && 0 < θr < π && 0 < a < π]
(* ConditionalExpression[(-1 + F) (-1 + 2 Sqrt[1 - Δp^2]), 2 a > π + 2 θr] *)


Both now are wrong! The correct solution is

Integrate[1/π (Δp^2 Sin[θ]^2 (1 - F))/(1 - Δp^2 Cos[θ]^2), {θ, 0, π},
Assumptions -> -1 < Δp < 1]
(* (-1 + F) (-1 + Sqrt[1 - Δp^2]) *)


Setting Ω τ - θr to zero is valid, because the integrand is periodic with period Pi. Incidentally, setting only Ω τ to zero also gives the correct result, although with an extraneous condition.

Integrate[1/π (Δp^2 Sin[θ - θr]^2 (1 - F))/(1 - Δp^2 Cos[θ - θr]^2), {θ, 0, π},
Assumptions -> -1 < Δp < 1 && 0 < θr < π]
(* ConditionalExpression[(-1 + F) (-1 + Sqrt[1 - Δp^2]), 2 θr < π] *)


Clearly, the strange behavior exhibited here and in the question represents one or more bugs.

These results were obtained with

\$Version
(* "10.4.1 for Microsoft Windows (64-bit) (April 11, 2016)" *)


and

(* "11.2.0 for Microsoft Windows (64-bit) (September 11, 2017)" *)


Incidentally, the newer version of Mathematica does not reproduce the ConditionalExpression in the question, instead giving,

(* ConditionalExpression[-I (-1 + F) (-I + Sqrt[-1 + Δp^2] + 2 Sqrt[-1 + Δp^2]
Floor[Arg[-1 + Δp^2]/(4 π)]),
Im[a] == Im[θr] && -(π/2) + Re[a] <= Re[θr] <= π/2 + Re[a]] *)


Identical behavior persists for

(* "11.3.0 for Microsoft Windows (64-bit) (March 7, 2018)" *)


except that it fails in a different way to reproduce the ConditionalExpression in the question, instead giving,

ConditionalExpression[1 - F, C ∈ Integers &&
Im[a] == Im[θr] && -(π/2) <= π C - Re[a] + Re[θr] <= π/2]

• Ok thank you for investigating this problem in detail. I never reported a bug to Wolfram, but it's clear that it might good a be idea to do so. I was not even aware that 11.1 came out by the way. – justmyfault Mar 19 '17 at 11:32
• Things seem to get wild in MMA 11.2.0 for the second expression in your answer. Integrate[1/π (Δp^2 Sin[a + θ - θr]^2 (1 - F))/(1 - Δp^2 Cos[a + θ - θr] 2), {θ, 0, π}, Assumptions -> -1 < Δp < 1 && 0 < θr < π && 0 < a < π] returns a very complicated ConditionalExpression with Cot, Tan and Log in it... – user58955 Sep 19 '17 at 15:25
• @user58955 Thanks, for checking whether this bug persisted. I shall take a close look at the matter next week, when I have time. – bbgodfrey Sep 19 '17 at 15:44
• @user58955 The expression in your comment is missing ^: Cos[a + θ - θr] 2 should be Cos[a + θ - θr]^2. In general, I find no difference between 11.1 and 11.2. – bbgodfrey Sep 27 '17 at 0:59