Is My 2D Gaussian Quadrature Algorithm Accurate?

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.
 
There are two things I don't understand about this problem. First, when finding the nth root of a number, there should in theory be n solutions. However, the formula produces n+1 roots. Here is how. The first root is simply ##\left(r\right)^{\left(\frac{1}{n}\right)}##. Then you multiply this first root by n additional expressions given by the formula, as you go through k=0,1,...n-1. So you end up with n+1 roots, which cannot be correct. Let me illustrate what I mean. For this...
Back
Top