% EXAMPLE3_35 Computations for Example 3.35. % clear eps0 Nmax M theta global eps0 Nmax M global theta % % Be sure to declare the variable theta global also in the routines % R_MOD and QUADRAT and make the changes indicated in the routines % QUADRAT and FEX35. % f1=' eta=%4.1f theta=%8.2e k=%2.1f eps0=%8.2e\n\n'; f2='%4.0f %4.0f %22.14e %12.4e %7.0f %7.0f\n'; fprintf('\n') disp(' M=2*N') %disp(' M=2*floor((N+1)/2)') %disp(' M=2') %disp(' M=0') Nm=10; eps0=100*eps; Nmax=200; eta=-1; theta=1e-4; %eta=-1; theta=1; k=.5; %k=1.5; %k=2.5; exact=.29051241701949266261; %exact=.46087845417799195532; %exact=1.1860735010757557836; %exact=1.3964418203491153395; %exact=2.6618732791071501381; %exact=7.6272560956534476326; % % The next values of exact relate to theta=1, theta=10, % and theta=100, respectively. % %exact=.383869768812140; %exact=.820188540209559; %exact=2.41632871284393; % % When theta is "large", increase eps0 to avoid an excessively % large value of Nmax that would be needed for convergence. Also % increase Nmax appropriately. For example, % eps0=1e6*eps; Nmax=800; % fprintf(f1,eta,theta,k,eps0) disp(' n m integral error Ncap kount') for N=1:Nm M=2*N; % M=2*floor((N+1)/2); % M=2; % M=0; [xw,Ncap,kount]=fermi_dirac(N,M,eta,theta,k,eps0,Nmax); int=sum(xw(:,2).*fex35(xw(:,1),eta,theta)); err=abs((int-exact)/exact); if ~all(N-[1 4 7 10]) fprintf(f2,N,M,int,err,Ncap,kount) end end