1
$\begingroup$

I would like to integrate the following,

$$f(x,y) = \frac{1}{(a - x)^2 + b (y-x^2)^2},$$

where $a>0$ and $b>0$ between definite bounds to obtain an analytic result ideally. I've tried the following,

Assuming[{a > 0, b > 0, xlower < xupper, ylower < yupper},
  Integrate[1 / ((a - x)^2 + b (y - x^2)^2), {x, xlower, xupper},
            {y, ylower, yupper}]]

But the calculation hangs and doesn't seem to return anything. It may be that an analytic result for such an integral does not exist and so I am unsurprised in a way.

I've also tried to numerically integrate this function for $a=1$ and $b=100$ (it's the inverted Rosenbrock function which is used to test optimisation algorithms and so I'm not surprised I'm running into trouble) but I seem to get hugely variable results dependent on algorithm settings.

The following returns 11.74... but this result is different to Python's default double integrator and so am not sure what to trust.

NIntegrate[1/((1 - x)^2 + 100 (y - x^2)^2), {x, -2, 4}, {y, -1, 6}, 
   AccuracyGoal -> 200000, MaxRecursion -> 100000, MaxPoints -> 10000000]

Any tips on how to calculate this integral either analytically or numerically would be well received!

Edit: I've tried Monte Carlo and Quasi-Monte Carlo approaches (and also adaptive approaches) and all methods seem to return different results.

For example,

NIntegrate[1/((1 - x)^2 + 100 (y - x^2)^2), {x, -2, 4}, {y, -1, 6}, 
   AccuracyGoal -> 200000, MaxRecursion -> 100000, 
   MaxPoints -> 10000000, Method -> "MonteCarlo"]

yields 5.19... but

NIntegrate[1/((1 - x)^2 + 100 (y - x^2)^2), {x, -2, 4}, {y, -1, 6}, 
   AccuracyGoal -> 200000, MaxRecursion -> 100000, 
   MaxPoints -> 10000000, Method -> "AdaptiveQuasiMonteCarlo"]

yields 1.89...

This integral is obviously pretty pathological and is likely why it is used to test optimisers!

Edit 2: I see that the function when $a=1$ and $b=100$ has a point $(x,y)=(1,1)$ where $f\rightarrow\infty$ and so perhaps the integral is not normalisable...?

Improve this question
$\endgroup$

2 Answers 2

Reset to default
1
$\begingroup$

I don't know if your work requires you to integrate over discontinuities, but Mathematica can provide a symbolic result for your double integral if you restrict your parameters and domain to avoid them:

y1 = Simplify[Integrate[1/((a - x)^2 + b (y - x^2)^2), y]];
y2 = Simplify[(y1 /. y -> yupper) - (y1 /. y -> ylower)];
x1 = Simplify[Integrate[y2, x]];
x2 = Simplify[(x1 /. x -> xupper) - (x1 /. x -> xlower)];

x2 (expressed with root objects) is quite long, so I won't show its output here; it has a leaf count of 5903 (determined with LeafCount@x2. You can convert the root objects to radicals using ToRadicals@x2, but this makes the expression much longer: leaf count = 96,183.

Note, however, that the above relies on the Fundamental Theorem of Calculus, and is thus only valid for values of a, b, ylower, yupper, xlower, and xupper over which both y1 and x1 are continuous. [You can check this by plotting, being careful to use a MaxRecursion sufficiently high to avoid missing discontinuities.]

Based on my own plots, the above conditions seem to be satisfied with, for instance: a=1, b=100, ylower=10, yupper=20, xlower=3, xupper=20.*

So:

a = 1
b = 100
y1 = Simplify[Integrate[1/((a - x)^2 + b (y - x^2)^2), y]];
y2 = Simplify[(y1 /. y -> 10) - (y1 /. y -> 20)];
x1 = Simplify[Integrate[y2, x]];
x2 = Simplify[(x1 /. x -> 3) - (x1 /. x -> 20)];
N[x2, 1000];
Precision[%]
NumberForm[%%, 6]

1000.

0.144867 + (0.*10^-1000) I

And once you've used the above check to find values that avoid discontinuities, Mathematica works normally:

int = Integrate[1/((1 - x)^2 + 100 (y - x^2)^2), {x, 3, 20}, {y, 10, 20}]
N[int, 1000]
Precision[%]
NumberForm[%%, 6]

1000

0.144867 + (0.*10^-1000) I

Clarification added based on discussion in Comments Section: Here I've used arbitrary-precision arithmetic to reliably numericize the symbolic result. With six significant figures, I could have obtained the same real component with N[x2, 6]. However, I chose N[x2, 1000] to show that the imaginary part can be made arbitrarily small with a sufficiently high-precision calculation (more specifically, it is always smaller than the least significant figure), indicating it is an artifact of the conversion of the symbolic result to numerical form, and can thus be ignored.

You can compare the above with the Monte Carlo results from NIntegrate. Since these are Monte Carlo calculations, each run will of course yield a different result. Thus, here, using WorkingPrecision->1000 simply means that particular Monte Carlo run was calculated to a precision of 1000. It doesn't mean the result itself has that precision relative to the true value.

For example, running the following code three times, I get values of 0.144902, 0.144897, and 0.144942:

NIntegrate[1/((1 - x)^2 + 100 (y - x^2)^2), {x, 3, 20}, {y, 10, 20}, 
PrecisionGoal -> Automatic, MaxRecursion -> 100000, 
MaxPoints -> 10000000, WorkingPrecision -> 1000, 
Method -> "AdaptiveQuasiMonteCarlo"];
Precision[%]
NumberForm[%%, 6]
Improve this answer
$\endgroup$
6
  • $\begingroup$ Extreme precision is not required. With the symbolic definition of x2: xval = x2 /. {a -> 1, b -> 100, yupper -> 10, ylower -> 20, xupper -> 3, xlower -> 20} // N[#, 20] & // Chop evaluates to 0.14486739879335968388 Or a reasonable approximation is xval2 = x2 /. {a -> 1, b -> 100, yupper -> 10, ylower -> 20, xupper -> 3, xlower -> 20} // RootApproximant $\endgroup$
    – Bob Hanlon
    Commented Feb 16, 2019 at 2:49
  • $\begingroup$ @BobHanlon Well, sure, but one could make the same argument that what you wrote is likewise unnecessarily precise. Indeed, all that's needed for comparison to the numerical results is N[x2, 6]. I just chose N[x2, 1000] so that the OP could be reassured that the imaginary component could be made arbitrarily small, and thus could be considered an artifact that can be ignored, enabling direct comparison to the numerical results. $\endgroup$
    – theorist
    Commented Feb 16, 2019 at 4:52
  • $\begingroup$ @BobHanlon ...continued: The interesting question for me, and I'd like to hear your take, is this: Normally, I consider a numericization of the symbolic integration (if one can get a symbolic integration) to represent the reference. Thus, if there is a small discrepancy between that and what NIntegrate gives, my default assumption is that NIntegrate's result contains a small calculation error. But when the symbolic result is this large (leaf count in the 1000s), is it possible that it's the numericization of the symbolic result that has the error (here I'm referring to the real component)? $\endgroup$
    – theorist
    Commented Feb 16, 2019 at 5:04
  • $\begingroup$ As long as you use arbitrary precision (rather than machine precision) then Mathematica will track the precision and Precision should be fairly good at telling you the effect of the large number of calculations. $\endgroup$
    – Bob Hanlon
    Commented Feb 16, 2019 at 5:09
  • 1
    $\begingroup$ Only the "numericization of the symbolic result." $\endgroup$
    – Bob Hanlon
    Commented Feb 16, 2019 at 5:19
0
$\begingroup$

This integral diverges since there is a pole at $(x,y)=(a,a^2)$ where $f\rightarrow\infty$.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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