# Incorrect definite integral

I have a definite integral that depends on a parameter. The integrand is nonsingular when $h$ is real and $> 1$. When I evaluate the integral with assumption $h > 1$, then Mathematica returns a result. When I evaluate the integral without assumptions, and then refine the result with the assumption, the result is different. Numerical integration for a specific value of $h$ suggests that the first result is correct and the second one is incorrect. Maybe I am just being stupid today, but how can this be?

itg[q_] :=  Cos[q]^2 ( 1 + 1 / (Cos[q]^2 + (Sin[q] + 2 h)^2 ))

Limit[itg[q], h -> Infinity]
(* Cos[q]^2 *)

Integrate[ Cos[q]^2, {q, 0, 2 Pi}]
(* π *)

ans1 = Integrate[ itg[q], {q, 0, 2 Pi}, Assumptions -> h > 1]
(* π + π/(4 h^2) *)

ans2 =  Evaluate[ Integrate[ itg[q], {q, 0, 2 Pi}]]
(* ConditionalExpression[π + 1/8 (4 + 1/h^2) π,
Re[1/h + 4 h] > 4 || Re[1/h + 4 h] < -4 ||
1/h + 4 h ∉ Reals] *)

Refine[ ans2, Assumptions -> h > 1]
(* π + 1/8 (4 + 1/h^2) π *)

NIntegrate[ itg[q] /. h -> 3, {q, 0, 2 Pi}]
(* 3.22886 *)

ans1 /. h -> 3.
(* 3.22886 *)

ans2 /. h -> 3.
(* 4.75602 *)


• Please provide Mathematica code that can be copied by other people, instead of screenshots. People will be more willing to help you if they don't have to retype everything you did. – J. M.'s torpor Mar 9 '18 at 22:51
• Thank you, I have added the code – mekanix Mar 9 '18 at 23:23
• @mekanix This function itg looks fascinating. What kind of physical/mathematical problem gave rise to the function? Thanks. – Vixillator Mar 12 '18 at 10:19
• This integral evaluates the added mass due to acceleration of a circular cylinder adjacent to an infinite plane through an ideal fluid. The integration is for motion parallel to the infinite plane, i.e. the "1,1" component of the added mass tensor. The potential is for a cylinder moving with unit velocity in parallel to the wall plus its image. The parameter h is the distance to the plane, relative to the radius of the circle – mekanix Mar 12 '18 at 16:38
• The subject of what is happening here has been discussed previously in this Wolfram Blog entry. – J. M.'s torpor Mar 14 '18 at 7:51

This is a bug in Integrate. Mathematica does handle the branch cut wrong, if you don't give conditions for h.

With conditions for h, you get the right integral.

itg[q_, h_] = Cos[q]^2 (1 + 1/(Cos[q]^2 + (Sin[q] + 2 h)^2));

int1[h_] = Integrate[itg[q, h], {q, 0, 2 π}, Assumptions -> -1/2 < h < 1/2]

(*   ConditionalExpression[2 π, h != 0]  *)

int1[0] = Integrate[itg[q, 0], {q, 0, 2 π}]

(*   2 π   *)

int2[h_] = Integrate[itg[q, h], {q, 0, 2 π}, Assumptions -> h >= 1/2]

(*   π + π/(4 h^2)   *)

int3[h_] = Integrate[itg[q, h], {q, 0, 2 π}, Assumptions -> h <= -1/2]

(*   π + π/(4 h^2)   *)

int[h_] =
Piecewise[{{π + π/(4 h^2),
h <= -1/2}, {2 π, -1/2 < h < 1/2}}, π + π/(4 h^2)];

ans2[h_] = Integrate[itg[q, h], {q, 0, 2 π}]

(*   ConditionalExpression[
1/8 (12 + 1/h^2) π, (1/h + 4 h \[NotElement] Reals ||
Re[1/h + 4 h] >= 4 || Re[1/h + 4 h] <= -4) && h != -(1/2) &&
2 h != 1]    *)


Compare with NIntegrate and int[h] . ans2[h] is wrong.

nint[h_] := NIntegrate[itg[q, h], {q, 0, 2 π}, MaxRecursion -> 30]

{Plot[{nint[h], ans2[h]}, {h, -2, 2}, PlotRange -> {0, 7},
PlotStyle -> {Blue, Red}],
Plot[int[h], {h, -2, 2}, PlotRange -> {0, 7}],
Plot[ri[h], {h, -2, 2}, PlotRange -> {0, 7}]}


