Editing Gaussian quadrature (section)

=== Fundamental theorem ===
Let {{mvar|p<sub>n</sub>}} be a nontrivial polynomial of degree {{mvar|n}} such that
<math display="block">\int_a^b \omega(x) \, x^k p_n(x) \, dx = 0, \quad \text{for all } k = 0, 1, \ldots, n - 1.</math>

Note that this will be true for all the orthogonal polynomials above, because each {{mvar|p<sub>n</sub>}} is constructed to be orthogonal to the other polynomials {{mvar|p<sub>j</sub>}} for {{math|''j''<''n''}}, and {{math|''x''<sup>''k''</sup>}} is in the span of that set.

If we pick the {{mvar|n}} nodes {{mvar|x<sub>i</sub>}} to be the zeros of {{mvar|p<sub>n</sub>}}, then there exist {{mvar|n}} weights {{mvar|w<sub>i</sub>}} which make the Gaussian quadrature computed integral exact for all polynomials {{math|''h''(''x'')}} of degree {{math|2''n'' − 1}} or less. Furthermore, all these nodes {{mvar|x<sub>i</sub>}} will lie in the open interval {{math|(''a'', ''b'')}}.<ref>{{harvnb|Stoer|Bulirsch|2002|pp=172–175}}</ref>

To prove the first part of this claim, let {{math|''h''(''x'')}} be any polynomial of degree {{math|2''n'' − 1}} or less. Divide it by the orthogonal polynomial {{mvar|p<sub>n</sub>}} to get
<math display="block"> h(x) = p_n(x) \, q(x) + r(x). </math>
where {{math|''q''(''x'')}} is the quotient, of degree {{math|''n'' − 1}} or less (because the sum of its degree and that of the divisor {{mvar|p<sub>n</sub>}} must equal that of the dividend), and {{math|''r''(''x'')}} is the remainder, also of degree {{math|''n'' − 1}} or less (because the degree of the remainder is always less than that of the divisor). Since {{mvar|p<sub>n</sub>}} is by assumption orthogonal to all monomials of degree less than {{mvar|n}}, it must be orthogonal to the quotient {{math|''q''(''x'')}}. Therefore
<math display="block"> \int_a^b \omega(x)\,h(x)\,dx = \int_a^b \omega(x)\,\big( \, p_n(x) q(x) + r(x) \, \big)\,dx = \int_a^b \omega(x)\,r(x)\,dx. </math>

Since the remainder {{math|''r''(''x'')}} is of degree {{math|''n'' − 1}} or less, we can interpolate it exactly using {{mvar|n}} interpolation points with [[Lagrange polynomials]] {{math|''l''<sub>''i''</sub>(''x'')}}, where
<math display="block"> l_i(x) = \prod _{j \ne i} \frac{x-x_j}{x_i-x_j}. </math>

We have
<math display="block"> r(x) = \sum_{i=1}^n l_i(x) \, r(x_i). </math>

Then its integral will equal
<math display="block"> \int_a^b \omega(x)\,r(x)\,dx = \int_a^b \omega(x) \, \sum_{i=1}^n l_i(x) \, r(x_i) \, dx = \sum_{i=1}^n \, r(x_i) \, \int_a^b \omega(x) \, l_i(x) \, dx = \sum_{i=1}^n \, r(x_i) \, w_i, </math>

where {{math|''w''<sub>''i''</sub>}}, the weight associated with the node {{math|''x''<sub>''i''</sub>}}, is defined to equal the weighted integral of {{math|''l''<sub>''i''</sub>(''x'')}} (see below for other formulas for the weights). But all the {{mvar|x<sub>i</sub>}} are roots of {{mvar|p<sub>n</sub>}}, so the division formula above tells us that
<math display="block"> h(x_i) = p_n(x_i) \, q(x_i) + r(x_i) = r(x_i), </math>
for all {{mvar|i}}. Thus we finally have
<math display="block"> \int_a^b \omega(x)\,h(x)\,dx = \int_a^b \omega(x) \, r(x) \, dx = \sum_{i=1}^n w_i \, r(x_i) = \sum_{i=1}^n w_i \, h(x_i). </math>

This proves that for any polynomial {{math|''h''(''x'')}} of degree {{math|2''n'' − 1}} or less, its integral is given exactly by the Gaussian quadrature sum.

