Associated Legendre polynomials are not orthogonal

21 次查看（过去 30 天）
Peter Harrington
Peter Harrington 2014-9-7
评论： Peter Harrington 2014-9-8
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

请先登录，再回答此问题。

回答（3 个）

Roger Stafford
Roger Stafford 2014-9-7
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
the cyclist 2014-9-7
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
Peter Harrington 2014-9-7
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
  3 个评论
显示 1更早的评论
Roger Stafford
Roger Stafford 2014-9-7
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
Peter Harrington 2014-9-8
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Polynomials 的更多信息

标签

产品

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by