High accuracy root finder of Legendre polynomials' derivatives?

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?

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
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
else
end if
do k = 1, n+1
end do
good = good*2*(m+beta+1)/(m+alpha+beta+2)
end do
end program Lobattotest

• In Mathematica, $$NSolve[JacobiP[38, 1, 1, x] == 0, x, WorkingPrecision -> 30]$$ – Moo Feb 6 '18 at 2:14
• @Moo Thanks. I also saw a website that gets the nodes and weights for Gauss-Legendre-Lobatto quadrature. – user5713492 Feb 6 '18 at 3:05

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$$

• Except you want the roots of the Jacobi polynomial $P_{38}^{(1,1)}(x)$, not the Legendre polynomial $P_{38}(x)$! – user5713492 Feb 5 '18 at 10:22
• Or in Matlab: syms x; vpa(solve(legendreP(38,x)),20). Or see jacobip if that's what is required. – horchler Feb 5 '18 at 16:18