Gaussian Quadrature, Double Integral in python

Question

I am trying to solve an integral, such that

$$F_z = G\sigma z\int\int_{-L/2}^{L/2} \frac{dxdy} {(x^2 + y^2 + z^2)^{3/2}}$$

I tried to appraoch the problem in this way,

$$F(x,y) = \int_{-L/2}^{L/2} \frac{dx} {(x^2 + y^2 + z^2)^{3/2}}$$

and then

$$F_z = G\sigma z \int_{-L/2}^{L/2} F(x,y) dy$$

But in the coding part I learned gaussian quadrature for only x component, however in the function there are 2 components, so when I try to apply the gaussian quadrature I dont know what the y value should be. When I take the integral they should be treated as constants I guess ?

from numpy import linspace, cos, pi, tan, ones, copy

sigma = 100 #kg/m^2 sheet density
G = 6.674 * 10**-11 #m^3kg^-1s^-2
L = 10 # m
N = 100 # step size

def gaussxw(N):

    # Initial approximation to roots of the Legendre polynomial
    a = linspace(3,4*N-1,N)/(4*N+2)
    x = cos(pi*a+1/(8*N*N*tan(a)))

    # Find roots using Newton's method
    epsilon = 1e-15
    delta = 1.0
    while delta>epsilon:
        p0 = ones(N,float)
        p1 = copy(x)
        for k in range(1,N):
            p0,p1 = p1,((2*k+1)*x*p1-k*p0)/(k+1)
        dp = (N+1)*(p0-x*p1)/(1-x*x)
        dx = p1/dp
        x -= dx
        delta = max(abs(dx))

    # Calculate the weights
    w = 2*(N+1)*(N+1)/(N*N*(1-x*x)*dp*dp)

    return x,w

def f(x):
    return (x**2 + y**2 + z**2)**(-3/2) 

def gaussian(a, b):
    x, w = gaussxw(N)
    xp = 0.5 * (b - a) * x + 0.5 * (b + a)
    wp = 0.5 * (b - a) * w
    s = 0
    for k in range(N):
        s += wp[k] * f(xp[k])
    return s

z = 1
Fz = G * sigma * z * gaussian(-L/2, L/2)

The problem is in the above code what should be the y value ? (which thats basically what I am asking) or should I change the gaussxw to gaussyw and try that ? (but I am also not sure how to do that)

Answer 1

$\displaystyle\int_{-L/2}^{L/2}\frac{dx}{(x^2+y^2+z^2)^{3/2}}$ is not a function of $x$. If we call it $F(y)$, then when you compute $\displaystyle\int_{-L/2}^{L/2}F(y)\,dy$ using the Gaussian quadrature, as a subproblem you have to compute $F(y)$ at given $y$. As $F(y)$ is represented by an integral, you can compute it using the quadrature again. That is, if the quadrature you use for a single integral is $\int_{-L/2}^{L/2}f(x)\,dx\approx\sum_{i=1}^{n}w_i f(x_i)$, then for a double integral it is used like this: $$\int_{-L/2}^{L/2}\int_{-L/2}^{L/2}f(x,y)\,dx\,dy\approx\sum_{i=1}^{n}\sum_{j=1}^{n}w_i w_j f(x_i,x_j).$$

A side note: the integral has a closed form. Denoting $a=L/2$ and using polar coordinates, we have \begin{align} z\iint\limits_{[-a,a]^2}\frac{dx\,dy}{(x^2+y^2+z^2)^{3/2}}&=8z\int\limits_{0}^{\pi/4}\int\limits_{0}^{a/\!\cos\phi}\frac{r\,dr\,d\phi}{(r^2+z^2)^{3/2}}\\&=8z\int_0^{\pi/4}\left(\frac{1}{z}-\frac{\cos\phi}{\sqrt{a^2+z^2\cos^2\phi}}\right)d\phi\\\color{gray}{[w=z\sin\phi]}&=2\pi-8\int_0^{z/\sqrt{2}}\frac{dw}{\sqrt{a^2+z^2-w^2}}\\&=2\pi-8\arcsin\frac{z}{\sqrt{2(a^2+z^2)}}. \end{align}

Answer 2

Addendum: if you are just interested in asymptotics for small or large values of $L$ and $z$, the following inequalities are fairly simple but useful. The square $Q(L)$ centered at the origin with side length $L$ contains a circle with radius $\frac{L}{2}$ and its is contained in a circle with radius $\frac{L}{\sqrt{2}}$. It follows that $$\begin{eqnarray*} I(L,z) &=& \iint_{Q(L)}\frac{d\mu}{(x^2+y^2+z^2)^{3/2}}\leq \int_{0}^{L/\sqrt{2}}\frac{2\pi\rho}{(\rho^2+z^2)^{3/2}}\,d\rho\\&=&2\pi\left(\frac{1}{z}-\frac{1}{\sqrt{z^2+\frac{L^2}{2}}}\right)=\frac{\pi L^2}{z\sqrt{z^2+L^2/2}\left(z+\sqrt{z^2+L^2/2}\right)}\end{eqnarray*}$$ as well as $$\begin{eqnarray*} I(L,z) &=& \iint_{Q(L)}\frac{d\mu}{(x^2+y^2+z^2)^{3/2}}\geq \int_{0}^{L/2}\frac{2\pi\rho}{(\rho^2+z^2)^{3/2}}\,d\rho\\&=&2\pi\left(\frac{1}{z}-\frac{1}{\sqrt{z^2+L^2/4}}\right)=\frac{\tfrac{\pi}{2} L^2}{z\sqrt{z^2+L^2/4}\left(z+\sqrt{z^2+L^2/4}\right)}\end{eqnarray*}$$ and the actual value of $I(Q,z)$ is expected to be reasonably close to the integral of $\frac{1}{(x^2+y^2+z^2)^{3/2}}$ over a circle centered at the origin with the same area as $Q(L)$, $$ I(L,z)\approx \frac{2L^2}{z\sqrt{z^2+L^2/\pi}\left(z+\sqrt{z^2+L^2/\pi}\right)}. $$ Given the closed form for $I(L,z)$ provided by metamorphy, it is interesting to check that the last approximation is indeed a fairly good one.