6

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?

bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
Melquiades
  • 163
  • 5

2 Answers2

6

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.

bbgodfrey
  • 61,439
  • 17
  • 89
  • 156
5

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}})}$$

enter image description here

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*)
mattiav27
  • 6,677
  • 3
  • 28
  • 64