Is My 2D Gaussian Quadrature Algorithm Accurate?

Click For Summary
The discussion centers on the accuracy of a 2D Gaussian quadrature algorithm used to evaluate the integral of a specific function. The calculated result from the Python implementation is significantly different from Wolfram's answer, raising concerns about the algorithm's correctness. Participants confirm the use of appropriate weights and grid points, but there is uncertainty about potential subtle errors in the setup. The conversation emphasizes the importance of verifying the theoretical foundations of the chosen points and weights in Gaussian quadrature. Overall, the accuracy of the algorithm remains in question, suggesting a need for further review and validation.
Shackleford
Messages
1,649
Reaction score
2
Homework Statement
Approximate the following integral using two-point Gaussian Quadrature for double integrals.
Relevant Equations
## I = \sum_{i=0}^{n}\sum_{j=0}^{n} w_i w_j f(x_i, y_j)
##
## \int_{-1}^{1} \int_{-1}^{1} e^{-(x^2 + y^2)} cos(2π (x^2 + y^2)\,dx\,dy ##

## I = \int_{-1}^{1} \int_{-1}^{1}f(x,y) \,dx\,dy \approx \sum_{i=0}^{n}\sum_{j=0}^{n} w_i w_j f(x_i, y_j) ##

## = w_0 w_0 f(x_0, y_0) + w_0 w_1 f(x_0, y_1) + w_1 w_0 f(x_1, y_0) + w_1 w_1 f(x_1, y_1) ##

## w_0 = 1, w_1 = 1 ##
## x_0 = 1/\sqrt 3, x_1 = -1/\sqrt 3, y_0 = 1/\sqrt 3, y_1 = 1/\sqrt 3 ##

Of course, the grid points are

## (1/\sqrt 3, 1/\sqrt 3), (1/\sqrt 3, -1/\sqrt 3), (-1/\sqrt 3, 1/\sqrt 3), (-1/\sqrt 3, -1/\sqrt 3). ##

Wolfram's answer is 0.111328. My crude code in Python is -1.02683...

def GaussQuad2D(x_point_dict, y_point_dict, weight_dict, num_points):
summation = []
for e in np.arange(num_points):
w = weight_dict[e]
x = x_point_dict[e]
for j in np.arange(num_points):
y = y_point_dict[j]
f = w * np.exp(-(x**2 + y**2)) * np.cos(2*np.pi*(x**2 + y**2))
summation.append(f)
return sum(summation)

x_point_dict = {0:(-1/np.sqrt(3)), 1:(1/np.sqrt(3))}
y_point_dict = {0:(-1/np.sqrt(3)), 1:(1/np.sqrt(3))}
weight_dict = {0:(1), 1:(1)}
answer = GaussQuad2D(x_point_dict, y_point_dict, weight_dict, 2)

There's a nested loop in there.
Screen Shot 2019-09-20 at 10.11.31 PM.png
 
Physics news on Phys.org
If I understand your formulas, here is what I got with Maple:
maplesnip.JPG

As you can see, it agrees with your Python. Much simpler in Maple though.
 
LCKurtz said:
If I understand your formulas, here is what I got with Maple:
View attachment 249926
As you can see, it agrees with your Python. Much simpler in Maple though.

Thanks for the reply. I've included a screenshot of the problem so that there's no ambiguity. The analytical/graphical approach certainly has its advantages.

prob.png
 
Just so you know, I didn't look up the theory of 2D Gaussian quadrature, so I am just assuming you know your formula is correct. What I mean is about your choices for the points and the weights. So, if your setup happens to be wrong, so is mine.
 
LCKurtz said:
Just so you know, I didn't look up the theory of 2D Gaussian quadrature, so I am just assuming you know your formula is correct. What I mean is about your choices for the points and the weights. So, if your setup happens to be wrong, so is mine.

As far as I understand it, I'm using the correct weights and points. After that, it's a fairly straightforward algorithm. Someone with a bit of expertise might notice a subtle detail that I'm overlooking.
 
First, I tried to show that ##f_n## converges uniformly on ##[0,2\pi]##, which is true since ##f_n \rightarrow 0## for ##n \rightarrow \infty## and ##\sigma_n=\mathrm{sup}\left| \frac{\sin\left(\frac{n^2}{n+\frac 15}x\right)}{n^{x^2-3x+3}} \right| \leq \frac{1}{|n^{x^2-3x+3}|} \leq \frac{1}{n^{\frac 34}}\rightarrow 0##. I can't use neither Leibnitz's test nor Abel's test. For Dirichlet's test I would need to show, that ##\sin\left(\frac{n^2}{n+\frac 15}x \right)## has partialy bounded sums...