function lindh_integ_maxp2B % This is the calculation of the Lindh Gaussian - Eq. (5) with speed % quadrature points. The convergence is shown in Fig. 3.5 format long e hold off % Strings for plotting yerr=[]; str1='-ok'; str2='-sk'; str3='-^k'; str4='-dk'; alfa=4; scale=1/sqrt(alfa); symbol=[str1 str2 str3 str4]; % Specific ell values strell=[6 8 10 14]; xm=[]; kmax=0; for k=1:1:4 % Scan through ell values ell=strell(k); xm=[]; yerr=[]; for m=2:1:10 % Scan through (low) quadrature orders xm=[xm m]; ptwt=maxwellp2(m); x=ptwt(:,1); w=ptwt(:,2); n1=length(x); bigw=w./(x.^2.*exp(-x.*x)); ws=scale*bigw; xs=scale*x; fx=xs.^(ell+2).*exp(-alfa*xs.^2); %fx=10*xs.^2.*exp(-10*xs.^2); s=2*alfa^((ell+3)/2)/gamma((ell+3)/2)*sum(ws.*fx); err=log10(abs(s-1.d0)/1.d0); yerr=[yerr err]; end kmin=kmax+1; kmax=kmin+2; plot(xm,yerr,symbol(kmin:kmax),'linewidth',1.6,'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([2 10 -20 0]) set(gca,'FontSize',32) set(gca,'Ytick',[-20:5:0],'linewidth',1.6) set(gca,'Xtick',[2:2:10],'linewidth',1.6) end leg=legend('$\ell=6$','$\ell=8$','$\ell=10$','$\ell=14$','Location','NorthEast'); set(leg,'Interpreter','latex','fontsize',32); str={'$\alpha = 4; s=1/\sqrt{\alpha} = 0.5$'}; text(3,-18,str,'Interpreter','latex','fontsize',32) set(gcf, 'Units','centimeters','Papersize',[36,36]) set(gcf, 'Units','centimeters','Position',[3 3 24 20]) % Speed points and weights for p = 1 function ptwt = maxwellp2(m) load abmaxp2.dat; n=100; a=abmaxp2(1:m,1); b=abmaxp2(1:m,2); rtb=sqrt(b); rtb(m)=[]; t=diag(rtb,-1)+diag(a)+diag(rtb,1); [f,lambda]=eig(t); pt=diag(lambda); wt=sqrt(pi)*f(1,:).^2/4.d0; ptwt=[pt,wt'];