To prove the second part of the claim, consider the factored form of the polynomial {{math|''p''<sub>''n''</sub>}}. Any complex conjugate roots will yield a quadratic factor that is either strictly positive or strictly negative over the entire real line. Any factors for roots outside the interval from {{mvar|a}} to {{mvar|b}} will not change sign over that interval. Finally, for factors corresponding to roots {{mvar|x<sub>i</sub>}} inside the interval from {{mvar|a}} to {{mvar|b}} that are of odd multiplicity, multiply {{math|''p''<sub>''n''</sub>}} by one more factor to make a new polynomial
<math display="block"> p_n(x) \, \prod_i (x - x_i). </math>

This polynomial cannot change sign over the interval from {{mvar|a}} to {{mvar|b}} because all its roots there are now of even multiplicity. So the integral
<math display="block"> \int_a^b p_n(x) \, \left( \prod_i (x - x_i) \right) \, \omega(x) \, dx \ne 0, </math>
since the weight function {{math|''ω''(''x'')}} is always non-negative. But {{math|''p''<sub>''n''</sub>}} is orthogonal to all polynomials of degree {{math|''n'' − 1}} or less, so the degree of the product
<math display="block"> \prod_i (x - x_i) </math>
must be at least {{mvar|n}}. Therefore {{math|''p''<sub>''n''</sub>}} has {{mvar|n}} distinct roots, all real, in the interval from {{mvar|a}} to {{mvar|b}}.

==== General formula for the weights ====
The weights can be expressed as

{{NumBlk|:|<math>w_{i} = \frac{a_{n}}{a_{n-1}} \frac{\int_{a}^{b} \omega(x) p_{n-1}(x)^2 dx}{p'_{n}(x_{i}) p_{n-1}(x_{i})}</math>|{{EquationRef|1}}}}

where <math>a_{k}</math> is the coefficient of <math>x^{k}</math> in <math>p_{k}(x)</math>. To prove this, note that using [[Lagrange interpolation]] one can express {{math|''r''(''x'')}} in terms of <math>r(x_{i})</math> as
<math display="block">r(x) = \sum_{i=1}^{n} r(x_{i}) \prod_{\begin{smallmatrix} 1 \leq j \leq n \\ j \neq i \end{smallmatrix}}\frac{x-x_{j}}{x_{i}-x_{j}}</math>
because {{math|''r''(''x'')}} has degree less than {{mvar|n}} and is thus fixed by the values it attains at {{mvar|n}} different points. Multiplying both sides by {{math|''ω''(''x'')}} and integrating from {{mvar|a}} to {{mvar|b}} yields
<math display="block">\int_{a}^{b}\omega(x)r(x)dx = \sum_{i=1}^{n} r(x_{i}) \int_{a}^{b}\omega(x)\prod_{\begin{smallmatrix} 1 \leq j \leq n \\ j \neq i \end{smallmatrix}} \frac{x-x_{j}}{x_{i}-x_{j}}dx</math>

The weights {{mvar|w<sub>i</sub>}} are thus given by
<math display="block">w_{i} = \int_{a}^{b}\omega(x)\prod_{\begin{smallmatrix}1\leq j\leq n\\j\neq i\end{smallmatrix}}\frac{x-x_{j}}{x_{i}-x_{j}}dx</math>

This integral expression for <math>w_{i}</math> can be expressed in terms of the orthogonal polynomials <math>p_{n}(x)</math> and <math>p_{n-1}(x)</math> as follows.

We can write
<math display="block"> \prod_{\begin{smallmatrix} 1 \leq j \leq n \\ j \neq i \end{smallmatrix}} \left(x-x_{j}\right) = \frac{\prod_{1\leq j\leq n} \left(x - x_{j}\right)}{x-x_{i}} = \frac{p_{n}(x)}{a_{n}\left(x-x_{i}\right)}</math>

where <math>a_{n}</math> is the coefficient of <math>x^n</math> in <math>p_{n}(x)</math>. Taking the limit of {{mvar|x}} to <math>x_{i}</math> yields using L'Hôpital's rule
<math display="block"> \prod_{\begin{smallmatrix} 1 \leq j \leq n \\ j \neq i \end{smallmatrix}} \left(x_{i}-x_{j}\right) = \frac{p'_{n}(x_{i})}{a_{n}}</math>

We can thus write the integral expression for the weights as

{{NumBlk|:|<math>w_{i} = \frac{1}{p'_{n}(x_{i})}\int_{a}^{b}\omega(x)\frac{p_{n}(x)}{x-x_{i}}dx</math>|{{EquationRef|2}}}}

In the integrand, writing
<math display="block">\frac{1}{x-x_i} = \frac{1 - \left(\frac{x}{x_i}\right)^{k}}{x - x_i} + \left(\frac{x}{x_i}\right)^{k} \frac{1}{x - x_i}</math>

yields
<math display="block">\int_a^b\omega(x)\frac{x^kp_n(x)}{x-x_i}dx = x_i^k \int_{a}^{b}\omega(x)\frac{p_n(x)}{x-x_i}dx</math>

provided <math>k \leq n</math>, because
<math display="block">\frac{1-\left(\frac{x}{x_{i}}\right)^{k}}{x-x_{i}}</math>
is a polynomial of degree {{math|''k'' − 1}} which is then orthogonal to <math>p_{n}(x)</math>. So, if {{math|''q''(''x'')}} is a polynomial of at most nth degree we have
<math display="block">\int_{a}^{b}\omega(x)\frac{p_{n}(x)}{x-x_{i}} dx = \frac{1}{q(x_{i})} \int_{a}^{b} \omega(x)\frac{q(x) p_n(x)}{x-x_{i}}dx </math>

We can evaluate the integral on the right hand side for <math>q(x) = p_{n-1}(x)</math> as follows. Because <math>\frac{p_{n}(x)}{x-x_{i}}</math> is a polynomial of degree {{math|''n'' − 1}}, we have
<math display="block">\frac{p_{n}(x)}{x-x_{i}} = a_{n}x^{n-1} + s(x)</math>
where {{math|''s''(''x'')}} is a polynomial of degree <math>n - 2</math>. Since {{math|''s''(''x'')}} is orthogonal to <math>p_{n-1}(x)</math> we have
<math display="block">\int_{a}^{b}\omega(x)\frac{p_{n}(x)}{x-x_{i}}dx=\frac{a_{n}}{p_{n-1}(x_{i})} \int_{a}^{b}\omega(x)p_{n-1}(x)x^{n-1}dx </math>

We can then write
<math display="block">x^{n-1} = \left(x^{n-1} - \frac{p_{n-1}(x)}{a_{n-1}}\right) + \frac{p_{n-1}(x)}{a_{n-1}}</math>

The term in the brackets is a polynomial of degree <math>n - 2</math>, which is therefore orthogonal to <math>p_{n-1}(x)</math>. The integral can thus be written as
<math display="block">\int_{a}^{b}\omega(x)\frac{p_{n}(x)}{x-x_{i}}dx = \frac{a_{n}}{a_{n-1} p_{n-1}(x_{i})} \int_{a}^{b}\omega(x) p_{n-1}(x)^{2} dx </math>

According to equation ({{EquationNote|2}}), the weights are obtained by dividing this by <math>p'_{n}(x_{i})</math> and that yields the expression in equation ({{EquationNote|1}}).

<math>w_{i}</math> can also be expressed in terms of the orthogonal polynomials <math>p_{n}(x)</math> and now <math>p_{n+1}(x)</math>. In the 3-term recurrence relation <math>p_{n+1}(x_{i}) = (a) p_{n}(x_{i}) + (b) p_{n-1}(x_{i})</math> the term with <math>p_{n}(x_{i})</math> vanishes, so <math>p_{n-1}(x_{i})</math> in Eq. (1) can be replaced by <math display="inline">\frac{1}{b} p_{n+1} \left(x_i\right)</math>.

====Proof that the weights are positive====
Consider the following polynomial of degree <math>2n - 2</math>
<math display="block">f(x) = \prod_{\begin{smallmatrix} 1 \leq j \leq n \\ j \neq i \end{smallmatrix}}\frac{\left(x - x_j\right)^2}{\left(x_i - x_j\right)^2}</math>
where, as above, the {{mvar|x<sub>j</sub>}} are the roots of the polynomial <math>p_{n}(x)</math>. 
Clearly <math>f(x_j) = \delta_{ij}</math>. Since the degree of <math>f(x)</math> is less than <math>2n - 1</math>, the Gaussian quadrature formula involving the weights and nodes obtained from <math>p_{n}(x)</math> applies. Since <math>f(x_{j}) = 0</math> for {{mvar|j}} not equal to {{mvar|i}}, we have
<math display="block">\int_{a}^{b}\omega(x)f(x)dx=\sum_{j=1}^{n}w_{j}f(x_{j}) = \sum_{j=1}^{n} \delta_{ij} w_j = w_{i} > 0.</math>

Since both <math>\omega(x)</math> and <math>f(x)</math> are non-negative functions, it follows that <math>w_{i} > 0</math>.