Do indefinite integration, to see where the error arises.

mint[q_, h_] = Integrate[itg[q, h], q, Assumptions -> 0 <= q <= 2 π]

(*   (q + 12 h^2 q + (2 - 8 h^2) ArcTan[(Cos[q/2] + 2 h Sin[q/2])/(
2 h Cos[q/2] + Sin[q/2])] + 4 h Cos[q] + 4 h^2 Sin[2 q])/(16 h^2)   *)


Manipulate shows a jump due to the branch cut of ArcTan, with position depending on h. Mathematica handles this incorrectly.

Manipulate[Plot[mint[q, h], {q, 0, 2 π}], {h, -2, 2}]


For comparison, do indefinite integration with the rule-based-integrator of Albert D. Rich:

rint[q_, h_] = Int[itg[q, h], q]

(*   q/2 + 1/16 (4 + 1/h^2) q +
1/8 (4 - 1/h^2) ArcTan[(4 h + (1 + 4 h^2) Tan[q/2])/(1 - 4 h^2)] +
Cos[q]/(4 h) + 1/2 Cos[q] Sin[q]   *)


Branch cut is now always at q == π and can be eliminated with limit calculation.

Manipulate[Plot[rint[q, h], {q, 0, 2 π}], {h, -2, 2}]

ri[h_] = rint[2 π, h] - Limit[rint[q, h], q -> π, Direction -> -1] +
Limit[rint[q, h], q -> π, Direction -> 1] - rint[0, h] //
FullSimplify

(*   ((1 - Sqrt[(1 + 4 h^2)^2/(-1 + 4 h^2)^2] +
8 h^2 (2 + Sqrt[(1 + 4 h^2)^2/(-1 + 4 h^2)^2] -
2 h^2 (-3 +
Sqrt[(1 + 4 h^2)^2/(-1 + 4 h^2)^2]))) π)/(8 (h^2 + 4 h^4)
)   *)


Applying conditions to this general integral simplifies. ri[h] is shown in the picture above.

ri[h] // FullSimplify[#, h >= 1/2 || h <= -1/2] &

(*   π + π/(4 h^2)   *)

ri[h] // FullSimplify[#, -1/2 <= h <= 1/2] &

(*   2 π   *)


Since the other answers have already explained what is going on, let me supply an expression for the antiderivative that is continuous over the real line (with the assumption that h >= 1):

sol[h_, q_] := (q + 4 h^2 q + 2 h Cos[q] (1 + 2 h Sin[q]))/(8 h^2) -
1/8 (4 - 1/h^2) ArcTan[Cos[q]/(2 h + Sin[q])]


Verify that it is an antiderivative:

D[sol[h, q], q] == Cos[q]^2 (1 + 1/(Cos[q]^2 + (Sin[q] + 2 h)^2)) // Simplify
True


The definite integral can now be evaluated in the usual way:

sol[h, 2 π] - sol[h, 0] // Simplify
π + π/(4 h^2)

• Do you remember how you got this continuous antiderivative? Why is ArcTan[(Cos[q] + 2*h*Sin[q])/(2*h*Cos[q] + Sin[q])] == q + ArcTan[Cos[2*q]/(2*h + Sin[2*q])] for 0<q<pi/2 ? – Andreas Feb 24 at 19:36
• @Andreas, if memory serves, it was a matter of adding q - ArcTan[Tan[q]] and then applying the customary addition formulae. – J. M.'s torpor Feb 24 at 19:52
• crazy idea! I wouldn't have thought of that... – Andreas Feb 24 at 20:26

Sometimes it helps to plot the function before starting to work with it. The given integrand is a function that has two singularities at {q, h}={Pi/2,-1/2} and {q, h}={3 Pi/2,1/2}, as can be seen from the plots and code below. With an indefinite integration (ans2), the integral includes the two singularities plus the rest of the domain. But when we integrate with the assumption of h>1, the integration domain does not include the two singularities, and therefore the magnitude of the definite integral (ans1) will be less than ans2. The two singularities introduce a positive bias to the integral function, and the difference (ans2-ans1) is Pi/2 as h approaches infinity, but this difference is less than Pi/2 for small values of h (say h=3), as is shown in the last plot below.

(* INTEGRAND *)
itg[q_]:=Cos[q]^2 (1+1/(Cos[q]^2+(Sin[q]+2 h)^2))
Plot3D[itg[q],{q,0,2 Pi},{h,-1,1}]


ContourPlot[itg[q],{q,0,2 Pi},{h,-5,5}]


