Orthogonal polynomials

Orthogonal polynomials is specific set of orthogonal functions, used, in particular, for construction of formulas of numerical integration. Many widely used formulas for the numerical integration appear at consideration of specific orthogonal polynomials. In particular, Gauss-Hermite quadrature, Gauss-Legendre quadrature and Gauss-Legendre quadrature appear to be special case of such an application.

Aiming this goal, the description is loaded to TORI. However, any beautiful construction is expected to have other applications too.

For the set of orthogonal polynomials \(P_m\), \(m=0,1,2, ..\) at the interval \((a,b)\) with weight function \(\rho\) is assumed relation of orthogonality

\(\displaystyle \int_a^b \rho(x)\, P_m(x)\, P_n(x) \, \mathrm d x =0~\) for \(~n\ne m\).

Usially, function \(\rho\) is supposed to he real and positive at the interval \((a,b)\). As for \(a\) and \(b\), they are supposed to be real number, \(a<b\); some of them may be also infinite, \(a=-\infty\), or \(b=\infty\), or both.

Coefficients

The polynomial is determined with its coefficients \(s\),

\(\displaystyle P_m(x)=\sum_{n=0}^{n} s_{m,n} x^n\)

that need to be chosen in such a way, that the orthogonality above holds. Usually it is supposed, that the highest coefficient is just unity, \(s_{m,m}=1\). Then, the recurrent relation takes place,

\(\displaystyle P_{m+1}(x)= (x\!-\!u_m) P_m(x) - v_m P_{m-1}(x)\)

where

\(\displaystyle u_{m}=\frac {\int_a^b \rho(x)\, x\, P_m(x)^2\, \mathrm d x} {\int_a^b \rho(x)\, P_m(x)^2\, \mathrm d x}~\) , \(~\displaystyle v_{m}=\frac {\int_a^b \rho(x)\, x\, P_m(x)\, P_{m-1}(x)\, \mathrm d x} {\int_a^b \rho(x)\, P_{m-1}(x)^2\, \mathrm d x}\)

Denominators in the last expressions can be interpreted as squares of norm of the corresponding polynomials; and it is assumed, that \(v_0\!=\!0\).

Example

As an example, the weight function \(\rho(x)=\exp(-x^2)\) is considered here, with \(a\!=\! 0\) and \(b\!=\!\infty\). Namely for this case, the numerical integration is required for evaluation of the contour integral for function naga. Denote them as Halfmite polynomials, as the interval \((a,b)=(0,\infty)\) is a half of interval \((a,b)=(-\infty,\infty)\) that leads to the Hermite polynomials.

In Mathematica, construction of the first 5 Halfmite polynomials can be implemented as follows:

R[x_] = Exp[-x^2]
P0[x_] = 1
u0 = Integrate[R[x] x P0[x]^2, {x,0,Infinity}]/Integrate[R[x] P0[x]^2, {x,0,Infinity}]
P1[x_] = x - u0
u1 = Integrate[R[x] x P1[x]^2, {x,0,Infinity}]/Integrate[R[x] P1[x]^2, {x,0,Infinity}]
v1 = Integrate[R[x] x P1[x] P0[x], {x,0,Infinity}]/Integrate[R[x] P0[x]^2, {x,0,Infinity}]
P2[x_] = (x - u1) P1[x] - v1 P0[x]
u2 = Integrate[R[x] x P2[x]^2, {x,0,Infinity}]/Integrate[R[x] P2[x]^2, {x,0,Infinity}]
v2 = Integrate[R[x] x P2[x] P1[x], {x,0,Infinity}]/Integrate[R[x] P1[x]^2, {x,0,Infinity}]
P3[x_] = (x - u2) P2[x] - v2 P1[x]
u3 = Integrate[R[x] x P3[x]^2, {x,0,Infinity}]/Integrate[R[x] P3[x]^2, {x,0,Infinity}]
v3 = Integrate[R[x] x P3[x] P2[x], {x,0,Infinity}]/Integrate[R[x] P2[x]^2, {x,0,Infinity}]
P4[x_] = (x - u3) P3[x] - v3 P2[x]
u4 = Integrate[R[x] x P4[x]^2, {x,0,Infinity}]/Integrate[R[x] P4[x]^2, {x,0,Infinity}]
v4 = Integrate[R[x] x P4[x] P3[x], {x,0,Infinity}]/Integrate[R[x] P3[x]^2, {x,0,Infinity}]
P5[x_] = (x - u4) P4[x] - v4 P3[x]
r[x_] = Exp[-x^2/2];
Plot[{P1[x] r[x], P2[x] r[x], P3[x] r[x], P4[x] r[x]}, {x,0,4}, PlotRange->All, AspectRatio->Automatic]

