Problem 10.6

Nasser Abbasi

 

OUTPUT

The program will return an array that contain all roots for a legendre polynomial or order n. So root m is the element m in the array. Comparing this to table 10.7, the results agree.

 

These are example runs for number of different order n.

 

» clear all

» r=nma_legendreRoot(1)

 

r =

 

     0

 

» r=nma_legendreRoot(2)

 

r =

 

   -0.5774    0.5774

 

» r=nma_legendreRoot(3)

 

r =

 

   -0.7746         0    0.7746

 

» r=nma_legendreRoot(4)

 

r =

 

   -0.8611   -0.3400    0.3400    0.8611

 

» r=nma_legendreRoot(5)

 

r =

 

   -0.9062   -0.5385         0    0.5385    0.9062

 

» r=nma_legendreRoot(6)

 

r =

 

   -0.9325   -0.6612   -0.2386    0.2386    0.6612    0.9325

 

» r=nma_legendreRoot(7)

 

r =

 

   -0.9491   -0.7415   -0.4058         0    0.4058    0.7415    0.9491

 

 

» r=nma_legendreRoot(8)

 

r =

 

   -0.9603   -0.7967   -0.5255   -0.1834    0.1834    0.5255    0.7967   0.9603

 

» r=nma_legendreRoot(9)

 

r =

 

   -0.9682   -0.8360   -0.6134   -0.3243         0    0.3243    0.6134    0.8360    0.9682

 

» r=nma_legendreRoot(10)

 

r =

 

 

   -0.9739   -0.8651   -0.6794   -0.4334   -0.1489    0.1489    0.4334   0.6794    0.8651    0.9739

 

» r=nma_legendreRoot(11)

 

r =

 

    -0.9782   -0.8871   -0.7302   -0.5191   -0.2695         0    0.2695   0.5191    0.7302    0.8871    0.9782

 

» r=nma_legendreRoot(12)

 

r =

 

   -0.9816   -0.9041   -0.7699   -0.5873   -0.3678   -0.1252    0.1252   0.3678    0.5873    0.7699    0.9041    0.9816

 

» r=nma_legendreRoot(13)

 

r =

 

 

 -0.9842   -0.9176   -0.8016   -0.6423   -0.4485   -0.2305         0     0.2305    0.4485    0.6423    0.8016    0.9176    0.9842

 

 

Source code

 

function roots_ = nma_legendreRoot(n)

%function nma_legendreRoot(n)

%

%Finds all roots for the legendre polynomial Pn(x).

%

%INPUT:

%  n: the order of the legendre polynomial

%

%OUTPUT

%  roots: a vector that contains all roots for this order

%         the length of the vector is the number of roots.

%

%by Nasser Abbbasi, HW 10.6

 

% Design:

%       for n=1, root 1 is zero.

%       roots are all real.

%       roots are interwindes. i.e.

%

%        first root for first order, r(1,1) =  0

%

%    so  first root for second order,  r(2,1) < r(1,1)

%        second root for second order, r(2,2) > r(1,1)

%

%    so  first root for third order,  r(3,1)   < r(2,1)

%        second root for third order, r(3,2)   > r(2,1) and less than r(2,2)

%        third root for third order,  r(3,3)   > r(2,2)

%

%    so  first root for 4th order,    r(4,1)   < r(3,1)

%        second root for 4th order,   r(4,2)   > r(3,1) and less than r(3,2)

%        third root for 4th order,    r(4,3)   > r(3,2) and less than r(3,3)

%        fourth root for 4th order,   4(4,4)   > r(3,3)

%

% i.e. we have a tree as shown below, where the number of nodes

% at each level is the order n of the polynomial, and the nodes

% are the roots for that polynomial order.

%

%

%                       x              n=1

%                   x       x          n=2

%               x       x      x       n=3

%           x       x       x    x     n=4

%

% So at each level, I make an array of guesses, which I'll use

% to find the roots at this level. At the next level, I'll use

% the roots found at the above level as the guess.

% This continues untill last level (leaves of the tree) is reached.

% 

% To find the roots from the guesses, Use Newton method.

% To evaluate Pn(x), use the recusrion formula 10.15 in the book,

% which is allready implemented in the function legndr.

%

% To evaluate P'n(x), use the formula given in the problem 10.6.

%


 

roots_=[];

 

if(n<=0)

   error('order of legendre polynomial must be greater than zero');

end

 

if(n==1)

   roots_(1)=0;

   return;

end

 

EPS=0.01;  % some small increment to add to x as we get close to the root.

 

r(1) = 0; 

maxNumberOfIter = 10;  % for netwon, stop after this many iters unless got a root.

 

for(i=2:n)

    % do the edge guesses first

    guess(1)   = r(1)-EPS;   

    guess(i)   = r(end)+EPS;

    %

    % now do the internal roots guesses, take half way between

    % above level 2 adjacent roots. see tree diagram above.

    %

    for(k=2:i-1)

        m       = abs( r(k) - r(k-1) )/2;

        guess(k)= m + r(k-1);

    end

    %

    % now that we have the guess for this level, find the

    % roots at this level. Then use these for the next level lower

    % roots guesses.

    % Use Newton method  Xn+1 = Xn - P/P’

    %

    for(N = 1:length(guess))   % find all roots from guesses

        x = guess(N);

        for(k = 1:maxNumberOfIter)

           pn   = legndr(i,x);  % supplied by prof.

           % note: pn(1) is p0, pn(2) is p1, etc...

           pd   = derivativeP(i,x,pn(i+1),pn(i));

           newX = x - (pn(i+1)/pd);             

           if(newX == x )

               x = newX;

               break;

           end

           x = newX;

        end

        thisLevelRoots(N) = x;

    end

    r = thisLevelRoots;  % for use by next iteration

end

 

%

% Ok, done, return the roots

%

roots_ = r;

return;

 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%

% Function that find P' for a given n,x

%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function pd = derivativeP(n,x,pn,pn_1)

 

   pd = (n*x*pn - n*pn_1)/(x^2-1);

 

return;