Plotting Legendre Polynomials - Getting different results for own method

Question

I'm trying to plot the legendre polynomials, defined as:

P0(x) = 1
P1(x) = x
Pn+1(x) = ((2n+1)/(n+1)) * x * Pn(x) - (n / (n+1)) * Pn-1(x)

I've done it the easy slow way, and I've done it the direct, a little more complicated way. Both result in a similar figure, however, not quite the same. The Amplitudes are different. Here is the code with the respect figure (Notice that I adjust the subscript n+1 of the definition to n):

xi = linspace(-1,1,500);
n = 10;

Method 1:

Pn = cell(n+1,1);
Pn{1} = @ (x) 1;
Pn{2} = @ (x) x;
for i=3:(n+1)
    Pn{i} = @ (x) ((2*(i-1)+1)/(i)).*x.*Pn{i-1}(x) - ((i-1)/i) .* Pn{i-2}(x);
end
plot(xi,Pn{1}(xi),'--r',xi,Pn{2}(xi),'--g',xi,Pn{3}(xi),'--b',...
    xi,Pn{4}(xi),'--m',xi,Pn{5}(xi),'--c',xi,Pn{6}(xi),'--k');
legend('P0','P1','P2','P3','P4','P5');

Figure 1:

Method 1

Method 2:

%Notice here that the results of Pj get stored into YI(j+1)
YI = zeros(length(xi),6);
YI(:,1) = ones(size(xi))';
YI(:,2) = xi';
for i=3:6;
    Pn1 = 1;
    Pn2 = xi;
    for j=2:(i-1) 
        Pn3 = ((2*(j-1)+1)/j) .* xi .* Pn2 - ((j-1) / j) .* Pn1;
        Pn1 = Pn2;
        Pn2 = Pn3;
    end
    YI(:,i) = Pn3';
end
figure('Name','direct method');
plot(xi,YI(:,1)','--r',xi, YI(:,2)', '--g', xi, YI(:,3)', '--b', ...
    xi,YI(:,4)','--m', xi,YI(:,5)', '--c', xi,YI(:,6)', '--k');

Figure 2:

Method 2

This is weird, to say the least. For method 1 I'm calculating all the polynomials up to P11, but only utilizing the first 6 to plot. Does anybody know what is going on?

And remember! DO NOT use i and j as variable names in Matlab. They are the complex unit and you may screw up the complex mathematics of Matlab. Use ii, jj or any other variable name you want, but no i or j. You make Matlab slower (it neds to sort out variable names more complexly) and destroy the math in any built in function using complex numbers — Ander Biguri, Jun 06 '14 at 09:37

Amro · Accepted Answer · 2014-06-06T14:32:56.687

2

Method#2 can be made much simpler:

X = linspace(-1,1,500); X = X(:);
N = 10;

Y = zeros(numel(X),N);
Y(:,1) = 1;
Y(:,2) = X;
for n=1:(N-1)
    Y(:,n+2) = ((2*n+1) .* X .* Y(:,n+1) - n .* Y(:,n)) / (n+1);
end

figure, plot(X, Y(:,1:6))
legend(num2str((1:6)'-1,'P_%d(x)'))
xlabel('x'), ylabel('P_n(x)'), title('Legendre Polynomials')

legendre_polynomials

This is equivalent to the plot shown on the Wikipedia page.

EDIT:

I had an off-by-one bug in the array indexing; MATLAB uses 1-based indexing, but the formula is defined in a 0-based manner. It is fixed now, sorry for the confusion ;)

formula

To confirm P(n=2,x=0) should be -1/2:

>> interp1(X, Y(:,3), 0)
ans =
   -0.5000

edited Jun 06 '14 at 14:32

answered Jun 06 '14 at 09:32

Amro

123,847
25
243
454

Thanks for the simplification, but that doesn't really answer my question. What is exactly wrong with my method? – spurra Jun 06 '14 at 12:21
I've went through the approach with a test value of 0, inputting it into P2. P2(t) is defined as 3/2 * t^2 - 1/2. My approach gives -0.5, which is correct, and yours does -2/3, which is obviously not according to above definition. However, all the legendre polynomials look like your curve. I am really confused now, I have no idea what is going on. – spurra Jun 06 '14 at 12:57
@BananaCode: the problem is with the second for-loop, it should be `for j=3:i` instead of `for j=2:(i-1)`. As you may have noticed now, your code is computing `P_n` by repeatedly starting from `P_0` going all the way up... Thanks to [dynamic programming](https://en.wikipedia.org/wiki/Dynamic_programming) method, you dont have to compute those values over and over again, just reuse the computation from the previous (two) iteration – Amro Jun 06 '14 at 12:57
But your values are wrong according to the definition above! Like I said, take P2 for instance. P2(0) should equal -0.5 but yours does -2/3. I have no doubt that you are correct, but this issue is confusing me. – spurra Jun 06 '14 at 13:06
Ok now I'm sure. You have your indexes wrong @Amro. My method 1 was wrong, and method 2 correct. But thanks for your simpler version! – spurra Jun 06 '14 at 13:35
@BananaCode: oh you are absolutely right! The code is now fixed. – Amro Jun 06 '14 at 14:34

Plotting Legendre Polynomials - Getting different results for own method

1 Answers1

EDIT: