% EXAMPLE3_65 Computations for Example 3.65. % f1='%10.0f %20.15f %12.4e %12.4e\n'; f2=' timing=%6.3f\n'; fprintf('\n') Sexact=1-log(2); % % Gauss-Laguerre % N=50; disp(' n S err can ') ab=r_laguerre(N); tic for n=1:N qtmax=0; xw=gauss(n,ab); gl=0; for k=1:n t=xw(k,1); qterm=xw(k,2)/(exp(t)+1); if abs(qterm)>qtmax, qtmax=abs(qterm); end gl=gl+xw(k,2)/(exp(t)+1); end can=qtmax/abs(gl); err=abs((gl-Sexact)/Sexact); fprintf(f1,n,gl,err,can) end time=toc; fprintf(f2,time) fprintf('\n') % % Rational Gauss-Laguerre % global eps0 Nmax M mc mp iq idelta irout AB Z ab0 eps0=100*eps; Nmax=200; ab0=r_laguerre(Nmax); mc=1; mp=0; iq=1; idelta=2; irout=1; AB=[0 Inf]; N=10; tic for n=1:N M=2*floor((n+1)/2); for m=1:2:M-1 Z(m,1)=-2/((m+1)*pi*i); Z(m,2)=1; Z(m+1,1)=2/((m+1)*pi*i); Z(m+1,2)=1; end [abmod,Ncap,kount]=r_mod(n,ab0); xw=gauss_rational(n,abmod); rgl=0; qtm=0; for k=1:n t=xw(k,1); qterm=xw(k,2)/(exp(t)+1); if abs(qterm)>qtm, qtm=abs(qterm); end rgl=rgl+xw(k,2)/(exp(t)+1); end can=qtm/abs(rgl); err=abs((rgl-Sexact)/Sexact); fprintf(f1,n,rgl,err,can) end time=toc; fprintf(f2,time) fprintf('\n') % % Gauss-Fermi % N=20; load -ascii abfermi1; ab=abfermi1(1:N,:); tic for n=1:N qtmax=0; xw=gauss(n,ab); ge=0; for k=1:n t=xw(k,1); qterm=xw(k,2)*exp(-t); if abs(qterm)>qtmax, qtmax=abs(qterm); end ge=ge+xw(k,2)*exp(-t); end can=qtmax/abs(ge); err=abs((ge-Sexact)/Sexact); fprintf(f1,n,ge,err,can) end time=toc; fprintf(f2,time)