Radau := proc (mu0, alpha, beta, n, x1, x, w) # # The procedure Radau computes the Gauss-Radau quadrature rule # # b n + 1 # / ----- # | \ # | omega(x) f(x) dx = ) w[k] f(x[k]) # | / # / ----- # a k = 1 # # where the abscissa x1 on the boundary is 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 # beta[k], k=1..n # n # x1 # Output: # x[k], k=1..n+1 # w[k], k=1..n+1 # # # 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 denotes the prescribed abscissa (x1=a or x1=b). # # x[k] denotes the k. abscissa. # w[k] denotes the weight at the abscissa x[k]. # local alphan1tilde, D, e, i, Q, T, y; # set up the symmetric tridiagonal matrix T = Tn - x1 I T := array (1..n, 1..n, symmetric, [[0$n]$n]); for i from 1 to n do T[i,i] := alpha[i] - x1; od; for i from 1 to n-1 do T[i,i+1] := beta[i]; T[i+1,i] := beta[i]; od; # solve (Tn - x1 I) y = en for y[n] e := linalg[vector] (n, 0); e[n] := 1; y := linalg[linsolve] (T, e); # compute alphan+1tilde alphan1tilde := simplify (x1 + beta[n]^2 * y[n]); # set up the symmetric tridiagonal matrix T = Tn+1 T := array (1..n+1, 1..n+1, symmetric, [[0$n+1]$n+1]); for i from 1 to n do T[i,i] := alpha[i]; od; T[n+1,n+1] := alphan1tilde; for i from 1 to n do T[i,i+1] := beta[i]; T[i+1,i] := beta[i]; od; # compute the eigenvalue decomposition Tn+1 = QDQ' D := evalf (Eigenvals (T, Q)); # abscissas x := array (1..n+1, [seq (D[i], i=1..n+1)]); # weights w := array (1..n+1, [seq (evalf (mu0) * Q[1,i]^2, i=1..n+1)]); RETURN (NULL); end: