C ALGORITHM 840, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 31, NO. 1, March, 2005, P. 149--165. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Matlab/ # Matlab/Dp/ # Matlab/Dp/Drivers/ # Matlab/Dp/Drivers/Res # Matlab/Dp/Drivers/TP_Lobptsts_MAIN.m # Matlab/Dp/Drivers/TP_LobptswtsPERT.m # Matlab/Dp/Src/ # Matlab/Dp/Src/PN_Legendre_vectN.m # Matlab/Dp/Src/PNx_Legendre_vectN.m # Matlab/Dp/Src/PNxx_Legendre_vectN.m # Matlab/Dp/Src/TP_LobattoQuad.m # Matlab/Dp/Src/TP_Lobattoresid.m # Matlab/Dp/Src/TP_NewtonJAC.m # Matlab/Dp/Src/TP_eigprolODE.m # This archive created: Thu May 12 16:45:55 2005 export PATH; PATH=/bin:$PATH if test ! -d 'Matlab' then mkdir 'Matlab' fi cd 'Matlab' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test -f 'Res' then echo shar: will not over-write existing file "'Res'" else cat << "SHAR_EOF" > 'Res' Reference solution: N=6, c=5 xprol = -1.0000 -0.7456 -0.2705 0.2705 0.7456 1.0000 wprol = 0.0756 0.3939 0.5305 0.5305 0.3939 0.0756 Deriv = -4.8506 6.2863 -1.8652 1.1596 -1.1654 0.7225 -2.5163 1.1613 1.7431 -0.8619 0.7941 -0.4747 1.1400 -2.6451 0.3814 1.8843 -1.3105 0.7150 -0.7150 1.3105 -1.8843 -0.3814 2.6451 -1.1400 0.4747 -0.7941 0.8619 -1.7431 -1.1613 2.5163 -0.7225 1.1654 -1.1596 1.8652 -6.2863 4.8506 Deriv2d = 2.0605 -12.9281 13.8273 -10.1868 10.9431 -6.9792 13.5207 -22.5485 11.0554 -3.3876 2.5002 -1.3490 -1.9454 8.7203 -13.0121 8.1305 -3.2059 1.3910 1.3910 -3.2059 8.1305 -13.0121 8.7203 -1.9454 -1.3490 2.5002 -3.3876 11.0554 -22.5485 13.5207 -6.9792 10.9431 -10.1868 13.8273 -12.9281 2.0605 Eigenvalue check: c=5, first six eigenmodes [eigenvalues ] chi_array = 4.1951 12.9117 20.1769 26.5874 33.8971 43.3590 Eigenfunctions B_array = -0.8688 0.0000 -0.4790 -0.0000 0.1256 0.0000 0.0000 0.8425 0.0000 -0.5321 0.0000 0.0842 0.4823 -0.0000 -0.7611 -0.0000 0.4307 0.0000 0.0000 -0.5261 -0.0000 -0.7792 -0.0000 0.3388 -0.1114 -0.0000 0.4299 -0.0000 0.8514 -0.0000 0.0000 0.1149 -0.0000 0.3275 -0.0000 0.9074 0.0141 0.0000 -0.0802 0.0000 -0.2698 -0.0000 -0.0000 -0.0131 0.0000 -0.0497 0.0000 -0.2326 -0.0011 -0.0000 0.0079 -0.0000 0.0337 -0.0000 0.0000 0.0009 -0.0000 0.0041 -0.0000 0.0248 0.0001 0.0000 -0.0005 0.0000 -0.0024 0.0000 -0.0000 -0.0000 0.0000 -0.0002 0.0000 -0.0015 -0.0000 0.0000 0.0000 0.0000 0.0001 -0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0001 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 -0.0000 0.0000 -0.0000 0.0000 -0.0000 0 0 0 0 0 0 0 0 0 0 0 0 Reference solution: N=13, c=0.5 c*(13)=9.8175 xprol = Columns 1 through 6 -1.0000e+00 -9.4951e-01 -8.3613e-01 -6.7122e-01 -4.6808e-01 -2.4011e-01 Columns 7 through 12 0 2.4011e-01 4.6808e-01 6.7122e-01 8.3613e-01 9.4951e-01 Column 13 1.0000e+00 wprol = Columns 1 through 6 1.3955e-02 8.3590e-02 1.4127e-01 1.8630e-01 2.1774e-01 2.3609e-01 Columns 7 through 12 2.4211e-01 2.3609e-01 2.1774e-01 1.8630e-01 1.4127e-01 8.3590e-02 Column 13 1.3955e-02 ************************************** Reference solution: N=20, c= 0.99 c*(20) xprol = Columns 1 through 6 -1.0000e+00 -9.7286e-01 -9.1293e-01 -8.2742e-01 -7.2322e-01 -6.0578e-01 Columns 7 through 12 -4.7913e-01 -3.4621e-01 -2.0923e-01 -6.9989e-02 6.9989e-02 2.0923e-01 Columns 13 through 18 3.4621e-01 4.7913e-01 6.0578e-01 7.2322e-01 8.2742e-01 9.1293e-01 Columns 19 through 20 9.7286e-01 1.0000e+00 wprol = Columns 1 through 6 7.5508e-03 4.4672e-02 7.3949e-02 9.5900e-02 1.1160e-01 1.2261e-01 Columns 7 through 12 1.3020e-01 1.3527e-01 1.3838e-01 1.3986e-01 1.3986e-01 1.3838e-01 Columns 13 through 18 1.3527e-01 1.3020e-01 1.2261e-01 1.1160e-01 9.5900e-02 7.3949e-02 Columns 19 through 20 4.4672e-02 7.5508e-03 SHAR_EOF fi # end of overwriting check if test -f 'TP_Lobptsts_MAIN.m' then echo shar: will not over-write existing file "'TP_Lobptsts_MAIN.m'" else cat << "SHAR_EOF" > 'TP_Lobptsts_MAIN.m' % TP_Lobptsts_MAIN.m % input: N = number of grid points % c = bandwidth parameter c % output: xprol = vector of length(1,N) containing the prolate-Lobatto grid points % wprol = vector of length(1,N) containing the prolate-Lobatto quadrature weights % Deriv(i,j) = d C_{j}/dx(x_i) % Deriv2(i,j) = d**2 C_{j}/dx**2 (x_i) clf,clear N=input('The number of interpolation points N=') cstar= (pi/2) * (( N-1) + 0.5) c=input('bandwidth parameter c=') tic [xprol,wprol,Deriv,Deriv2d]=TP_LobptswtsPERT(N,c); CPUtime = toc SHAR_EOF fi # end of overwriting check if test -f 'TP_LobptswtsPERT.m' then echo shar: will not over-write existing file "'TP_LobptswtsPERT.m'" else cat << "SHAR_EOF" > 'TP_LobptswtsPERT.m' function [xprol,wprol,Deriv,Deriv2d]=TP_LobptswtsPERT(N,c) % is called by Main program; % input: N = number of grid points % c = bandwidth parameter c % output: xprol = vector of length(1,N) containing the prolate-Lobatto grid points % wprol = vector of length(1,N) containing the prolate-Lobatto quadrature weights % Deriv(i,j) = d C_{j}/dx(x_i) % Deriv2(i,j) = d**2 C_{j}/dx**2 (x_i) cstar= (pi/2) * (( N-1) + 0.5) if c > cstar disp('WARNING! c > cstar! Function returns NULL result') else % Approximate Legendre-Lobatto points (i. e., prolate points and weights for c=0) ta=pi*( 0:(1/(N-1)):1); xguess = - cos(ta); wguess = ( pi/N) * (1- xguess.^2).^(1/2); wguess(1)= 2/( (N-1)*N); wguess(N)=2/( (N-1)*N); if N <= 12 % direct call % [xprol,wprol,Deriv,Deriv2d]= TP_Lobptwtder2dGUESS(N,c,xguess,wguess); [xprol,wprol]=TP_LobattoQuad(N,c,xguess,wguess); else [xL,wL]=TP_LobattoQuad(N,0,xguess,wguess); [xprol1,wprol1]=TP_LobattoQuad(N,cstar/10,xguess,wguess); xguess= xL + 100*(xprol1 - xL) * (c/cstar)^2; wguess=wL + 100*(wprol1 - wL) * (c/cstar)^2; [xprol,wprol]=TP_LobattoQuad(N,c,xguess,wguess); end % if end % if on whether c is too large or not % ******************************************************************* % Next compute derivative matrix. [chi_array,B_array]=TP_eigprolODE(c,N-1); % B_array contains the NORMALIZED Legendre coefficients of the % NORMALIZED eigenfunctions. The j-th column is the % eigenfunction corresponding the chi_array(j), i. e., % psi_{j-1}(x; c). [psicoll,dummy]=size(B_array) ; PN2d=PN_Legendre_vectN(xprol,psicoll - 1) ; % evaluate normalized Legendre % at points xproll; PN2d(ii,j)=P_{j-1}(xprol(ii)). PSI_Gram= PN2d * B_array; % PSI_Gram(i, j) = psi( x_i, j ); matrix Psi in notes. PNx2d=PNx_Legendre_vectN(xprol,psicoll-1); PNxx2d=PNxx_Legendre_vectN(xprol,psicoll-1); PSI_x= PNx2d * B_array; % matrix Psi^{(x)} in notes % PSI_x(i,j) = d psi/dx(x_i, j) PSI_xx= PNxx2d * B_array; % matrix Psi^{(xx)} in notes % PSI_xx(i,j) = d**2 psi/dx**2(x_i, j) C_matrix=inv(PSI_Gram); % Matrix C in the noted. Deriv = PSI_x * C_matrix; % First derivatives of prolate CARDINAL basis. Deriv2d = PSI_xx * C_matrix; % Second derivatives of prolate CARDINAL basis. SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test -f 'PN_Legendre_vectN.m' then echo shar: will not over-write existing file "'PN_Legendre_vectN.m'" else cat << "SHAR_EOF" > 'PN_Legendre_vectN.m' function PNa=PN_Legendre_vectN(x,n) %This evaluates the NORMALIZED, ORTHOGONORMAL Legendre polynomials P_{n)(x) % for all degrees up to and including n. % x may be either a scalar or a vector % Input: x= scalar or vector of grid points where Legendre polynomials % are to be evaluated. % n=degree of highest Legendre polynomial needed. % % Output: PNa is a size(x) times (n+1) array % Example: let x = [0 0.3 0.9], n=3. Then the output will be the 3 x 4 array % PNa= | P_0(x(1)=0) P_{1}(x(1)) P_{2}(x(1)) P_{3}(x(1)) | % | P_0(x(2)=0.3) P_{1}(x(2)) P_{2}(x(2)) P_{3}(x(2)) | % | P_0(x(3)=0.9) P_{1}(x(3)) P_{2}(x(3)) P_{3}(x(3)) | PNa=zeros(length(x),n+1); PNa(:,1)=1; if n > 0 PNa(:,2)=x'; end % if if n > 1 for j=1:(n-1) PNa(:,j+2) = (1/(j+1))*( (2*j+1)*x' .*PNa(:,j+1) - j*PNa(:,j)); end % j end % if for j=1:(n+1) PNa(:,j)=PNa(:,j) * sqrt( j - 1/2); end % j SHAR_EOF fi # end of overwriting check if test -f 'PNx_Legendre_vectN.m' then echo shar: will not over-write existing file "'PNx_Legendre_vectN.m'" else cat << "SHAR_EOF" > 'PNx_Legendre_vectN.m' function PNxa=PNx_Legendre_vectN(x,n) %This evaluates the FIRST DERIVATIVE of the % NORMALIZED, ORTHOGONORMAL Legendre polynomials P_{n)(x) % for all degrees up to and including N. % x may be either a scalar or a vector % Input: x= scalar or vector of grid points where Legendre polynomials % are to be evaluated. % n=degree of highest Legendre polynomial needed. % % Output: PNxa is a size(x) times (n+1) array % Example: let x = [0 0.3 0.9], n=3. Then the output will be the 3 x 4 array % PNxa= | dP_0/dx(x(1)=0) dP_{1}/dx(x(1)) dP_{2}/dx(x(1)) dP_{3}/dx(x(1)) | % | dP_0/dx(x(2)=0.3) dP_{1}/dx(x(2)) dP_{2}/dx(x(2)) dP_{3}/dx(x(2)) | % | dP_0/dx(x(3)=0.9) dP_{1}/dx(x(3)) dP_{2}/dx(x(3)) dP_{3}/dx(x(3)) | % Checked by LegendreFirstDeriv.maple; % Recurrence is that for Gegenbauer polynomials of order 3/2 PNxa=zeros(length(x),n); PNxa(:,1)=0; if n > 0 PNxa(:,2)=1 ; end % if if n > 1 PNxa(:,3)=3*x'; end % if if n > 2 for j=1:(n-2) PNxa(:,j+3) = (1/(j+1)) * ( 2*(j+3/2)*x' .*PNxa(:,j+2) - (j+2)*PNxa(:,j+1) ) ; end % j end % if for j=1:(n+1) PNxa(:,j)=PNxa(:,j) * sqrt( j-1/2); end % j SHAR_EOF fi # end of overwriting check if test -f 'PNxx_Legendre_vectN.m' then echo shar: will not over-write existing file "'PNxx_Legendre_vectN.m'" else cat << "SHAR_EOF" > 'PNxx_Legendre_vectN.m' function PNxxa=PNxx_Legendre_vectN(x,n) %This evaluates the SECOND DERIVATIVE of the % NORMALIZED, ORTHOGONORMAL Legendre polynomials P_{n)(x) % for all degrees up to and including N. % x may be either a scalar or a vector % Input: x= scalar or vector of grid points where Legendre polynomials % are to be evaluated. % n=degree of highest Legendre polynomial needed. % % Output: PNxxa is a size(x) times (n+1) array % Example: let x = [0.1 0.3 ], n=3. Then the output will be the 2 x 4 array % PNxxa= | P_{0,xx}(x(1)=0.1) P_{1,xx}(x(1)) P_{2,xx}(x(1)) P_{3,xx}(x(1)) | % | P_{0,xx}(x(2)=0.3) P_{1,xx}(x(2)) P_{2,xx}(x(2)) P_{3,xx}(x(2)) | % P"_{n+3} = (1/(n+1)) * { 2 (n+5/2)*x*P"_{n+2} - (n+4) P"_{n+1} } % This is the recurrence for Gegenbauer of degree (n+1) and order (5/2) % on pg. 503 of my book. % Checked by LegendreSecondDeriv.maple; PNxxa=zeros(length(x),n); PNxxa(:,1)=0; if n > 0, PNxxa(:,2)=0; end % if if n > 1, PNxxa(:,3)=3; end % if if n > 2, PNxxa(:,4)=15*x' ; end % if if n > 3 for j=1:(n-3) PNxxa(:,j+4) = (1/(j+1)) * ( 2*(j+5/2)*x' .*PNxxa(:,j+3) - (j+4)*PNxxa(:,j+2) ) ; end % j end % if for j=1:(n+1), PNxxa(:,j)=PNxxa(:,j) * sqrt( j-1/2); end % j SHAR_EOF fi # end of overwriting check if test -f 'TP_LobattoQuad.m' then echo shar: will not over-write existing file "'TP_LobattoQuad.m'" else cat << "SHAR_EOF" > 'TP_LobattoQuad.m' function [xprol,wprol]=TP_LobattoQuad(N,c,xguess,wguess); global exact_integral psicoeffsymm; % TP_LobattoQuad computes the prolate-Lobatto-Gaussian quadrature points and weights % input: N = number of collocation points % c = bandwidth parameter % xguess --- vector of length(1,N) containing first guess for the Lobatto-grid points % wguess -- vector of length(1,N) containing first guess for Lobatto-prolate quadrature weights. % % output: xprol is vector of size(1,N) with the final prolate-Lobatto points % wprol is vector of size(1,N), the corresponding quadrature points M = (N/2) - 1 - (1/2)*rem(N,2); %if mod(N,2)==0 % npts is even %M= N/2 - 1 %else % npts is ODD %M= (N-3)/2; %end % if % M is the number of non-negative quadrature abscissas % xandw0 = vector of length(N-1) containing a first guess for the unknowns in % the Newton iteration. The first M entries are the grid points on the % interval 0 < x_k < 1. The remaining entries are the quadrature weights % associated with the nonnegative collocation points, i. e., weights for % x=1 and, if N is odd, x=0 are included if mod(N,2)==0 % N is even xnonneg=xguess(M+2:2*M+1); wnonneg=wguess(M+2:N); else xnonneg=xguess(M+3:N-1); wnonneg=wguess(M+2:N); end % if xandw0 = [xnonneg wnonneg]; % step one: compute the exact integral of each prolate function % of degree 0, 1, ... (2*npts-1) % Because of symmetry, we only use the integrals of EVEN prolate functions. nmax=2*N; % degree of higest prolate function that will be integrated exactly. [chi_array,B_array]=TP_eigprolODE(c,nmax); [psirow,psicol]=size(B_array) ; psicoeffsymm = B_array(1:2:psirow,1:2:(2*(N-1)-1) ) ; % Because the integral of overline(P_{j}) on [-1,1] is the inner product of % P_{j} with itself and the Legendre are orthogonal, the integral % is just the constant in the prolate normalized Legendre series exact_integral = sqrt(2) * psicoeffsymm(1,1:N-1); % apply a Newtoniteration until the quadrature is exact NewtMax=20; nalpha=4; % nalpha=0; JacMod=1; % [xandw,residratio]=TP_Newton(N-1,nalpha,NewtMax,JacMod,xandw0,lambda,'TP_Lobattoresid'); xandw=TP_NewtonJac(xandw0,c,'TP_Lobattoresid'); % process roots: xandw is a vector of dimension (N-1) contains only nonnegative x % and quadrature weights if mod(N,2)==0 % N is even xprol=[-1 -xandw(M:-1:1) xandw(1:M) 1] ; wprol=[xandw(N-1:-1:M+1) xandw(M+1:N-1)] ; else % N is ODD xprol=[-1 -xandw(M:-1:1) 0 xandw(1:M) 1]; wprol=[xandw(2*M+2:-1:M+2) xandw(M+1:2*M+2)]; end % if % xprol, wprol are each N-dimensional vectors containing % all grid points and weights, respectively. SHAR_EOF fi # end of overwriting check if test -f 'TP_Lobattoresid.m' then echo shar: will not over-write existing file "'TP_Lobattoresid.m'" else cat << "SHAR_EOF" > 'TP_Lobattoresid.m' function resid=TP_Lobattoresid(xprolwprol, c) global exact_integral psicoeffsymm ; global xproltemp; N=length(xprolwprol) + 1; resid = exact_integral; if mod(N,2)==1 M=length(xprolwprol)/2 - 1; % One grid point is always % at x=1 for Lobatto grid % M is the number of non-negative, non-unit grid points xprol=xprolwprol(1:M); wprol=xprolwprol(M+1:2*M+2); else M = (N-2)/2; % N is even xprol=xprolwprol(1:M); wprol=xprolwprol(M+1:2*M+1); end % if % First step: compute grid point values of the Legendre polynomials. [NEIG,psicoeffcoll]= size(psicoeffsymm); if mod(N,2)==1 xprolfull = [0 xprol 1]; % else % N even % disp('N is even !!!!!!!!!!!!!!!') xprolfull=[xprol 1] ; % xprolfull includes the fixed grid pt., x=1 end % if PN2dquad=PN_Legendre_vectN(xprolfull,2*NEIG); % evaluate normalized Legendre % at points xprolfull; P+N2d(ii,j)=P_{j-1}(xprolfull(ii)). PN2dquadsymm=PN2dquad(:,1:2:2*NEIG); % grid point values of symetric functions only % compute grid point values of the prolate functions at the quadrature points psi_j = PN2dquadsymm * psicoeffsymm; %psi_j = psi( x_i, j ) if mod(N,2)==0 approx_integral=2* psi_j' * wprol' ; else % if N is odd, then xprol(1)=0, and we must not double-count it. %approx_integral=2* (psi_j(:,2:N) )' * (wprol (2:N))' + psi_j(:,1)*wprol(1); approx_integral=zeros(N-1,1); for j=1:N-1 approx_integral(j)=psi_j(1,j)*wprol(1); % disp([' j=',int2str(j),' ',num2str(approx_integral(j)),' psi(j,1)=',num2str(psi_j(1,j)), 'wprol(1)=',num2str(wprol(1)) ]) for k=2:length(wprol) approx_integral(j)= approx_integral(j) + 2*psi_j(k,j) * wprol(k); % disp([num2str(k),' ',num2str(approx_integral(j)),num2str(2*psi_j(k,j)*wprol(k)),' psi(j,k)=',num2str(psi_j(k,j)), 'wprol(k)=',num2str(wprol(k)) ]) end % k end % j end % if resid= exact_integral' - approx_integral; SHAR_EOF fi # end of overwriting check if test -f 'TP_NewtonJAC.m' then echo shar: will not over-write existing file "'TP_NewtonJAC.m'" else cat << "SHAR_EOF" > 'TP_NewtonJAC.m' function xandw=TP_NewtonJac(xandw0,c,Res); % input: xandw=vector of (N-1) unknowns, first guess for start % of march. This vector stores both grid points and quadrature weights % Res = string variable with name of residual function, not % including the termination .m. % (Res must return a vector of % length N-1, but as a function of the arguments (xandw,lambda) % where xandw is a ROW vector and c is a scalar % rhist(1:NewtMax) stores the norms of the residuals at each iteration. % The numbers in rhist must decrease rapidly if the iteration is converging. % rhist is graphed at each Newton iteration so that one can follow the % convergence in ``real time''. % % Output: xandw stores the (N-1) unknowns. By symmetry, and using % known values (two of the grid points are always 1 and -1, and x=0 % is also always a grid point if N is odd), one can deduce the total % of N grid points and N quadrature weights from xandw. This % expansion to N grid points and N weights is done in the % calling program global exact_integral psicoeffsymm; Nmin1=length(xandw0) N=Nmin1 + 1; % N is the number of grid pts; Nmin1 is the number of unknowns if mod(N,2)==0 % npts is even M= N/2 - 1 else % npts is ODD M= (N-3)/2; end % if % Parameters controlling the Newton iteration % nalpha= number of stages in powers-of-two minimum % residual % NewtMax is the maximum number of Newton iterations to try. % JacMod: Jacobian matrix is recomputed every JacMod steps % JacMod=1 for classical Newton method % JacMod > NewtMax computes the Jacobian just once NewtMax=20; % maximum number of Newton iterations nalpha=4; % smallest value of relaxation parameter gamma tested is % 2**(nalpha) JacMod=1; % Jacmod==1 means that the Jacobian is updated at each % iteration Jacobian=zeros(Nmin1,Nmin1); xandw=xandw0; rhist=0; % Compute Jacobian matrix; deltanorm=1.E50; haha=1.E-5*max(abs(xandw)); if haha==0, haha=1.E-5, end % if anorm=1.E50; deltanorm=1.E50; for inewt=1:NewtMax % $$$$$$$$$$$$$$$ BEGIN NEWTONS if deltanorm > 1.E-10*anorm residual=feval(Res,xandw,c); % if rem(inewt-1,JacMod)==0 % Jacobiancheck .... compute by finite differences in unknowns % for jcol=1:Nmin1 % atemp=xandw; % atemp=n-dim column vector % atemp(jcol)=xandw(jcol)+haha; % Jaccol=(feval(Res,atemp,c)- residual) / haha; % Jacobiancheck(1:Nmin1,jcol)=Jaccol; % end if mod(N,2)==1 xunk=xandw(1:M); wunk=xandw(M+1:2*M+2); else M = (N-2)/2; % N is even xunk=xandw(1:M); wunk=xandw(M+1:2*M+1); end % if % M is the number of interior grid points on one-half % of the interval, always excluding x=1 and x=0. % First step: compute grid point values of the Legendre polynomials. [NEIG,psicoeffcoll]= size(psicoeffsymm); if mod(N,2)==1 xunkfull = [0 xunk 1]; % else % N even % disp('N is even !!!!!!!!!!!!!!!') xunkfull=[xunk 1] ; % xunkfull includes the fixed grid pt., x=1 end % if PN2dquad=PN_Legendre_vectN(xunkfull,2*NEIG); % evaluate normalized Legendre % at points xunkfull; P+N2d(ii,j)=P_{j-1}(xunkfull(ii)). PN2dquadsymm=PN2dquad(:,1:2:2*NEIG); % grid point values of symetric functions only % compute grid point values of the prolate functions at the quadrature points psi_j = PN2dquadsymm * psicoeffsymm; %psi_j = psi( x_i, j ) PNx2dquad=PNx_Legendre_vectN(xunkfull,2*NEIG); % evaluate normalized Legendre % at points xunkfull; PNx2d(ii,j)=P_{j-1}(xunkfull(ii)). PNx2dquadsymm=PNx2dquad(:,1:2:2*NEIG); % grid point values of symetric functions only % compute grid point values of the prolate functions at the quadrature points psix_j = PNx2dquadsymm * psicoeffsymm; %psix_j = d psi/dx( x_i, j ) iparity=rem(N,2); for jcol=1:M Jacobian(:,jcol) = - wunk(jcol+iparity) * psix_j(jcol+iparity,:)'; end % jcol for jcol=M+1:Nmin1 Jacobian(:,jcol) = - psi_j(jcol-M,:)'; end % jcol Jacobian = 2*Jacobian; % symmetry correction if mod(N,2) == 1 % N odd, avoid double-counting the weight for x=0 Jacobian(:,M+1)=0.5 * Jacobian(:,M+1); end % if % end % if to compute Jacobian %disp('Jacobian') %Jacobian %Jacobiancheck % Jacdiff=Jacobian - Jacobiancheck rnorma(inewt)=max(abs(residual)); rnorm=rnorma(inewt); anorm=max(abs(xandw)); delta_a=Jacobian\residual; % delta_a is n-dim column vector dnorma(inewt)=max(abs(delta_a)); deltanorm=dnorma(inewt); minrnorm=1.E50; for ialph=0:nalpha % alphaline=(1+0.005*i)/(2^ialph); alphaline= 1 / (2^ialph); atemp1=alphaline* delta_a'; % n-dimensional column vectors atemp= xandw - atemp1; residual=feval(Res,atemp,c)'; rnorma(ialph+1)=max(abs(residual)); if rnorma(ialph+1) < minrnorm iminnorm=ialph; underrel=alphaline; minrnorm=rnorma(ialph+1); end % if end % ialpha loop disp([' inewt and stuff',' c=',num2str(c)]) inewt, minrnorm, underrel rhist(inewt)=minrnorm; gammaa(inewt)=underrel; rhist,gammaa atemp1=underrel* delta_a'; % atemp1 is n-dim column vector xandw=xandw - atemp1; iita(inewt)=inewt; Zinewt_underrel_minrnorm=[inewt underrel minrnorm]; semilogy(1:inewt,rhist,'g-') set(gca,'FontSize',18) title(['inewt=',int2str(inewt),' c=',num2str(c)]) ylabel('Residual') xlabel('Iteration number') drawnow end % Convergence test if block end % of Newton loop, inewt SHAR_EOF fi # end of overwriting check if test -f 'TP_eigprolODE.m' then echo shar: will not over-write existing file "'TP_eigprolODE.m'" else cat << "SHAR_EOF" > 'TP_eigprolODE.m' function [chi_array,B_array]=TP_eigprolODE(c,nmax); % Computes the prolate spheroidal functions of order zero. % These are normalized so that int_{-1}^{1} [\psi_{j}(x)]**2 dx=1 % Input: c=bandwidth % nmax=degree of highest prolate function desired. % [number of prolate functions computed is (nmax+1)] % output: chi_array contains the eigenvalues of the prolate % differential equation chi, sorted so the % smallest is first % B_array contains the NORMALIZED Legendre coefficients of the % NORMALIZED eigenfunctions. The j-th column is the % eigenfunction corresponding the chi_array(j), i. e., % psi_{j-1}(x; c). neig=max(2*nmax+2,30); N=nmax+1; % matrix element from Eq.(62), pg. 814, of Xiao, Rokhlin and Yarvin % Inverse Problems (2002) maindiag=ones(neig,1); sup=ones(neig,1); sub=ones(neig,1); for j=1:neig % k=2*j - 2; % for symmetric modes k=j-1; % for general case maindiag(j)=k*(k+1) + c*c*( 2*k*(k+1)-1) / ( (2*k+3)*(2*k-1) ) ; sup(j)= c*c*(k+1)*(k+2)/ ( (2*k+3)*sqrt( (2*k+1)*(2*k+5)) ); sub(j)=c*c*(k-1)*(k)/ ( (2*k-1)*sqrt( (2*k-3)*(2*k+1)) ); end % AA=zeros(neig,neig); for j=3:neig-2 AA(j,j)=maindiag(j); AA(j,j+2)=sup(j); AA(j,j-2)=sub(j); end % j AA(1,1)=maindiag(1); AA(1,3)=sup(1); AA(2,2)=maindiag(2); AA(2,4)=sup(2); AA(neig,neig)=maindiag(neig); AA(neig-2,neig)=sub(neig); AA(neig-1,neig-1)=maindiag(neig-1); AA(neig-3,neig-1)=sub(neig-1); [VV,DD]=eig(full(AA)); [chi_arrayfull,iindex]=sort(diag(DD)); chi_array=chi_arrayfull(1:N); B_array=zeros(neig,N); for j=1:N B_array_unnorm(:,j)=VV(:,iindex(j)); end % j % The coefficients are in terms of normalized Legendre functions, % but we now need to normalize the eigenfunctions so as to have % unit L2 norm on [-1, 1]. for j=1:N psi_squared(j) = B_array_unnorm(:,j)' * B_array_unnorm(:,j); end % j psinormfact = ones(1,N) ./ sqrt(psi_squared); for j=1:N B_array(:,j) = B_array_unnorm(:,j) * psinormfact(j); end % j SHAR_EOF fi # end of overwriting check cd .. cd .. cd .. # End of shell archive exit 0