Algorithm for orthogonalizing polynomials with specific inner product

Question

I am attempting to generate a as big as possible collection of orthogonal polynomials $p_1, p_2, …, p_n$, $\left\langle p_i, q_i\right \rangle = \delta_{ij}$ where the inner product is with respect to a specific weight function.

It is well known that this can be done using the Gram-Schmidt process but my problem is that for large enough $n$ the resulting polynomials lose their orthogonality due to numerical instability. I have also implemented the Modified Gram-Schmidt algorithm but for around $n=30$ these also fail to be orthogonal.

My question is if there is a better method for generating these polynomials? I'm looking for a link or reference or description of the method.

Answer 1

I can't tell from what you said what you implemented, this is an example called [Chebfun]. 1. It'd be better if you posted your code. Note this is some code for the Legendgre polynomials with pseudo-code.

x = chebfun('x',[-1 1]);
w = exp(pi*x);
N = 5;
P = OrthPoly(w,N);

    function P = OrthPoly(w,N)
        if isnumeric(w), w = chebfun(w,[-1 1]); end
        d = w.ends;                     % the domain
        x = chebfun('x',d);             % linear chebfun
        P = chebfun(1./sqrt(sum(w)),d); % the constant (normalised)
        for k = 1:N;
            xk = x.*P(:,k);
            P(:,k+1) = xk;
            for j = 1:k       % Subtract out the components
                C = sum(w.*xk.*P(:,j));
                P(:,k+1) = P(:,k+1) - C*P(:,j);
            end
            P(:,k+1) = P(:,k+1)./sqrt(sum(w.*P(:,k+1).^2)); % normalise
        end
    end

It's actually the lanczos iteration..

Note this is from Trefethen..