Zeros of the 5th halfmite polumomial are solutions of equation \(P_5(x)\!=\!0\). Their approximations are:

0.100242151968215596318612501691,
0.482813966046200700004307982394,
1.06094982152571708116749191432,
1.77972941852026130183228102279,
2.66976035608765641367980178954

These numbers can be considered as nodes of the Gauss-Halfmite quadrature of 5th order. The corresponding scheme is supposed to integrate exactly the approximation of integrand with product of polynomial of 9th order to the Gaussian exponential.

Coefficients of the Halfmite polynomials appear as rational combinations of \(\pi\) and \(\sqrt{\pi}\); and for the 5th Haldmite polynomial still can be calculated exactly in real time;

\(P_5(x)= x^5 +\) \(\! \frac{\left(-4096+3000 \pi -540 \pi^2\right) x^4}{4 \sqrt{\pi } \left(656-435 \pi+72 \pi ^2\right)} +\) \(\! \frac{\left(-14080 \sqrt{\pi}+9006 \pi^{3/2}-1440 \pi^{5/2}\right) x^3}{4\sqrt{\pi} \left(656-435 \pi +72 \pi^2\right)} +\) \(\! \frac{\left(16384-11436 \pi +1980 \pi^2\right) x^2}{4 \sqrt{\pi } \left(656-435 \pi+72 \pi^2\right)} +\) \(\! \frac{\left(11904 \sqrt{\pi}-7182 \pi ^{3/2}+1080 \pi ^{5/2}\right) x}{4\sqrt{\pi } \left(656-435 \pi +72 \pi ^2\right)}+\frac{-8192+5124 \pi -801 \pi ^2}{4 \sqrt{\pi } \left(656-435 \pi +72 \pi ^2\right)} \)

For the highest Halfmite polynomials, the precise analytical calculus should be implemented.

Quadrature formula

For the \(N\)-point quadrature formula, based on the orthogonal polynomials, nodes \(x_n\) are chosen at zeros of the Nth polynomial,

\(P_{N}(x_n)=0~\)

The weights are

\(w_n=\displaystyle \frac{1}{P_{N-1}(x_n)\, P_N^{~\prime}(x_n)} \int_a^b{\rho(x)\, P_{N-1}(x)^2\, \mathrm d x}\)

In particular, for the example above,

References

http://www.mmcs.sfedu.ru/_old/docmanupload/doc_download/337---q-q-- С. Н. Овчинникова. Численное интегрирование. Методическое пособие ФЕДЕРАЛЬНОЕ АГЕНТСТВО ПО ОБРАЗОВАНИЮ Федеральное государственное образовательное учреждение высшего профессионального образования «ЮЖНЫЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ», 2007 год. (In Russian. Warning! The deduction looks fine and detailed there, but there are many misprints in formulas, see p.18. They should be revised before the implementation.)

https://en.wikipedia.org/wiki/Orthogonal_polynomials

Keywords

Gauss-Hermite quadrature, Gauss-Legendre quadrature, Gauss-Legendre quadrature, Numerical integration,