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}]




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


+0.0000600275+0.000314159 I.

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

  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$ – bbgodfrey Feb 9 '15 at 0:19
  • 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

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 *)

(* 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]


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

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}]
(*0.0000600275 - 0.000314159 I*)

and for the numerical integration:

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

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.