Unclear choice of polynomial in Gaussian Quadrature method

integrationnumerical methodsquadrature

Recently, I have been reading a paper regarding the evaluation of an integral using the Gaussian quadrature method but I cannot understand the rationality behind the weight and orthogonal polynomials the author has chosen.

$\begin{gathered}
\int\limits_{{s_1}}^{{s_2}} {f\left( x \right)dx} \approx \frac{{\left( {{s_2} – {s_1}} \right)}}{2}\sum\limits_{s = 1}^{\operatorname{Max} \,\,sample} {\frac{{\pi \sqrt {1 – {v^2}} }}{{\operatorname{Max} \,\,sample}}f\left( u \right)} \hfill \\
v = \cos \left( {\frac{{\left( {2s – 1} \right)\pi }}{{2 \times \operatorname{Max} \,\,sample}}} \right) \hfill \\
u = \frac{{\left( {{s_2} – {s_1}} \right)}}{2}v + \frac{{\left( {{s_2} + {s_1}} \right)}}{2} \hfill \\
\end{gathered} $

I understand that $u = \frac{{\left( {{s_2} – {s_1}} \right)}}{2}v + \frac{{\left( {{s_2} + {s_1}} \right)}}{2}$ is some smart change of variable to move the integration region to $\left[ { – 1,1} \right]$ but I do not know why a trigonometric function such as cos is chosen ?

From my domain knowledge, in this application the function $f\left( x \right)$ is something that look like a Gaussian bell-shaped. Does this knowledge justify the author's choice ?

Please help me understand this !

Thank you

Best Answer

That's Gauss-Chebyshev quadrature and will be the most accurate if $f(x)$ has branch points like $(x-s_1)^{n-1/2}$ and $(s_2-x)^{m-1/2}$ for nonnegative integers $n$ and $m$ at the endpoints. If $f(x)$ is regular at one or both endpoints, not as great but still can be pretty accurate. For regular $f(x)$ Gauss-Legendre quadrature would be more accurate although it requires more effort to set up. Clenshaw-Curtis quadrature can be easier to set up than Gauss-Legendre in that it uses the straightforward Gauss-Chebyshev quadrature to expand the function $f(x)$ in Chebyshev polynomials and then uses their known integrals. It would take more information to decide whether the author is using an optimal quadrature strategy or something quick and dirty.

EDIT: The issue is that the error in Gaussian quadrature is equal to $$\left(\frac{s_2-s_1}L\right)^{2n+1}C_nf^{(2n)}(\xi)$$ where $n$ is the number of data points, $L$ is the length of the model interval, $2$ for Gauss-Chebyshev quadrature, $s_1<\xi<s_2$, and $C_n$ is a constant depending on the quadrature algorithm and the number of data points. If $f(x)$ looks like $\sqrt{x-s_1}$ or $\sqrt{s_2-x}$ at one of its endpoints then the derivatives of $f(x)$ are unbounded so the error could be anything. Gauss-Chebyshev uses a weight function $w(x)=(1-x^2)^{-1/2}$ so this creates an effective $g(x)=\sqrt{(x-s_1)(s_2-x)}f(x)$ and this is the function whose derivatives should all remain finite at the endpoints if we want to impose a bound on the error of the quadrature as given.

As I mentioned, Gauss-Chebyshev quadrature may even give good results if $g(x)$ has unbounded derivatives but other methods may yield better results with fewer data points.

Related Solutions

Modified the Gauss–Laguerre quadrature in the case of $\int\limits_{x=0}^{+\infty}{{e^{-ax}}f(x)}dx$

The Gauss-Laguerre formula is for integrals of the form $\int_0^{+\infty}e^{-x} f(x) dx$. In your case, you just need to perform a change of variable $y = ax$ $$ \int_0^{+\infty} e^{-ax} f(x) dx = \frac 1a \int_0^{+\infty} e^{-y} \underbrace{f(y/a)}_{= \tilde f(y)} \,dy, $$ and apply the usual rule for $\tilde f$, i.e. $$ \int_0^{+\infty}e^{-ax} f(x) dx \approx \frac 1a \sum_{i=1}^n w_i f\left(\frac{x_i}{a}\right). $$
The order of the points is irrelevant... as long as you assign the correct weights.

