Associated Legendre polynomials are not orthogonal

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

Thanks in advance,

Pete

Answer 1

[Roger Stafford]

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.

Answer 2

[the cyclist]

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.

Answer 3

[Peter Harrington]

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

[John D'Errico]

Um, I think the disconnect is in what you are looking for.

There is a difference between a Legendre polynomial and an associated Legendre polynomial. The first is related to the second in terms of a derivative formula.

https://en.wikipedia.org/wiki/Associated_Legendre_polynomials http://en.wikipedia.org/wiki/Legendre_polynomials

And while you ASK for an associated Legendre poly in the title of your post, you act as if you are looking for a simple Legendre polynomial. In fact, the link you post here is for exactly that.

So which is it? My guess is you really don't want the associated ones. Or maybe you do, and don't know it.

So for example, using my own sympoly toolbox to generate the basic Legendre polynomials, I get:

p1 = orthpoly(1,'leg')
p1 =
    x
p2 = orthpoly(2,'leg')
p2 =
    -0.5 + 1.5*x^2

defint(p1*p2,'x',[-1,1])
ans =
    0

defint(p1*p1,'x',[-1,1])
ans =
    0.66667

[Roger Stafford]

I agree with what John has said, but I will expand upon his remarks a little. Peter, you have stated, "For the Legendre polynomials orthogonality requires the weighting function x = 1." That is a true statement but only as applied to Associated Legendre polynomials of the same order, m, and different degrees, L1 ~= L2. With different orders m1 ~= m2 and the same degree, L, which is what you were testing in your original question, these polynomials need a weighting function of 1/sqrt(1-x^2) to achieve orthogonality. That is also a true statement.

That is where your tests with y*y' went wrong. You were comparing these polynomials, all with degree 5, and differing orders 0 through 5. That is what matlab's call

   y = legendre(5, x);

gives you. y will have six rows corresponding to orders 0 to 5, and each row an Associated Legendre polynomial of degree 5. The only one that is a Legendre polynomial of the kind referred to in

   http://en.wikipedia.org/wiki/Legendre_polynomials

is that in the first row with order zero.

To test for orthogonality with a weighting factor of 1, you need to call on 'legendre' a number of times with different degrees and compare rows having the same order - that is, equal row numbers.

   x = linspace(-1,1,1000);
   y5 = legendre(5,x); % Choose degrees 5 & 6
   y6 = legendre(6,x); 
   y6(7,:) = []; % Discard order 6 since y5 doesn't have that order
   sum(y5.*y6,2)

You should get a 6-element column vector of "zeros" to the accuracy of your approximation of the six definite integrals. At least they should be very small compared with, say, sqrt(sum(y5.*y5,2).*sum(y6.*y6,2)).

[Peter Harrington]

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.