function radint_maxp2_vsN_b % Maxwell (p = 2) quadrature for the radial function, Eq. (3.10) % with f(r) given by the sum of 3 Gaussians, Eq. (3.15) - Fig. 3.6(B) format long e hold off myfile = fopen('radint_sp.dat', 'wt'); exact=(1.d0+1.d0/sqrt(10.d0)+1.d0/sqrt(100.d0))*sqrt(pi)/4.d0; %exact=(1.d0+1.d0/sqrt(10.d0))*sqrt(pi)/4.d0; % nucrxspeed2 is versus the scale factor for 3 values of N - fixed b yerr=[]; int=[]; str1='-ok'; str2='-^k'; str3='-sk';str4='-dk'; symbol=[str1 str2 str3 str4]; % Scale factors %sc=[.5 .6 .7 .8]; sc=[.3 .4 .45 .5]; kmax=0; %This is the small loop through 4 scales yns=[]; for k=1:1:4 kmin=kmax+1; kmax=kmin+2; nx=[]; int=[]; yerr=[]; scale=sc(k); % This is the larger loop through N % for m=4:1:16 for m=4:1:24 nx=[nx m]; % fprintf(myfile,'%8.2f\n',scale); ptwt=abmaxwellp2(m); x=ptwt(:,1); w=ptwt(:,2); n1=length(x); s=0.; bigw=w./(x.^2.*exp(-x.*x)); ws=scale*bigw; xs=scale*x; fx=xs.^2.*(exp(-xs.^2)+10*exp(-10*xs.^2)+100*exp(-100*xs.^2)); %fx=xs.^2.*(exp(-xs.^2)+10*exp(-10*xs.^2)); s=sum(ws.*fx); int=[int s]; exact=(1.d0+1.d0/sqrt(10.d0)+1.d0/sqrt(100.d0))*sqrt(pi)/4.d0; %exact=(1.d0+1.d0/sqrt(10.d0))*sqrt(pi)/4.d0; err=log10(abs(s-exact)/exact); fprintf('%i %16.8f\n',m,s-exact) yerr=[yerr err]; end pcoeff=polyfit(nx,yerr,1) for kkk=1:2 fprintf(myfile,'%13.5f %13.5\n',pcoeff(kkk)); end pause yn=polyval(pcoeff,nx); yns=[yns,yn']; plot(nx,yerr,symbol(kmin:kmax),'markersize',10,'markerfacecolor','k') hold on xlabel('${\rm N}$','Interpreter','latex','fontsize',32) ylabel('$\log_{10}[{\rm Relative}\;\; {\rm Error}]$','Interpreter','latex','fontsize',32) axis([4 24 -7 0]) set(gca,'FontSize',36) set(gca,'Ytick',[-7:2:0],'linewidth',1.6) set(gca,'Xtick',[4:4:24],'linewidth',1.6) end plot(nx,line,'--k','linewidth',1.6) %str={'N = 20'}; %text(0.6,.55,str,'Interpreter','latex','fontsize',24) %title('fx=x^2*(exp(-x^2)+10*exp(-10*x^2)+100*exp(-100*x^2)','fontsize',14) legend('s=0.3','s=0.4','s=0.45','s=0.5','Location','SouthWest') %str={'(A)'}; str={'(B)'}; text(20,-1,str,'Interpreter','Latex','Fontsize',36) plot(nx,yns(:,1),'--k','linewidth',1.6) plot(nx,yns(:,2),'--k','linewidth',1.6) plot(nx,yns(:,3),'--k','linewidth',1.6) plot(nx,yns(:,4),'--k','linewidth',1.6) set(gcf, 'Units','centimeters','Papersize',[36,36]) set(gcf, 'Units','centimeters','Position',[3 3 24 20]) % Maxwell points and weights for p = 2 function ptwt = abmaxwellp2(m) load abspd.dat; n=90; a=abspd(1:m,1); b=abspd(1:m,2); rtb=sqrt(b); rtb(m)=[]; t=diag(rtb,-1)+diag(a)+diag(rtb,1); [f,lambda]=eig(t); pt=diag(lambda); for i=1:m wt(i)=sqrt(pi)*f(1,i)^2/4.d0; end ptwt=[pt,wt'];