# NIntegrate and Integrate giving different results for ill-behaved function

I'm trying to integrate the following function with Mathematica 8:

$$I(a,b)= \int_0^1 \mathrm{d}x\int_0^1\mathrm{d}y \,\theta(1-x-y) \frac{1}{x a^2-y(1-y)b^2},$$ where $\theta$ is the Heaviside function. However, I find different results with Integrate or NIntegrate and I don't understand why.

More specifically, for a=100 and b=90:

NIntegrate[HeavisideTheta[1 - x - y]/(x a^2 - y (1 - y) b^2), {x, 0, 1}, {y, 0, 1}]


gives

0.0000927294,


while

Integrate[HeavisideTheta[1 - x - y]/( x a^2 - y (1 - y) b^2), {x, 0, 1}, {y, 0, 1}, PrincipalValue -> True]


gives

+0.0000600275+0.000314159 I.


What is the correct result? Why does Integrate give a complex result?

• Various closely related questions have been asked several times, I remember these questions since I answered them: Why does Integrate declare a convergent integral divergent?, Bug in mathematica analytic integration? and Symbolic integration error. These should clarify the issue. Feb 9, 2015 at 0:22
• V10.0.2.0 behaves differently. NIntegrate gives failure-to-converge error messages before giving the answer 0.0000927294. Adding the option AccuracyGoal -> 20 gives 0. as the answer with no error messages. Integrate, on the other hand, returns only a partially integrated solution. Feb 9, 2015 at 5:11
• We deal with an improper double integral: look at the result of ContourPlot[{(x 100^2 - y (1 - y) 90^2) == 0, 1 - x - y == 0}, {y, 0, 1}, {x, 0, 1}]. It's unclear whether the integral under consideration txists at all. @bbgodfrey counts the iterated integral. Oct 31, 2020 at 6:53

Reversing the order of integration produces a solution:

ans= Integrate[HeavisideTheta[1 - x - y]/(x 100^2 - y (1 - y) 90^2), {y, 0, 1}, {x, 0, 1},
PrincipalValue -> True]
(* Log[(100*10^(38/81))/(81*19^(19/81))]/10000 *)

N[ans]
(* 0.000060027526501455836 *)


Solutions of this sort are what I would expect based on outlining a pencil-and-paper derivation.

Attempting to solve the integral numerically produces error messages. For instance,

NIntegrate[HeavisideTheta[1 - x - y]/(x 100^2 - y (1 - y) 90^2), {y, 0, 1}, {x, 0, 1},
MinRecursion -> 30, MaxRecursion -> 60]


produces

NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small. >>
NIntegrate::eincr: The global error of the strategy GlobalAdaptive has increased more than 2000 times. The global error is expected to decrease monotonically after a number of integrand evaluations. Suspect one of the following: the working precision is insufficient for the specified precision goal; the integrand is highly oscillatory or it is not a (piecewise) smooth function; or the true value of the integral is 0. Increasing the value of the GlobalAdaptive option MaxErrorIncreases might lead to a convergent numerical integration. NIntegrate obtained 0.0000591404833446268 and 0.0003356281358114434 for the integral and error estimates. >>
(* 0.0000591404833446268 *)


Of course,NIntegrate has many options, and one or more of them may produce an acceptable answer.

TL;DR Use HeavisideTheta's properties before integration.

This is my strategy.

First the HeavisideTheta gives you the following integration limits:

1. $$0\leq y \leq 1-x \qquad \& \qquad 0\leq x \leq 1$$

2. $$0\leq x \leq 1-y \qquad \& \qquad 0\leq y \leq 1$$

In both cases I used Integrate first then NIntegrate.

In the first case I could not integrate numerically beacuse of a divergence in the result of the symbolical integration; in fact I had:

Integrate[1/(100^2 x - y (1 - y) 90^2), {y, 0, 1 - x}]

Plot[%,{x,0,1}, PlotRange -> Full]


$$\dfrac{(\arctan(\frac{9}{\sqrt{-81 + 400 x}}) + \arctan({\frac{9 (1 - 2 x)}{\sqrt{-81 + 400 x}})}}{(\frac{450}{\sqrt{-81 + 400 x}})}$$

Maybe someone can come up with a solution also for this case.

In the second case I obtained for the symbolical integration:

Integrate[1/(100^2 x - y (1 - y) 90^2), {x, 0, 1 - y}]
Integrate[%, {y, 0, 1}]
N@%
(*0.0000600275 - 0.000314159 I*)


and for the numerical integration:

NIntegrate[
Integrate[1/(100^2 x - y (1 - y) 90^2), {x, 0, 1 - y}], {y, 0, 1}]
(*0.0000600275 - 0.000314159 I*)