# Algorithm for orthogonalizing polynomials with specific inner product

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.

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..