# Integrate returns unexpected result

Consider the following function $$g(x,y):= \frac{1}{( (1+y)^2+x^2 )( 1+ax^2y^2 )^2}$$, where I assume that $y\geq 0$ and $a\in (0,1]$ is a parameter. When I try to evaluate the integral $\int g(x,y)\,dx$ using Mathematica I get the following result:

$$\frac{\frac{axy^2(-1 + ay^2(1 + y)^2)}{(1 + ax^2y^2)} + \sqrt{a}y(-3 + ay^2(1 + y)^2)\arctan(\sqrt{a}xy) + \frac{2\arctan(\frac{x}{1 + y} )}{ 1 + y} }{2(-1 + ay^2(1 + y)^2)^2}$$ This result cannot be the correct primitive of $g(\cdot,y)$, since on the one hand the denominator is always has a positive root but the numerator is positive, and on the other hand we have that $$0\leq g(x,y)\leq \frac{1}{1+x^2 }$$ which implies that $g(\cdot,y)$ is Riemann integrable for all $y\geq 0$.

Could someone please tele me how to obtain a correct answer to the problem? This would be very much appreciated!

Code:

Integrate[1/(((1+y)^2+x^2)*(1+a*x^2*y^2)^2),x]


Output:

((a x y^2 (-1 + a y^2 (1 + y)^2))/(1 + a x^2 y^2) + Sqrt[a] y (-3 + a y^2 (1 + y)^2) ArcTan[Sqrt[a] x y] + ( 2 ArcTan[x/(1 + y)])/(1 + y))/(2 (-1 + a y^2 (1 + y)^2)^2)

-
Please post you expressions as Mathematica code, not in TeX. (Use backticks instead of dollar signs.) It makes it easier for those who would help you with your problem. –  Michael E2 Aug 8 '14 at 20:49
please provide the mathematica command you used? –  chris Aug 8 '14 at 20:50
Sure! Here's the code: Integrate[1/(((1+y)^2+x^2)*(1+a*x^2*y^2)^2),x] –  user19057 Aug 8 '14 at 20:53
Please also post the output you get in Mathematica form so it's easy to compare. The output I get in v10 looks like what you posted, and D[%,x] gives the expected result. –  Szabolcs Aug 8 '14 at 21:18
@Szabolcs same with v9 –  chris Aug 8 '14 at 21:22

I said in a comment that one could "factor" out ((-1 + a y^2 (1 + y)^2)^2 from the numerator. What I meant was that numerator is $O(((-1 + a y^2 (1 + y)^2)^2)$. The integral is a complicated expression, so the easiest way to examine it, it seemed to me, is to look at the coefficients of the power series expansion about the real roots of the denominator.

integral = Integrate[1/(((1 + y)^2 + x^2)*(1 + a*x^2*y^2)^2), x];


There are two real roots, if 0 < a <= 1:

Solve[(-1 + a y^2 (1 + y)^2) == 0, y]
(*
{{y -> 1/2 (-1 - Sqrt[-4 Sqrt[a] + a]/Sqrt[a])},
{y -> 1/2 (-1 + Sqrt[-4 Sqrt[a] + a]/Sqrt[a])},
{y -> 1/2 (-1 - Sqrt[4 Sqrt[a] + a]/Sqrt[a])},
{y -> 1/2 (-1 + Sqrt[4 Sqrt[a] + a]/Sqrt[a])}}
*)


We see below that the first two coefficients of the series expansions about the two real roots are zero. Thus the integral is bounded for -b < y < b for any positive real b.

root = 1/2 (-1 - Sqrt[4 Sqrt[a] + a]/Sqrt[a]);
ParallelMap[
FullSimplify,
CoefficientList[
Normal@Series[Numerator[integral], {y, root, 2}] /. y -> u + root,
u]
]
(*
{0, 0, 1/2 (4 + Sqrt[a]) Sqrt[a] (<..>)}
*)

root = 1/2 (-1 + Sqrt[4 Sqrt[a] + a]/Sqrt[a]);
ParallelMap[
FullSimplify,
CoefficientList[
Normal@Series[Numerator[integral], {y, root, 2}] /. y -> u + root,
u]
]
(*
{0, 0, (1/(2 + Sqrt[a] + Sqrt[4 Sqrt[a] + a])) 4 (4 + Sqrt[a]) a (<..>)}
*)


Since as has been observed in the comments,

D[integral, x] // Simplify
(* 1/((1 + a x^2 y^2)^2 (x^2 + (1 + y)^2)) *)


it seems to me that the result integral is correct.

-

These singularities in $g(\cdot,y)$ are all removable. There are eight singular values of y, depending only on a, and you can take the limit as $y\rightarrow s$ for each singular value $s$ and obtain limit for all but finitely many values of $x$. Those are also removable by taking a limit. The function is still continuous, smooth, bounded, etc.

The values of $y$ of concern can be found to be:

SP = ReplaceAll[y, Solve[2 (-1 + a y^2 (1 + y)^2)^2 == 0, y]]


This will return a list of eight singular values of y. You can take the limit at each value:

Map[Limit[g[x, y, a], y -> #] &, SP]


and in each case, you obtain an expression that exists except for finitely many values of x. At those values of x, again, you can take limits as x approaches the singular value and obtain the value at this removable singularity.

The reason this limit exists, for reference, is because at any of these singular values, you will have 0 on the top of that big fraction too. You could even use L'Hôpital's rule to obtain the value -- but I'll let Mathematica take care of it (see above). You can verify this with:

G = Numerator[Integrate[1/(((1+y)^2+x^2)*(1+a*x^2*y^2)^2),x]];
Map[ReplaceAll[G,y->#]&,SP]


This gives you nothing but a list of 0s, confirming that these singularities might be removable. Taking limits (as above) proves that they are.

And, as many have stated, you can verify this integral by taking the derivative and everything works out. Mathematica generally treats removable singularities like they are irrelevant -- if it could eliminate them, it would, but here it's not easy to do that with elementary functions.

-