I'm trying to perform nth-Order Legendre-Gauss Quadrature on an integral of mine. The algorithm I've written for the quadrature works on simple integrals, like the following.



However, the problem I'd like to solve is particularly ugly and the results from quadrature don't converge no matter how high in n I go (I've gone to 256th order!). The results are (seemingly) random across several orders of magnitude. I'm wondering if anyone has any insight into this issue and might have ideas on how to attack it.

The integral is used to calculate the intensity in x-ray diffraction of a lattice of atoms. The lattice is perfect except for a collection of dislocations. There are five types of dislocations detailed below, indexed by , which have an associated Burger's vector, , and gives a displacement vector on the lattice. Only the z-component of the displacement vector, here defined as is needed. Note that is a physical parameter (a constant for the lattice) called Poisson's ratio, are also constants.

Screw Dislocations
Type 1: \vec{B}^{(1)} = <0 \ 0 \ c>
Type 2: \vec{B}^{(2)} = <0 \ 0 \ -c>


Edge Dislocations
Type 3: \vec{B}^{(3)} = < a \ 0 \ 0>
Type 4: \vec{B}^{(4)} = < \frac{-a}{2} \ \frac{\sqrt{3}a}{2} \ 0>
Type 5: \vec{B}^{(5)} = < \frac{-a}{2} \ \frac{-\sqrt{3}a}{2} \ 0>


Note that
Now, the total (z-component of the) displacement vector in the lattice is



is the number of -type disolcations, and is the intersection of the th -type dislocation with the sample surface.

Finally, we can get to the integral. and are vectors that are constant in space, is the sample's thickness, and is the coherent area, which is simplified as a rectangle bounded by [x_{0},x_{1},y_{0},y_{1}].



Why doesn't the quadrature work on my integral?!! I've outlined Gauss Quadrature below.

Quadrature
An integral of the function can (in theory) be approximated by the following. \{\theta\} are the roots of the nth Legendre polynomial, and \{w\} are the associated weights.

I = \int_{x_{0}}^{x_{1}} \int_{y_{0}}^{y_{1}} \int{z_{0}}^{z_{1}} dxdydz f(x,y,z)
I \approx x_{a}y_{a}z_{a} \sum_{i,j,k}^{n} w_{i}w_{j}w_{k} f(x_{a}\theta_{i} + x_{b}, y_{a}\theta_{j} + y_{b}, z_{a}\theta_{k} + z_{b})

x_{a} = \frac{x_{1} - x_{0}}{2} \ , \ x_{b} = \frac{x_{1} + x_{0}}{2} \ , \ \mbox{etc.}