MATLAB Examples

Orthogonal polynomials via the Lanczos process

Pedro Gonnet, 1st November 2011

(Chebfun example approx/OrthPolysLanczos.m)

```format short e ```

The example approx/OrthPolys [2] discusses orthogonal polynomials, defined by the condition

```    b
/
| w(x) P_j(x) P_k(x) dx = < P_j, P_k > = 0 for all j not equal k
/
a```

As mentioned there, Chebfun has commands LEGPOLY, CHEBPOLY, LAGPOLY, HERMPOLY for computing some standard cases.

Here, instead of Gram-Schmidt, we construct some of these polynomials via a three-term recurrence relation

`  gamma(k+1)*P_k+1(x) = (x - beta(k))*P_k(x) - gamma(k)*P_k-1(x)`

where the coefficients gamma(k) and beta(k) are either known in advance or computed from

`  gamma(k) = < x*P_k-1 , P_k >,  beta(k) = < x*P_k , P_k >.`

Given any positive weight function w(x), we can construct both the polynomials P_k(x) and the recurrence coefficients beta(k) and gamma(k) using the Lanczos algorithm [1].

Whereas the Gram-Schmidt process requires O(n^2) evaluations of the inner product to compute the n first polynomials, the Lanczos process requires only two such evaluations per polynomial. Furthermore, the three-term recurrence coefficients can be used to construct a Gaussian quadrature rule for the given weight function w(x).

We start by initializing the parameters of this set of polynomials, e.g. the variable x in the interval [-1,1], the weight function in that same interval as well as the highest-degree polynomial we wish to construct.

```x = chebfun( 'x' , [-1,1] ); w = exp(pi*x); N = 5; ```

We then initialize the recursion by constructing the first two polynomials and their recurrence coefficients explicitly.

```P = chebfun( 1./sqrt(sum(w)) , domain(x) ); v = x.*P; beta(1) = sum(w.*v.*P); v = v - beta(1)*P; gamma(1) = sqrt(sum( w.*v.^2 )); P(:,2) = v / gamma(1); ```

The remaining polynomials and coefficients can be computed iteratively using the Lanczos process.

```for k=2:N v = x.*P(:,k) - gamma(k-1)*P(:,k-1); beta(k) = sum(w.*v.*P(:,k)); v = v - beta(k)*P(:,k); gamma(k) = sqrt(sum( w.*v.^2 )); P(:,k+1) = v / gamma(k); end ```

We can now plot these polynomials

```figure plot(P) title('Orthogonal polynomials on [-1,1] wrt w = exp(pi*x)'); ```

and verify that they are indeed orthogonal with respect to w.

```W = repmat(w,1,N+1); I = P'*(W.*P); err = norm(I-eye(N+1)) ```
```err = 2.0425e-14 ```

Since, along with the polynomials, we have also computed the coefficients beta and gamma of the three-term recurrence relation, we can also construct Gaussian quadrature rules with respect to the weight function w using the Golub-Welsch algorithm [3]. For N points, this is

```J = diag(beta) + diag(gamma(1:N-1),-1) + diag(gamma(1:N-1),+1) [V,D] = eig( J ); xi = diag( D ); wk = (2*V(1,:).^2)'; [ xi , wk ] ```
```J = 6.8543e-01 3.0631e-01 0 0 0 3.0631e-01 1.5836e-01 4.9306e-01 0 0 0 4.9306e-01 -2.4896e-02 5.1638e-01 0 0 0 5.1638e-01 -1.7956e-02 5.0738e-01 0 0 0 5.0738e-01 -6.5923e-03 ans = -8.2685e-01 8.3551e-03 -2.9103e-01 6.4620e-02 2.7149e-01 3.2433e-01 6.9835e-01 8.2708e-01 9.4239e-01 7.7562e-01 ```

We can verify that the nodes are indeed the roots of the Nth orthogonal polynomial:

```norm( xi - roots(P(:,N+1)) ) ```
```ans = 1.4347e-15 ```

References:

[3] G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature rules, Math. Comp. 23 (1969), 221--230.