(* 3D PARAMETRIC PLOT OF INTEGRAND *)
itgx := itg [q] Cos[q]; itgy := itg[q] Sin[q]; itgz := h;
ParametricPlot3D[{itgx, itgy, itgz}, {q, 0, 2 Pi}, {h, -2, 2}]


 (* SINGULARITY AT {q,h}={3 Pi/2, 1/2} *)
L1=Limit[itg[q],{h->1/2}]//Simplify;Plot[L1,{q,2.999 Pi/2,3.001 Pi/2},PlotStyle->Thickness[.01]]


 (* SINGULARITY AT {q,h}={Pi/2, -1/2} *)
L2=Limit[itg[q],{h->-1/2}]//Simplify;Plot[L2,{q,0.999 Pi/2,1.001Pi/2},PlotStyle->Thickness[.01]]


 (* INTEGRATE THE FUNCTION AND CHECK LIMITS *)
ans1=Integrate[itg[q],{q,0,2 Pi},Assumptions->h>1];
Limit[ans1,h->Infinity] (* Output: Pi *)
ans2=Evaluate[Integrate[itg[q],{q,0,2 Pi}]]//Normal;
Limit[ans2,h->Infinity] (* Output: 3 Pi/2 *)
ans2ans1=ans2-ans1//Simplify;
Limit[ans2ans1,h->Infinity]; (* Output: Pi/2 *)
Plot[ans2ans1,{h,.1,10},PlotStyle->Thickness[.01]]


Here is more on the two singularities analysis (see comments by @Akku14):

CASE: {q, h}={3 Pi/2,1/2}

(* CASE: {q, h}={3 Pi/2,1/2} *)
L1 = Limit[itg[q], {h -> 1/2}]
(* Out: Cos[q]^2 (1 + 1/(Cos[q]^2 + (1 + Sin[q])^2)) *)
L1 /. q -> 3 Pi/2
During evaluation of In[98]:= Power::infy: Infinite expression 1/0 encountered.
During evaluation of In[98]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered.
(* Out: Indeterminate *)


CASE: {q, h}={Pi/2,-1/2}

(* CASE: {q, h}={Pi/2,-1/2} *)
L2 = Limit[itg[q], {h -> -1/2}]
(* Out:  Cos[q]^2 (1 + 1/(Cos[q]^2 + (-1 + Sin[q])^2)) *)
L2 /. q -> Pi/2
During evaluation of In[100]:= Power::infy: Infinite expression 1/0 encountered.
During evaluation of In[100]:= Infinity::indet: Indeterminate expression 0 ComplexInfinity encountered.
(* Out: Indeterminate *)


ZOOM IN TO SINGULARITIES

(* ZOOM IN TO SINGULARITY AT {q,h}={3 Pi/2, 1/2} *)
zoomq = 0.9999999; q1 = zoomq (3 Pi/2); q2 = (3 Pi/2)/zoomq;
L1 = Limit[itg[q], {h -> 1/2}] // Simplify; Plot[L1, {q, q1, q2}, PlotRange-> Full, PlotPoints -> 200]


(* ZOOM IN TO SINGULARITY AT {q,h}={Pi/2, -1/2} *)
zoomq = 0.9999999; q1 = zoomq (Pi/2); q2 = (Pi/2)/zoomq;
L2 = Limit[itg[q], {h -> -1/2}] // Simplify; Plot[L2, {q, q1, q2}, PlotRange-> Full, PlotPoints -> 200]


The integrand itg is a bivariate function of the independent variable q (angle), and parameter variable h. The function itg has two point singularities (discontinuities) at {q,h}={Pi/2,-1/2) and {q,h}={3 Pi/2,1/2) where the integrand blows up to infinity. The definite integral w.r.t. q (from 0 to 2Pi) of this function takes a finite value at every value of h. However, because the integrand is discontinuous, as will be shown below it is not possible to compute a single analytic expression for the integral that is valid exactly for the entire span of parameter h, from -Infinity to +Infinity. Therefore, before the integration process, it is necessary to divide the domain of integration into three regions (A, B and C), after which we can compute an analytic expression for the integral that is valid for each region. It will be shown below that in fact the region h<-1/2 and h>1/2 are identical, which implies that only two expressions are required for the two integrals that suffice to cover all three regions. It is important to note that, as will be illustrated below, if we attempt to compute a single analytic expression for the definite in q, and indefinite in h integral over the entire domain (by disregarding the discontinuities), Mathematica will compute a simple average of the integral expressions from two regions, which means that this average (ans2) will not be exact in any of the individual regions. Technically this is not a bug with Mathematica but rather an outcome of Mathematica yielding the best average result from an ill-posed task. The issue here is that attempting to obtain a single analytic expression for the improper integral of a discontinuous integrand is not always possible, unless we are able to find a single integral expression that is able to account simultaneously for all the regions adjoining the discontinuities. What follows is the analysis code that illustrates the above:

(* The integrand itg: *)
itg[q_] := itg[q] = Cos[q]^2 (1 + 1/(Cos[q]^2 + (Sin[q] + 2 h)^2));

(* ZONE A {h>+1/2}: *)
ansA = ans1 = Integrate[itg[q], {q, 0, 2 Pi}, Assumptions -> h > 1/2] (* Note: Assumptions\[Rule]h>1 also yields the same result *)


OUT: [Pi] + [Pi]/(4 h^2)

(* ZONE B {h>-1/2 && h<+1/2}: *)
ansB = Integrate[itg[q], {q, 0, 2 Pi}, Assumptions -> {h > -1/2 && h < 1/2}] //Normal(* Note: Assumptions\[Rule]h>1/2 also yields the same result *)


OUT: 2 [Pi]

(* ZONE C {h<-1/2}: *)
ansC = Integrate[itg[q], {q, 0, 2 Pi}, Assumptions-> h < -1/2] (* Note: ansC is found to be identical to ansA(=ans1) *)


OUT: [Pi] + [Pi]/(4 h^2)

(* Now let us compute a single naive "average" expression for the integral that we hope will be valid over the entire domain of h from h=-Infinity to h=Infinity, ignoring the singularities *)
ansavg = (ansA + ansB)/2 // Simplify (* We do not need ansC because ansC = ansA *)


OUT: 1/8 (12 + 1/h^2) [Pi]

(* Now let us compute the indefinite integral over the entire domain, by ignoring the discontinuous nature of the integrand *)
ans2 = Evaluate[Integrate[itg[q], {q, 0, 2 Pi}]] // Normal // Simplify


OUT: 1/8 (12 + 1/h^2) [Pi]

(* Check to see if ansavg is the same as ans2: *)

ansavg == ans2


OUT: True (* They are identical! *)

It has been shown above that Mathematica computed the indefinite integral of a discontinuous function (itg) as a simple average between the integral expressions of the two regions adjoining the discontinuities (singularities).

• There are no singularities at {q, h}={Pi/2,-1/2} and {q, h}={3 Pi/2,1/2}. Singularities are points, where function value go to infinity. This is not the case here. Function value for h==1/2 approaching q->3 Pi/2 is 1 and 0 for all other h. Limit[itg[q, 1/2], q -> 3 Pi/2] . Therefore you don't get a bias to the integral. The difference for h>1/2 and h<1/2 comes from the normal integration, as my answer and NIntegrate show. – Akku14 Mar 11 '18 at 17:55
• Both the additional analysis in my post, and the plots (especially the plots), indicate that there are indeed two singularities of itg at {q, h}={Pi/2,-1/2} and {q, h}={3 Pi/2,1/2}. Your order of taking the limits is different from mine. For what it's worth, this "interchange of limits" issue appears to be a known major unsolved challenge of mathematical analysis.Re: en.wikipedia.org/wiki/Interchange_of_limiting_operations – Vixillator Mar 11 '18 at 19:27
• First: The function is not defined at these two points. But integration with exclusion of only one or two single points does not have any influence to the integral. Second you compare ans1 with ans2. Since ans2 is totaly wrong, this makes no sence and produces no further insight. – Akku14 Mar 11 '18 at 20:01
• Interchange of limits is not relevant here, because you only integrate over q at certain h. – Akku14 Mar 11 '18 at 20:06
• Doing a series expansion at q=Pi/2 and q=3 Pi/2 perhaps shows the singularities more clearly: Series[Cos[q]^2 (1 + 1/(Cos[q]^2 + (Sin[q] + 2 h)^2)), {q, Pi/2, 3}]//Normal Out: (1+1/(1+4 h+4 h^2)) (-([Pi]/2)+q)^2 Series[Cos[q]^2 (1 + 1/(Cos[q]^2 + (Sin[q] + 2 h)^2)), {q, 3 Pi/2, 3}]//Normal Out: (1+1/(1-4 h+4 h^2)) (-((3 [Pi])/2)+q)^2 – Vixillator Mar 11 '18 at 20:22