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);