Orthogonal polynomials via the Lanczos process
Pedro Gonnet, 1st November 2011
(Chebfun example approx/OrthPolysLanczos.m)
format short e
The example approx/OrthPolys  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 .
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 . 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
 G. H. Golub and J. H. Welsch, Calculation of Gauss quadrature rules, Math. Comp. 23 (1969), 221--230.