restart;
read `GQ.map`;
# order of the problem
N := 11;
# floating-point precision
Digits := 50;
# interval [a, b] and weight function omega
# (will result in Chebyshev polynomials of the first kind)
a := -1;
b := 1;
omega := x -> (1-x^2)^(-1/2);
# inner product iprod
iprod := (f, g, x) -> int (f * g * omega (x), x=a..b):
# define the three-term recurrence relationship of the
# polynomials pi[k](x)
alphahat := array (1..2*N-1):
betahat := array (1..2*N-1):
gammahat := array (1..2*N-2):
for k from 1 to 2*N-1 do
  alphahat[k] := k;
  betahat[k] := 1;
od:
for k from 1 to 2*N-2 do
  gammahat[k] := 2;
od:
# constant polynomial pi[0]
pi0 := 1:
# compute the modified moments nu[k]
pi := array (0..2*N-1):
nu := array (0..2*N-1):
pi[0] := pi0:
nu[0] := int (omega (x) * pi[0], x=a..b):
pi[1] := expand ((x * pi[0] -
                  alphahat[1] * pi[0]) / betahat[1]):
nu[1] := int (omega (x) * pi[1], x=a..b):
for k from 2 to 2*N-1 do
  pi[k] := expand ((x * pi[k-1] -
                    gammahat[k-1] * pi[k-2] -
                    alphahat[k] * pi[k-1]) / betahat[k]):
  nu[k] := int (omega (x) * pi[k], x=a..b):
od:
# moment mu0 = iprod (1, 1, x)
mu0 := nu[0] / pi0:
# compute orthogonal polynomials
SackDonovan (pi0, alphahat, betahat, gammahat, nu, N,
             'alpha', 'beta');
# Lanczos (iprod, 'alpha', 'beta', N);
print (map (combine, alpha));
# latex (map (combine, alpha));
# map (x -> evalf (evalf (x, 20), 5), alpha);
print (map (combine, beta));
# latex (map (combine, beta));
# map (x -> evalf (evalf (x, 20), 5), beta);
# compute the Gauss quadrature rule
lprint (`Gauss quadrature rule`);
Gauss (mu0, alpha, beta, N, 'x', 'w');
TestQuadratureRule (iprod, x, w, 2*N-1);
# compute the Gauss-Radau quadrature rule, x1 = a
lprint (`Gauss-Radau quadrature rule, x1 = a`);
Radau (mu0, alpha, beta, N-1, a, 'x', 'w');
TestQuadratureRule (iprod, x, w, 2*N-2);
# compute the Gauss-Radau quadrature rule, x1 = b
lprint (`Gauss-Radau quadrature rule, x1 = b`);
Radau (mu0, alpha, beta, N-1, b, 'x', 'w');
TestQuadratureRule (iprod, x, w, 2*N-2);
# compute the Gauss-Lobatto quadrature rule, x1 = a, x2 = b
lprint (`Gauss-Lobatto quadrature rule, x1 = a, x2 = b`);
Lobatto (mu0, alpha, beta, N-1, a, b, 'x', 'w');
TestQuadratureRule (iprod, x, w, 2*N-1);