JACOBI_ROOTS: Computes the n zeros of the Jacoby polynomial P_n^{(\alpha,\beta)}(x)
by Newton method and deflation process.
[x] = jacobi_roots(n,alpha,beta)
Input: n = polynomial degree
alpha,beta = parameters of Jacobi polynomial
Output: x = zeros of Jacobi polynomial (column array of size n)
flag = -1: The polynomial degree should be greater than 0
Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang,
"Spectral Methods. Fundamentals in Single Domains"
Springer Verlag, Berlin Heidelberg New York, 2006.0001 function [x,flag] = jacobi_roots(n,alpha,beta) 0002 % JACOBI_ROOTS: Computes the n zeros of the Jacoby polynomial P_n^{(\alpha,\beta)}(x) 0003 % by Newton method and deflation process. 0004 % 0005 % [x] = jacobi_roots(n,alpha,beta) 0006 % 0007 % Input: n = polynomial degree 0008 % alpha,beta = parameters of Jacobi polynomial 0009 % 0010 % Output: x = zeros of Jacobi polynomial (column array of size n) 0011 % flag = -1: The polynomial degree should be greater than 0 0012 % 0013 % Reference: CHQZ2 = C. Canuto, M.Y. Hussaini, A. Quarteroni, T.A. Zang, 0014 % "Spectral Methods. Fundamentals in Single Domains" 0015 % Springer Verlag, Berlin Heidelberg New York, 2006. 0016 0017 % Written by Paola Gervasio 0018 % $Date: 2007/04/01$ 0019 0020 0021 flag=0; 0022 if n < 1 0023 es='The polynomial degree should be greater than 0'; 0024 disp(es); flag = -1; x = []; 0025 return 0026 0027 else 0028 x=zeros(n,1); 0029 0030 x0=cos(pi/(2*n)); 0031 tol=1.e-14; kmax=15; 0032 for j=1:n 0033 diff=tol+1;kiter=0; 0034 while kiter <=kmax & diff>=tol 0035 [p,pd]=jacobi_eval(x0,n,alpha,beta); 0036 % deflation process q(x)=p(x)*(x-x_1)*... (x-x_{j-1}) 0037 % q(x)/q'(x)=p(x)/[p'(x)-p(x)*\sum_{i<j} 1/(x-x_i)] 0038 ss=sum(1./(x0-x(1:j-1))); 0039 x1=x0-p/(pd-ss*p); 0040 diff=abs(x1-x0); 0041 kiter=kiter+1; 0042 x0=x1; 0043 end 0044 x0=5.d-1*(x1+cos((2*(j+1)-1)*pi/(2*n))); 0045 x(j)=x1; 0046 end 0047 x=sort(x); 0048 end 0049 0050 0051 return