Lobatto := proc (mu0, alpha, beta, n, x1, x2, x, w)
#
# The procedure Lobatto computes the Gauss-Lobatto quadrature rule
#
# b n + 2
# / -----
# | \
# | omega(x) f(x) dx = ) w[k] f(x[k])
# | /
# / -----
# a k = 1
#
# where the abscissas x1 and x2 on the boundary are prescribed.
# The weight function omega(x) and the interval [a,b] need not be
# specified directly. Rather the calculation is based on the orthogonal
# polynomials p[k](x) with respect to the weight function omega(x) on
# the interval [a,b]. These polynomials must be specified by their three
# term recurrence relationship.
#
#
# Input:
# mu0
# alpha[k], k=1..n+1
# beta[k], k=1..n
# n
# x1
# x2
# Output:
# x[k], k=1..n+2
# w[k], k=1..n+2
#
#
# mu0 denotes the moment
#
# b
# /
# |
# mu0 = | omega(x) dx
# |
# /
# a
#
# The coefficients
# alpha[k], beta[k]
# define the sequence of orthogonal polynomials p[k](x) by the three
# term recurrence relationship
# x p[k-1] = beta[k-1] p[k-2] + alpha[k] p[k-1] + beta[k] p[k],
# where beta[0] = 0, p[-1] = 0, p[0] = sqrt(1/mu0).
#
# x1 and x2 denote the prescribed abscissas (x1=a and x2=b).
#
# x[k] denotes the k. abscissa.
# w[k] denotes the weight at the abscissa x[k].
#
local alphan2tilde, betan1tilde, D, e, i, Q, T, y, z;
# set up the symmetric tridiagonal matrix T = Tn+1 - x1 I
T := array (1..n+1, 1..n+1, symmetric, [[0$n+1]$n+1]);
for i from 1 to n+1 do
T[i,i] := alpha[i] - x1;
od;
for i from 1 to n do
T[i,i+1] := beta[i];
T[i+1,i] := beta[i];
od;
# solve (Tn+1 - x1 I) y = en+1 for y[n+1]
e := linalg[vector] (n+1, 0); e[n+1] := 1;
y := linalg[linsolve] (T, e);
# set up the symmetric tridiagonal matrix T = Tn+1 - x2 I
for i from 1 to n+1 do
T[i,i] := alpha[i] - x2;
od;
# solve (Tn+1 - x2 I) z = en+1 for z[n+1]
z := linalg[linsolve] (T, e);
# compute alphan+2tilde and betan+1tilde
alphan2tilde := simplify ((x2 * y[n+1] - x1*z[n+1]) /
(y[n+1] - z[n+1]));
betan1tilde := simplify (sqrt ((x2-x1) /
(y[n+1] - z[n+1])));
# set up the symmetric tridiagonal matrix T = Tn+2
T := array (1..n+2, 1..n+2, symmetric, [[0$n+2]$n+2]);
for i from 1 to n+1 do
T[i,i] := alpha[i];
od;
T[n+2,n+2] := alphan2tilde;
for i from 1 to n do
T[i,i+1] := beta[i];
T[i+1,i] := beta[i];
od;
T[n+1,n+2] := betan1tilde;
T[n+2,n+1] := betan1tilde;
# compute the eigenvalue decomposition Tn+2 = QDQ'
D := evalf (Eigenvals (T, Q));
# abscissas
x := array (1..n+2, [seq (D[i], i=1..n+2)]);
# weights
w := array (1..n+2,
[seq (evalf (mu0) * Q[1,i]^2, i=1..n+2)]);
RETURN (NULL);
end: