1
$\begingroup$

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!

Update about the original question:

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Count Iblis Nov 9 '15 at 20:22
  • $\begingroup$ Do you need to integrate $f(x)$ over $[-1, 1]$ or over $[-0.85, 0.85]$? $\endgroup$ – uranix Nov 9 '15 at 20:45
  • $\begingroup$ Thanks very much for the replies, would you kindly have a look at my updated part of the question please. $\endgroup$ – Lampard Nov 10 '15 at 0:56
2
$\begingroup$

Quadrature for polynomials

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. $$

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Lampard Nov 10 '15 at 0:58
  • $\begingroup$ 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? $\endgroup$ – Lampard Nov 11 '15 at 14:06
  • $\begingroup$ @Lampard I'll try another approach and will update my answer $\endgroup$ – uranix Nov 11 '15 at 16:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.