Wikipedia:Reference desk/Archives/Mathematics/2015 June 6

Welcome to the Wikipedia Mathematics Reference Desk Archives
The page you are currently viewing is an archive page. While you can leave answers for any questions shown below, please ask new questions on one of the current reference desk pages.

June 6

Gaussian quadrature

I am facing a task of integrating certain functions on 2-sphere. It is a 2-dimensional problem with one argument being angle φ (along the parallels of the 2-sphere) and the other angle θ along the meridians. I am trying to use Gauss-Legendre quadrature to integrate along the latter direction. The method requires an interval change and here I run into some conceptual problems. In Wikipedia article "Gaussian quadrature" the change of interval is presented as follows (I quote):

An integral over [a, b] must be changed into an integral over [−1, 1] before applying the Gaussian quadrature rule. This change of interval can be done in the following way:

$\int \limits _{a}^{b}{\mathfrak {f}}(x)~dx={\frac {b-a}{2}}\int \limits _{-1}^{1}{\mathfrak {f}}\left({\frac {b-a}{2}}x+{\frac {b+a}{2}}\right)~dx$ (1)

Applying the Gaussian quadrature rule then results in the following approximation:

$\int _{a}^{b}f(x)\,dx={\frac {b-a}{2}}\sum _{i=1}^{n}w_{i}f\left({\frac {b-a}{2}}x_{i}+{\frac {a+b}{2}}\right)$ (2)

I implement this formula in software. Abscissas and weights are given by many sources, e.g. here. For N = 4 they look like this:

X_1 = -0.861136  w_1 = 0.347854
X_2 = -0.339981  w_2 = 0.652145
X_3 =  0.339981  w_3 = 0.652145
X_4 =  0.861136  w_4 = 0.347854

Two points could be made:
$\bullet$ Weights w are always positive and vary in range: 0.0<w<1.0
$\bullet$ Argument x is always in the range -1.0<x<+1.0

What I don't understand is this.

I would expect the argument in parentheses in formula (2) on the right 

A=$\left({\frac {b-a}{2}}x_{i}+{\frac {a+b}{2}}\right)$

to vary in range -1.0<A<1.0. I don't think it is possible under any circumstances. Just assume the lower limit a equals 0. And it is also what my software tells me.

What am I missing here? Thanks for attention. --AboutFace 22 (talk) 18:28, 6 June 2015 (UTC)

The $x_{i}$ is in the range $[-1,1]$. You don't need A to be in this range as well, A is in the range $[a,b]$.

What happens basically: You take the original function; transform it to a different function g so that g in the range [-1, 1] corresponds to f in the range [a, b]; calculate the integral of g on [-1, 1]; then interpret the result in terms of the original function f.

The formula you quoted is already what results after the entire process was done for you. A value of $x_{i}$ for g maps back to the value of A for f. You're interested in the integral of f on [a, b], so of course you will have to evaluate f at points in the interval [a, b], not on the interval [-1, 1] (which could be completely meaningless for the original f). -- Meni Rosenfeld (talk) 19:32, 6 June 2015 (UTC)

(After edit conflict.) Let $A(t)=(b-a)t/2+(a+b)/2$. Then $A(-1)=a$ and $A(1)=b$. The transformation $A(t)$ goes from the interval $[-1,1]$ to the interval $[a,b]$. Obviously, this is what one wants if you are computing values of $f(A(t))$, since the domain of the function f is $[a,b]$.

However, this seems to be missing a larger issue, whether your numerical methodology is appropriate for integrating functions on spheres. Since the sphere is a symmetric space, I would think a better approach would be to exploit this. For example, for fewer than 30 sample points, simple averaging over the mid-edges of an icosahedron would be better than any adaptive approach to the problem (unless the function of interest is assumed to have some extra symmetries). If you absolutely insist on an adaptive method for the integral over the parallels, just a brute force "left-endpoint" approximation (equally spaced) will give better results than Gaussian quadrature because the circle is rotationally symmetric! Sławomir Biały (talk) 20:28, 6 June 2015 (UTC)

Thank you Meni & Slawomir. I will need more time to digest it all but now I want to make a few comments. First I don't understand what you, Slawomir, mean by "adaptive approach." I sample with 128 or 256 points. "the integral over the parallels" is a Fourier Transform and I have no problems there. I use this source. The Fourier Integral is given there by formula (9), and the Legendre transform by (10). And even there I have another question. The integrand in (10) has sin(θ) factor which disappears in the discrete expression (11). Why?

One thing of importance is that I don't have to worry about the expression (2) above. It has been settled. Many thinks, --AboutFace 22 (talk) 23:58, 6 June 2015 (UTC)