6
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ 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. $\endgroup$ – Artes Feb 9 '15 at 0:22
  • 1
    $\begingroup$ 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. $\endgroup$ – bbgodfrey Feb 9 '15 at 5:11
  • $\begingroup$ 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. $\endgroup$ – user64494 Oct 31 '20 at 6:53
6
$\begingroup$

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.

$\endgroup$
5
$\begingroup$

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*)
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.