% EXAMPLE2_29 Numerical illustration for Example 2.29. % % See also W. Gautschi, ``Algorithm 726: ORTHPOL -- a package of routines % for generating orthogonal polynomials and Gauss-type quadrature rules'', % ACM Trans. Math. Software 20 (1994), 21-62, Example 3.1. % f1='%9.3f %4.0f %20.15f %9.2e\n'; f2='%14.0f %20.15f %9.2e\n'; fprintf('\n') disp(' Recursion coefficients for an elliptic weight') disp(' function computed via Chebyshev moments') fprintf('\n') disp(' k2 n beta magb') load -ascii coeffell; bexact=coeffell; N=80; abm=r_jacobi(2*N-1,-1/2); for k2=[.1 .5 .9 .999] mom=mm_ell(N,k2); [ab,normsq]=chebyshev(N,mom,abm); if k2==.1 magb=abs(ab(:,2)-bexact(1:N))/eps; n=[1 2 6 12]; for k=1:4 if k==1 fprintf(f1,k2,n(k)-1,ab(n(k),2),magb(n(k))) else fprintf(f2,n(k)-1,ab(n(k),2),magb(n(k))) end end fprintf('\n') elseif k2==.5 magb=abs(ab(:,2)-bexact(N+1:2*N))/eps; n=[1 2 9 21]; for k=1:4 if k==1 fprintf(f1,k2,n(k)-1,ab(n(k),2),magb(n(k))) else fprintf(f2,n(k)-1,ab(n(k),2),magb(n(k))) end end fprintf('\n') elseif k2==.9 magb=abs(ab(:,2)-bexact(2*N+1:3*N))/eps; n=[1 2 20 44]; for k=1:4 if k==1 fprintf(f1,k2,n(k)-1,ab(n(k),2),magb(n(k))) else fprintf(f2,n(k)-1,ab(n(k),2),magb(n(k))) end end fprintf('\n') else magb=abs(ab(:,2)-bexact(3*N+1:4*N))/eps; n=[1 2 20 44 80]; for k=1:5 if k==1 fprintf(f1,k2,n(k)-1,ab(n(k),2),magb(n(k))) else fprintf(f2,n(k)-1,ab(n(k),2),magb(n(k))) end end end end