# Numerical integration using Gaussian quadrature, but with data at some nodes not available?

In order to calculate the integration of a function $f(x)$ from $[-1,1]$ using Gaussian quadrature, the values of $f(x)$ at the Gaussian nodes (in my case it's for a 6-node calculation, $f(x)$ at $x = -0.9325, -0.6612, -0.2386, 0.2386, 0.6612, 0.9325$) need to be known. But in my setup, I cannot obtain data points between $[-1, -0.85]$ or $[0.85,1]$, while the points between $(-0.85,0.85)$ can be measured quite densely(tens of data points can be obtained). My question is, since the data at two Gaussian nodes are not available, how could I perform the numerical integration (any quadrature would do) to obtain the desired results based on the available data please? Thanks a lot in advance for your help!

First of all Thanks for the replies @uranix @Count Iblis. My question actually originated from calculating spherical harmonic expansion coefficients up to 4th degree using 6th order Gaussian-Legendre quadrature, where theoretically if I could detect the values with polar angle (the angle between any unit vector and the z axis, hence ranging between 0 to 180 deg) equal to 21.2, 46.7, 76.2 deg etc (6 polar angles *12 equal-space azimuthal angles), I would be able to calculate the expansion coefficients. However, it turned out that the experiment can only detect the values with polar angle between ~35 to 145 deg, what should I do if I still want to calculate the SH coefficients please? What quadrature should I choose and what directions should I measure please? Many thanks again.

• You can do a transform that maps the points 1 and -1 to plus and minus infinity respectively (or you split the integral into a part from -1 to 0 and from 10 1 and treat them separately to map the -1 and 1 endpoints to infinity in the same way) . Put an extra parameter in the mapping giving you control of where 0.85 and -0.85 end up. Then you perform Gauss–Hermite quadrature, or in case of splitting up the original integral, you can do Gauss–Laguerre quadrature. – Count Iblis Nov 9 '15 at 20:22
• Do you need to integrate $f(x)$ over $[-1, 1]$ or over $[-0.85, 0.85]$? – uranix Nov 9 '15 at 20:45
• Thanks very much for the replies, would you kindly have a look at my updated part of the question please. – Lampard Nov 10 '15 at 0:56

Let's focus on quadratures with $2n$ nodes that integrate polynomials (up to some degree) precisely.

$$\int_{-1}^1 f(x) dx \approx \sum_{i=1}^{2n} w_i f(x_i), \quad |x_i| < \Delta = 0.85.$$

Due to symmetry one can take $x_{i} = x_{2n+1 - i}$, so the quadrature becomes $$\int_{-1}^1 f(x) dx \approx \sum_{i=1}^{n} w_i \left[f(x_i) + f(-x_i)\right], \quad 0 \leq x_i < \Delta = 0.85.$$

Assume the quadrature is precise for polynomials of degree up to $2p + 1$, then $$\frac{1}{2}\int_{-1}^1 x^{2k} dx = \frac{1}{2k+1} = \sum_{i=1}^{n} w_i x_i^{2k}, \qquad k = 0, 1, \dots, p.$$ Introduce $\xi_i = \frac{x_i^2}{\Delta^2} \in [0, 1]$. Then we have a problem for $\xi_i, w_i$ $$\sum_{i=1}^n w_i \xi_i^k = \frac{\Delta^{-2k}}{(2k+1)}, \quad k = 0, 1, \dots, p$$

Let's investigate if this problem admits solutions with strictly positive $w_i$.

### Case when all $w_i > 0$

If this problem has a solution then $$\mu(\xi) = \sum_{i=1}^n w_i \theta(\xi - \xi_i)$$ is a solution Hausdorff moment problem $$\int_0^1 \xi^k d\mu(\xi) = m_k, k = 0, 1, \dots$$ with $$m_k = \begin{cases} \displaystyle \frac{\Delta^{-2k}}{(2k+1)}, &k = 0, 1, \dots, p\\ \displaystyle \int_0^1 \xi^k d\mu(\xi), &k > p \end{cases},$$ thus $m_k$ should be completely monotonic. $$\begin{array}{c|c} k & m_k\\\hline 0 & 1\\ 1 & 0.461361\\ 2 & 0.383137\\ 3 & 0.378781\\ 4 & 0.407761\\ \end{array}$$ It can easily be seen that for $p > 3$ the sequence is not even monotonic, so there is no solution with positive $w_i$ for $p > 3$. The solution (one of) for $p = 3$ is simply the 4-point Gaussian quadrature rule $$w_{1,2} = \frac{18 \pm \sqrt{30}}{36}\\ \xi_{1,2} = \Delta^2 \left(\frac{3}{7} \mp \frac{2}{7}\sqrt\frac{6}{5}\right).$$

### Case when some $w_i < 0$

Actually now one can simply solve the moment system $$\sum_{i=1}^n w_i \xi_i^k = \frac{\Delta^{-2k}}{(2k+1)}, \quad k = 0, 1, \dots, p.$$ Note that it has $p + 1$ equation and $2n$ unknowns, so one may expect a solution to exist for $2n \geq p + 1$. The case $2n = p + 1$ give simply the Gaussian quadrature rules, for which $\xi_i$ lie outside $[0, 1]$ when $p > 3$. When $2n > p + 1$ there is no unique solution, so one may try to minimize the truncation error, proportional to $$\delta_{p+1} = \sum_{i=1}^n w_i \xi_i^{p+1} - \frac{\Delta^{-2p-2}}{(2p+3)}.$$

This leads us to a nonlinear programming problem \begin{aligned} \operatorname*{minimize}_{w_i, \xi_i} &&&\delta_{p+1}\\ \operatorname*{subject to} &&&0 \leq \xi_1\\ &&&\xi_1 \leq \xi_2\\ &&&\xi_2 \leq \xi_3\\ &&&\vdots\\ &&&\xi_{n-1} \leq \xi_n\\ &&&\xi_n \leq 1\\ &&&\delta_{k} = 0, \quad k = 0, 1, \dots, p \end{aligned}

Also one may try to minimize the absolute values of the weights, since the previous approach may lead to big negative weights \begin{aligned} \operatorname*{minimize}_{w_i, \xi_i} &&&\sum_{i=1}^n |w_i|\\ \operatorname*{subject to} &&&0 \leq \xi_1\\ &&&\xi_1 \leq \xi_2\\ &&&\xi_2 \leq \xi_3\\ &&&\vdots\\ &&&\xi_{n-1} \leq \xi_n\\ &&&\xi_n \leq 1\\ &&&\delta_{k} = 0, \quad k = 0, 1, \dots, p \end{aligned}

Using GAMS/BARON I've found some quadratures of this kind (you may also use NEOS server to solve GAMS models)

### For $p = 4$

The problem seems to be infeasible for $n \leq 3$. The $n = 4$ quadrature is \begin{aligned} \xi_{1,7}&=1.0000000000000000000\\ \xi_{2,6}&=0.85880349620038677732\\ \xi_{3,5}&=0.52084678062073430276\\ \xi_{4}&=0\\ x_{1,7}&=\pm 0.85000000000000000000\\ x_{2,6}&=\pm 0.78770903638639276855\\ x_{3,5}&=\pm 0.61344258003376366028\\ x_{4}&=0\\ w_{1,7}&=0.83593508805857758387\\ w_{2,6}&=-0.88789837944004855297\\ w_{3,5}&=0.74485659488394778770\\ w_{4}&=0.30710669649752318141\\ \end{aligned}

### For $p = 5$

The problem seems to be infeasible for $n \leq 3$. The $n = 4$ quadrature is \begin{aligned} \xi_{1,8}&=1.0000000000000000000\\ \xi_{2,7}&=0.95880851739522516021\\ \xi_{3,6}&=0.86778546156594895262\\ \xi_{4,5}&=0.090906828432672730508\\ x_{1,8}&=\pm 0.85000000000000000000\\ x_{2,7}&=\pm 0.83230953005360342218\\ x_{3,6}&=\pm 0.79181752694759044356\\ x_{4,5}&=\pm 0.25628145376247194956\\ w_{1,8}&=5.8486928393996124204\\ w_{2,7}&=-8.6490143669092829039\\ w_{3,6}&=3.2951628882441291627\\ w_{4,5}&=0.50515863926554132074\\ \end{aligned} I suggest not using this quadrature due to big negative weight $w_{2,7}$

The $n = 5$ quadrature is \begin{aligned} \xi_{1,10}&=1.0000000000000000000\\ \xi_{2,9}&=0.90991052436123162672\\ \xi_{3,8}&=0.67405321511028038136\\ \xi_{4,7}&=0.38251754788591463719\\ \xi_{5,6}&=0.14666023863496339182\\ x_{1,10}&=\pm 0.85000000000000000000\\ x_{2,9}&=\pm 0.81080845694343238659\\ x_{3,8}&=\pm 0.69785632326230130357\\ x_{4,7}&=\pm 0.52570802575914068441\\ x_{5,6}&=\pm 0.32551808308258552196\\ w_{1,10}&=1.8511750976308568442\\ w_{2,9}&=-2.6071695839739403098\\ w_{3,8}&=1.7504677140719977735\\ w_{4,7}&=-0.84051661731816264992\\ w_{5,6}&=0.84604338958924834195\\ \end{aligned}

Just in case if you would like to reproduce the results here is my Mathematica script.

## Quadrature for polynomials with weight

One of the possibilities is to build a special quadrature rule with weight $\omega(x) = \exp\left(-\frac{x^2}{2\sigma^2}\right)$.

The case $\sigma = \infty$ brings us to Gauss-Legendre quadrature. When $\sigma \to 0$ then the quadrature nodes group around $x = 0$.

To make a weighted quadrature, we need to build an orthogonal polynomial family with respect to the inner product $$\langle f,g\rangle = \int_{-1}^{1} f(x) g(x) \omega(x) dx.$$

Quick and dirty way is to make the Gramian matrix $$G_{ij} = \int_{-1}^{1} x^{i+j} \omega(x) dx\\$$ perform its Cholesky decomposition $$G = LL^\top\\ M = L^{-1}$$ and the polynomials are $$P_i(x) = \sum_{j=0}^{n} M_{ij} x^j$$ Note that the matrix $M$ actually depends on $\sigma$. The nodes of the quadrature are roots of $P_n(x)$ polynomial, so we want to find the smallest $\sigma$ for which $$P_n(0.85) = 0.$$ Thus $\pm 0.85$ will be the outermost nodes.

For $n = 6$ solving the equation numerically gives $\sigma = 0.431946$, this corresponds to the quadrature $$\int_{-1}^1 f(x) \omega(x) dx \approx \sum_{j=1}^6 w_j f(x_j)\\ x_{1,6} = \pm 0.8500000000\\ x_{2,5} = \pm 0.5142302766\\ x_{3,4} = \pm 0.1706206462\\ w_{1,6} = 0.006369511915\\ w_{2,5} = 0.0836555022 \\ w_{3,4} = 0.292371290.$$ Introducing $g(x) = f(x)\omega(x)$ $$\int_{-1}^1 g(x) dx \approx \sum_{j=1}^6 w_j \frac{g(x_j)}{\omega(x_j)} = \sum_{j=1}^6 \tilde w_j g(x_j)\\ \tilde w_{1,6} = 0.306100114\\ \tilde w_{2,5} = 0.345153691\\ \tilde w_{3,4} = 0.341741013.$$

Note that this quadrature is not quite good when $g(x)$ is a polynomial, but is exact when $\frac{g(x)}{\omega(x)}$ is a polynomial of 11-th degree.

Appendix. For $n = 10$ the same method gives $\sigma \approx 0.2565088$ and $$x_{1,10} = \pm 0.85000000\\ x_{2,9} = \pm 0.63553066\\ x_{3,8} = \pm 0.44273486\\ x_{4,7} = \pm 0.26175323\\ x_{5,6} = \pm 0.086652210\\ \tilde w_{1,10} = 0.22829400\\ \tilde w_{2,9} = 0.20160757\\ \tilde w_{3,8} = 0.18564405\\ \tilde w_{4,7} = 0.17726595\\ \tilde w_{5,6} = 0.17359644.$$

• Goodness me @uranix, so much useful information in the post, thank you so much! While I need some time to digest these slowly, would you please have a look at my updated question on spherical harmonics please？thanks again. – Lampard Nov 10 '15 at 0:58
• Hi @uranix I have tried the quadrature you generated, as you noted in your answer, it does not work that well when $g(x)$ is a polynomial, which means the spherical harmonic coefficients for $Y_{l0}$ are problematic as they correspond to the Legendre polynomials. How should I overcome this problem please? – Lampard Nov 11 '15 at 14:06
• @Lampard I'll try another approach and will update my answer – uranix Nov 11 '15 at 16:27