I've written some Matlab procedures that evaluate orthogonal polynomials, and as a sanity check I was trying to ensure that their dot product would be zero.
But, while I'm fairly sure there's not much that can go wrong, I'm finding myself with a slightly curious behaviour. My test is quite simple:
x = -1:.01:1;
for i0=0:9
v1 = poly(x, i0);
for i1=0:i0
v2 = poly(x,i1);
fprintf('%d, %d: %g\n', i0, i1, v1*v2');
end
end
(Note the dot product v1*v2' needs to be this way round because x is a horizontal vector.)
Now, to cut to the end of the story, I end up with values close to 0 (order of magnitude about 1e-15) for pairs of degrees that add up to an odd number (i.e., i0+i1=2k+1). When i0==i1 I expect the dot product not to be 0, but this also happens when i0+i1=2k, which I didn't expect.
To give you some more details, I initially did this test with Chebyshev polynomials of first kind. Now, they are orthogonal with respect to the weight
1 ./ sqrt(1-x.^2)
which goes to infinity when x goes to 1. So I thought that leaving this term out could be the cause of non-zero dot products.
But then, I did the same test with Legendre polynomials, and I get exactly the same result: when the sum of the degrees is even, the dot product is definitely far from 0 (order of magnitude 1e2).
One last detail, I used the trigonometric formula cos(n*acos(x)) to evaluate the Chebyshev polynomials, and I tried the recursive formula as well as one of the formulas involving the binomial coefficient to evaluate the Legendre polynomials.
Can anyone explain this odd (pun intended) behaviour?
You're being misled by symmetry. Both Chebyshev and Legendre polynomials are eigenfunctions of the parity operator, which means that they can all be classified as either odd or even functions. I guess the same goes for your custom orthogonal polynomials.
Due to this symmetry, if you multiply a polynomial P_n(x) by P_m(x), then the result will be an odd function if n+m is odd, and it will be even otherwise. You're computing sum_k P_n(x_k)*P_m(x_k) for a symmetric set of x_k values around the origin. This implies that for odd n+m you will always get zero. Try computing sum_k P_n(x_k)*Q_m(x_k) with P a Legendre, and Q a Chebyshev polynomial. My point is that for n+m=odd, the result doesn't tell you anything about orthogonality or the accuracy of your integration.
The problem is that probably you're not integrating accurately enough. These orthogonal polynomials defined on [-1,1] vary quite rapidly on their domain, especially close to the boundaries (x==+-1). Try increasing the points of your integration, using a non-equidistant mesh, or a proper integration using integral.
Final note: I'd advise you against calling your functions poly, since that's a MATLAB built-in. (And so is legendre.)