Lanczos := proc (iprod, alpha, beta, n) # # The procedure Lanczos computes the orthogonal polynomials # with respect to the inner product iprod. # # # Input: # iprod # n # Output: # alpha[k], k=1..n # beta[k], k=1..n-1 # # # The operator omega defines the inner product, e.g. # # b # / # | # iprod = (f,g,x) -> | f g dx # | # / # a # # # The procedure Lanczos computes the coefficients # alpha[k], beta[k] # of the three term recurrence relationship # x p[k-1] = beta[k-1] p[k-2] + alpha[k] p[k-1] + beta[k] p[k], # k=1..n # of the orthogonal polynomials p[k] with respect to the # weight function omega(x) on the interval [a,b], # where beta[0] = 0, p[-1] = 0. # local k, p, q, x; alpha := array (1..n); if n > 1 then beta := array (1..n-1) else beta := 'beta' fi; p[0] := collect (1 / sqrt (iprod (1, 1, x)), x, simplify); alpha[1] := normal (iprod (x*p[0], p[0], x)); q := collect ((x - alpha[1]) * p[0], x, simplify); for k from 2 to n do beta[k-1] := normal (sqrt (iprod (q, q, x))); p[k-1] := collect (q / beta[k-1], x, simplify); alpha[k] := normal (iprod (x*p[k-1], p[k-1], x)); q := collect ((x - alpha[k]) * p[k-1] - beta[k-1]*p[k-2], x, simplify); od; RETURN (NULL); end: