1

Let

$$f(x,y) = \exp \left(- \frac{1}{2}a x^2 - \frac{1}{2}c y^2 + bxy \right)$$

where $a,b,c\ge 0$. I want to integrate numerically:

$$\int_{x_0}^{x_1}\mathrm{d}x \int_{y_0}^{y_1}\mathrm{d}y \, f(x,y) x^ny^m$$

where $-\infty < x_0 < x_1 < \infty$, $-\infty < y_0 < y_1 < \infty$, and $n,m\in\{0,1,2\}$.

A naive method has problems when the peak of $f(x,y)$ is far from the rectangle of integration.

Is there a method that I can use?

a06e
  • 1,729
  • 15
  • 22

1 Answers1

2

I just tried the integration algorithm that comes with scipy.integrate for double integrals (dblquad), and it seems to work just fine for your problem

from __future__ import division, print_function 
import numpy as np
from scipy.integrate import dblquad

a = 2
b = 2
c = 1
fun = lambda x, y: np.exp(-a/2*x**2 - b/2*y**2 + c*x*y)*x**3
inte2 = dblquad(fun, 1, 4, lambda y: -5, lambda y: 8)
print(inte2)

and returns

(0.5814009878561697, 8.012117719083942e-09)

The first value is the integral and the second is an estimate of the error on the integral.

nicoguaro
  • 8,500
  • 6
  • 23
  • 49
  • Does this address the question? The question in particular asks about when the peak is far from the domain of integration, whereas your example does not seem to be peaked far outside the domain. – user14717 Mar 17 '18 at 10:03
  • @user14717, since I don't have a sense of what far away means in this question I just tried a random example. – nicoguaro Mar 17 '18 at 14:07