High accuracy root finder of Legendre polynomials' derivatives?

Question

I need GLL (Gauss-Legendre-Lobatto) nodes for the Legendre-Galerkin-NI spectral method. It requires me to find the roots of the derivatives of Legendre polynomials.

My Matlab program calculates the coefficients of the derivatives of the Legendre polynomials just fine, but the issue is finding their roots. Matlab's built in roots function works fine up to 21st degree polynomials, but when N=22 and the coefficients reach 4*10⁷ (and rest are very badly scaled), it starts giving me imaginary roots.

Does anyone have ideas on how to find the roots for N>21, or an alternate way to find the GLL nodes?

Answer 1

Use a CAS that allows arbitrary precision, e.g. Maple or Mathematica.

In Maple, for example, here are the roots of the Legendre polynomial of order $38$, to $20$ digits accuracy:

> Digits:= 20:
  fsolve(orthopoly[P](38, x);

$$- 0.99804993053568761981,\; - 0.98973945426638557194,\;- 0.97484632859015350764,\;- 0.95346633093352959567,\;- 0.92574133204858439682,\;- 0.89185573900463221679,\;- 0.85203502193236218886,\;- 0.80654416760531681555,\;- 0.75568590375397068074,\;- 0.69979868037918435591,\;- 0.63925441582968170718,\;- 0.57445602104780708113,\;- 0.50583471792793110324,\;- 0.43384716943237648437,\;- 0.35897244047943501326,\;- 0.28170880979016526136,\;- 0.20257045389211670320,\;- 0.12208402533786741987,\;- 0.040785147904578239913,\; 0.040785147904578239913,\; 0.12208402533786741987,\; 0.20257045389211670320,\; 0.28170880979016526136,\; 0.35897244047943501326,\; 0.43384716943237648437,\; 0.50583471792793110324,\; 0.57445602104780708113,\; 0.63925441582968170718,\; 0.69979868037918435591,\; 0.75568590375397068074,\; 0.80654416760531681555,\; 0.85203502193236218886,\; 0.89185573900463221679,\; 0.92574133204858439682,\; 0.95346633093352959567,\; 0.97484632859015350764,\; 0.98973945426638557194,\; 0.99804993053568761981 $$

Answer 2

With Matlab you need some kind of toolbox to get jacobiP; not having this, we look at some stuff about Jacobi Polynomials and we find that $$P_0^{(\alpha,\beta)}(x)=1$$ $$P_1^{(\alpha,\beta)}(x)=\frac12(\alpha+\beta+2)x+\frac12(\alpha-\beta)$$ And then $$P_{n+1}^{(\alpha,\beta)}(x)=\frac{(2n+\alpha+\beta+1)\left[(2n+\alpha+\beta)(2n+\alpha+\beta+2)x+(\alpha+\beta)(\alpha-\beta)\right]P_n^{(\alpha,\beta)}(x)+2(n+\alpha)(n+\beta)(2n+\alpha+\beta+2)P_{n-1}^{(\alpha,\beta)}(x)}{2(n+1)(n+\alpha+\beta+1)(2n+\alpha+\beta)}$$ So we can calculate $P_n^{(\alpha,\beta)}(x)$ in a stable fashion via this recurrence relation. To apply Newton-Raphson to find the roots, we need the derivatives $$P_n^{(\alpha,\beta)\prime}(x)=\frac{n\left[\alpha-\beta-(2n+\alpha+\beta)x\right]P_n^{(\alpha,\beta)}(x)+2(n+\alpha)(n+\beta)P_{n-1}^{(\alpha,\beta)}(x)}{(2n+\alpha+\beta)\left(1-x^2\right)}$$ I couldn't seem to persuade the asymptotic formulas to get me good initial values for the roots, so I gave up and just went halfway between the roots of the next lower order Jacobi polynomial. My Matlab file:

% JacobiZ.m

alpha = 1;
beta = 1;
x = -2*Jacobi(1,alpha,beta,-1)/(Jacobi(1,alpha,beta,1)- ...
    Jacobi(1,alpha,beta,-1))-1;
N = 38;

for n = 2:N,
    x = ([-1 x]+[x 1])/2;
    for k = 1:length(x),
        err = 1.0e6;
        while abs(err) > 1.0e-10,
            [yn,yn_1] = Jacobi(n,alpha,beta,x(k));
            yp = (n*(alpha-beta-(2*n+alpha+beta)*x(k))*yn+ ...
                2*(n+alpha)*(n+beta)*yn_1)/((2*n+alpha+beta)*(1-x(k)^2));
            err = yn/yp;
            x(k) = x(k)-err;
        end
    end
end
for k=1:length(x)
    fprintf('%19.16f %.6e\n',x(k), Jacobi(N,alpha,beta,x(k)));
end

And the function to evaluate Jacobi polynomials:

% Jacobi.m

function [yn,yn_1] = Jacobi(n,alpha,beta,x);
yn = ((alpha+beta+2)*x+(alpha-beta))*0.5;
yn_1 = 1;
for k = 1:n-1,
    temp = ((2*k+alpha+beta+1)*((2*k+alpha+beta)*(2*k+alpha+beta+2)*x+ ...
        (alpha-beta)*(alpha+beta))*yn- ...
        2*(k+alpha)*(k+beta)*(2*k+alpha+beta+2)*yn_1)/ ...
        (2*(k+1)*(k+alpha+beta+1)*(2*k+alpha+beta));
    yn_1 = yn;
    yn = temp;
end

And the output:

-0.9952979292443489 1.763563e-14
-0.9842662807175033 2.947670e-14
-0.9670100764879885 1.284882e-14
-0.9436397649436017 -6.676346e-15
-0.9143033396902095 5.038752e-15
-0.8791863434793398 6.172471e-15
-0.8385108227781064 -2.141470e-15
-0.7925339526015519 -6.298440e-16
-0.7415464191473844 -1.259688e-16
-0.6858705850843138 -5.668596e-16
-0.6258584527552575 3.149220e-16
-0.5618894392947227 -5.038752e-16
-0.4943679781252536 4.408908e-16
-0.4237209621555510 -1.889532e-16
-0.3503950449141809 1.259688e-16
-0.2748538167143244 6.298440e-17
-0.1975748737189108 1.574610e-17
-0.1190467984449711 -3.936525e-17
-0.0397660708021819 -1.180957e-17
 0.0397660708021819 -1.180957e-17
 0.1190467984449711 -3.936525e-17
 0.1975748737189108 1.574610e-17
 0.2748538167143244 6.298440e-17
 0.3503950449141809 1.259688e-16
 0.4237209621555510 -1.889532e-16
 0.4943679781252536 4.408908e-16
 0.5618894392947227 -5.038752e-16
 0.6258584527552575 3.149220e-16
 0.6858705850843138 -5.668596e-16
 0.7415464191473844 -1.259688e-16
 0.7925339526015519 -6.298440e-16
 0.8385108227781064 -2.141470e-15
 0.8791863434793398 6.172471e-15
 0.9143033396902095 5.038752e-15
 0.9436397649436017 -6.676346e-15
 0.9670100764879885 1.284882e-14
 0.9842662807175033 2.947670e-14
 0.9952979292443489 1.763563e-14

I'm a little more comfortable in Fortran, so I wrote a program that gets the nodes and weights for Gauss-Jacobi-Lobatto quadrature. For Gauss-Legendre-Lobatto quadrature, $\alpha=\beta=0$ of course.

! jacobi.f90
module jacobimod
   use ISO_FORTRAN_ENV, only: wp => REAL128
   implicit none
   contains
      subroutine Jacobi(n,alpha,beta,x,p1,p0)
         integer n
         real(wp) alpha, beta
         real(wp) x
         real(wp) p1, p0
         real(wp) temp
         integer k
         if(n < 0) then
            p1 = 0
            p0 = 0
         else if(n == 0) then
            p1 = 1
            p0 = 0
         else
            p1 = ((alpha+beta+2)*x+(alpha-beta))/2
            p0 = 1
            do k = 1, n-1
               temp = ((2*k+alpha+beta+1)*((2*k+alpha+beta)*(2*k+alpha+beta+2)*x+ &
                  (alpha-beta)*(alpha+beta))*p1- &
                  2*(k+alpha)*(k+beta)*(2*k+alpha+beta+2)*p0)/ &
                  (2*(k+1)*(k+alpha+beta+1)*(2*k+alpha+beta))
               p0 = p1
               p1 = temp
            end do
         end if
      end subroutine Jacobi
      subroutine Lobatto(n, alpha, beta, x, lambda)
         integer n
         real(wp) alpha, beta
         real(wp) x(0:n+1)
         real(wp) lambda(0:n+1)
         integer m, k
         real(wp) p0, p1
         real(wp) yp, err
         real(wp) b
         b = 1
         lambda(0) = 2**(alpha+beta+1)*gamma(alpha+1)* &
            gamma(beta+1)/gamma(alpha+beta+3)
         do m = 1, n
            b=b*(m+alpha+1)/(m+beta+1)
            lambda(0) = lambda(0)*m/(m+alpha+beta+2)
         end do
         x(0) = -1
         lambda(n+1) = lambda(0)*(beta+1)/b
         lambda(0) = lambda(0)*(alpha+1)*b
         call Jacobi(1,alpha+1,beta+1,real(-1,wp),x(1),p0)
         call Jacobi(1,alpha+1,beta+1,real(1,wp),p1,p0)
         x(1) = -2*x(1)/(p1-x(1))-1
         x(2:n+1) = 1
         b = 2**(alpha+beta+2)*(2*n+alpha+beta+2)*gamma(alpha+1)* &
            gamma(beta+1)/gamma(alpha+beta+3)
         do m = 1, n
            b = b*(m+alpha)*(m+beta)/(m*(m+alpha+beta+2))
         end do
         do m = 2, n
            x(1:m) = (x(0:m-1)+x(1:m))/2
            do k = 1, m
               err = 1
               do while(abs(err) > sqrt(epsilon(err))/100)
                  call Jacobi(m,alpha+1,beta+1,x(k),p1,p0)
                  yp = (m*(alpha-beta-(2*m*alpha+beta+2)*x(k))*p1+ &
                     2*(m+alpha+1)*(m+beta+1)*p0)/((2*m+alpha+beta+2)*(1-x(k)**2))
                  err = p1/yp
                  x(k) = x(k)-err
               end do
               call Jacobi(m,alpha+1,beta+1,x(k),p1,p0)
               yp = (m*(alpha-beta-(2*m*alpha+beta+2)*x(k))*p1+ &
                  2*(m+alpha+1)*(m+beta+1)*p0)/((2*m+alpha+beta+2)*(1-x(k)**2))
               lambda(k) = b/((1-x(k)**2)*yp*p0)
            end do
         end do
      end subroutine Lobatto
end module jacobimod

program Lobattotest
   use jacobimod
   implicit none
   integer n
   real(wp), allocatable :: x(:), lambda(:)
   real(wp) alpha, beta
   integer m, k
   real(wp) bad, good
   write(*,'(a)',advance='no') 'Please enter the value of alpha:> '
   read(*,*) alpha
   write(*,'(a)',advance='no') 'Please enter the value of beta:> '
   read(*,*) beta
   write(*,'(a)',advance='no') 'Please enter the value of n:> '
   read(*,*) n
   allocate(x(0:n+1),lambda(0:n+1))
   call Lobatto(n, alpha, beta, x, lambda)
   do m = 0, n+1
      write(*,*) x(m), lambda(m)
   end do
   good = 2**(alpha+beta+1)*gamma(beta+1)*gamma(alpha+1)/gamma(alpha+beta+2)
   do m = 0, 2*n+2
      if(m == 0) then
         bad = lambda(0)
      else
         bad = 0
      end if
      do k = 1, n+1
         bad = bad+lambda(k)*(1+x(k))**m
      end do
      write(*,*) m,good-bad
      good = good*2*(m+beta+1)/(m+alpha+beta+2)
   end do
end program Lobattotest