Gaussian Quadrature strategy for $\int\limits_0^s {f\left( x \right){e^{ – ax}}} dx$

In principle you can use the Gram-Schmidt procedure to derive the orthogonal polynomials - which will result in modified Laguerre polynomials, however these quickly become unwieldy $$ p_0(x)=1\\ p_1(x) = x - 1+\frac{s e^{-s}}{1-e^{-s}}\\ p_2(x) = x^2 - Ap_1(x) - B $$ where $$ A = \frac{\frac{s^{3} e^{s}}{e^{\left(2 \, s\right)} - e^{s}} + \frac{2 \, s^{2} e^{s}}{e^{\left(2 \, s\right)} - e^{s}} + \frac{4 \, s e^{s}}{e^{\left(2 \, s\right)} - e^{s}} - \frac{2 \, s}{e^{\left(2 \, s\right)} - e^{s}} - \frac{2 \, s}{e^{s} - 1} + \frac{4 \, e^{s}}{e^{\left(2 \, s\right)} - e^{s}} - \frac{4 \, e^{s}}{e^{s} - 1} - \frac{4}{e^{\left(2 \, s\right)} - e^{s}} + \frac{4}{e^{s} - 1}}{\frac{s^{2} e^{\left(2 \, s\right)}}{e^{\left(3 \, s\right)} - 2 \, e^{\left(2 \, s\right)} + e^{s}} - \frac{s^{2}}{e^{\left(2 \, s\right)} - 2 \, e^{s} + 1} + \frac{e^{\left(2 \, s\right)}}{e^{\left(3 \, s\right)} - 2 \, e^{\left(2 \, s\right)} + e^{s}} - \frac{e^{\left(2 \, s\right)}}{e^{\left(2 \, s\right)} - 2 \, e^{s} + 1} - \frac{2 \, e^{s}}{e^{\left(3 \, s\right)} - 2 \, e^{\left(2 \, s\right)} + e^{s}} + \frac{2 \, e^{s}}{e^{\left(2 \, s\right)} - 2 \, e^{s} + 1} + \frac{1}{e^{\left(3 \, s\right)} - 2 \, e^{\left(2 \, s\right)} + e^{s}} - \frac{1}{e^{\left(2 \, s\right)} - 2 \, e^{s} + 1}}\\ B = \frac{s^{2} e^{-s} + 2s e^{-s}}{e^{-s} - 1}+2 $$

By fixing $s$ this process can be done semi-numerically with a computer algebra system like Sagemath or Mathematica. For instance, for $s=0.4$ the polynomials are $$ p_0(x)=1\\ p_1(x)= x-0.186702087312105\\ p_2(x)=x^2 - 0.389329291486829x + 0.0246035818232814 $$ which has roots $x_i\in [0.0793791632266308,0.309950128260198]$.

The next step is to solve the tower equations to determine the weights $w_i$ $$ \sum_{i=0}^2 w_ip_n(x_i)=\delta_{n,0}0.3296799539643607. $$ The number multiplying the delta is the norm of $p_0(x)$. This leads to the matrix equation $$ \left(\begin{array}{rr} 1 & 1 \\ -0.107322924085475 & 0.123248040948093 \end{array}\right)\left(\begin{array}{c} w_0\\ w_1 \end{array}\right) =\left(\begin{array}{c} 0.3296799539643607 \\ 0 \end{array}\right) $$ which has solution $w_i\in \left[0.1762251741456236, 0.153454779818737\right]$.

I've checked this by writing a small routine in Sagemath

def gqexp(f):
   xi = vector([0.0793791632266308,0.309950128260198])
   wi = vector([0.1762251741456236,0.153454779818737])
   fncall = vector([f(x) for x in xi])
   return expand(wi.dot_product(fncall))

Best Answer

Related Solutions

Modified the Gauss–Laguerre quadrature in the case of $\int\limits_{x=0}^{+\infty}{{e^{-ax}}f(x)}dx$

Gaussian Quadrature strategy for $\int\limits_0^s {f\left( x \right){e^{ – ax}}} dx$

Related Question