# Associated Legendre polynomials are not orthogonal

11 views (last 30 days)
Peter Harrington on 7 Sep 2014
Hi:
I am seeking an orthogonal set of polynomials, so I was excited to see the matlab had the legendre function to generate the polynomials. However, they look nothing like the polynomials plotted in wikipedia nor do they obey the rules of orthogonality that make these polynomials attractive.
Here is my code:
x = linspace(-1, 1, 1000);
y = legendre(5, x);
y*y'
One would expect the off-diagonal elements to be zero or close to zero but this does not occur. I am curious if there is an analytical fix that would make the polynomials orthogonal.
Also does anyone have some code to generate these polynomials correctly (I mean disassociated ;).
Pete

Roger Stafford on 7 Sep 2014
As cyclist has noted, the factor 1/(1-x^2) is essential in performing an orthogonality test. The Associated Legendre "polynomials" for differing m values are only orthogonal when each function is divided by sqrt(1-x^2).
Also note that when you make this change, your y*y' approximation to an integral would then yield NaNs at the two endpoints because of a zero-divided-by-zero occurrence. In the case of m1 = 0 and m2 = 1, you will get actual square root singularities in the integrand at the two ends even though the integral is still finite. To get accurate results for the orthogonality test I would advise you to use for-loops to compute with matlab's numerical quadrature function, 'integral', which can handle singularities, rather than your crude y*y' method. Either that or make an appropriate change of variable so that points are packed closely together at the two endpoints.

the cyclist on 7 Sep 2014
I hand-coded a few of the polynomials, for L = 3. I don't know what plots you are looking at on Wikipedia -- this page for Associated Legendre polynomials doesn't have any plot -- but these look good to me.
L = 3;
s = 10000;
x = linspace(-1, 1, s);
y = legendre(L, x);
P_3_0 = (5*x.^3 - 3*x);
P_3_1 = (3/2)*(1-5*x.^2).*sqrt(1-x.^2);
P_3_2 = 15*x.*(1-x.^2);
P_3_3 = -15*(1-x.^2).^(3/2);
figure
subplot(2,1,1), plot(x,y)
set(gca,'XLim',[-1 1])
set(gca,'YLim',[min(y(:)) max(y(:))])
subplot(2,1,2), plot(x,[P_3_0; P_3_1; P_3_2; P_3_3])
set(gca,'XLim',[-1 1])
set(gca,'YLim',[min(P_3_3) max(P_3_2)])
As to your comment about orthogonality, you seem to have coded the condition for fixed M, but the output from MATLAB you are using is fixed L, for different M. If you look at the Wikipedia page I referenced, in the orthogonality section, you'll see what I mean. To do the orthogonality check on
legendre(L,x)
as you are doing, looks like you need the
1./(1-x^2)
term.

Peter Harrington on 7 Sep 2014
Thanks for you reply Cyclist and Roger. I am a little bit confused with your replies. I am basing everything off this wikipedia page, http://en.wikipedia.org/wiki/Legendre_polynomials.
For the Legendre polynomials orthogonality requires the weighting function x = 1. I think that you are confusing the Legendre polynomials for the Chebyshev polynomial of the first kind that are orthogonal with the weight function 1/sqrt(1-x^2).
Although the inner product test of orthogonality may seem crude, it is practical because it is how the polynomials will be used. The ideal result should converge as the number of points increase and as one can see this case does not occur.
I also coded my own legendre dissociated function whose plots match the wiki article that I cited above. It still does not achieve the orthogonal results that I am seeking.
The function is below.
function y = LegPoly(x, df)
% y = LegPoly(x, df)
% Peter Harrington Peter.Harrington(at)Ohio.edu, 7-Sep-14
% start algorithm here
[m, n] = size(x); if m < n x = x'; m = n; end
y = zeros(m, df+1); y(:, 1) = ones(m, 1); y(:, 2) = x;
for i = 3:df+1 j = i-1; y(:, i) = ((2*j+1).*x.*y(:, j) - j*y(:, j-1))/i; end
Peter Harrington on 8 Sep 2014
Thanks again for your patience. I started with the associated polynomials, because that was the Matlab function, but then I switched to the normal legendre polynomials, which I coded up myself.
I think my problem was confusion between degree and order. I think that both your explanations have helped me resolve my confusion.