function legendre_norm_expmap % This calculates the norm or orthogonality of two Maxwell polynomials format long e hold off % alpha_n and beta_n for the Maxwell polynomials, p = 2. load abspd.dat; n=90; alf=abspd(:,1); b=abspd(:,2); rtb=sqrt(b); xn=[]; acc1=[]; acc2=[]; % Start NQUAD LOOP m=2; for nquad=5:2:60 xn=[xn nquad]; rtbx=rtb; rtbx(nquad:n-1,:)=[]; %delete the bottom nquad to n rows alfx=alf; alfx(nquad+1:n,:)=[]; %Legendre quadrature points and weights ptwt=legptwt(nquad); pt=ptwt(:,1); wt=ptwt(:,2); rmax=8; % Boyd1 exponential map r=-log((1-pt)/2); wtfcn=(r.*r).*exp(-r.*r); srwtfcn=sqrt(wtfcn); % Calculation of the Maxwell Polynomials by recurrence gam=sqrt(pi/4); srgam=sqrt(gam); rtpi=sqrt(pi); b1=sqrt(4/(rtpi*(3*pi/8-1))); a1=rtpi*b1/2; p1=(2/sqrt(rtpi))*ones(1,nquad); p2=a1*r'-b1*ones(1,nquad); q1=p2; for n=2:6 p3=(r'.*p2-alf(n)*p2)/rtb(n) - p1*rtb(n-1)/rtb(n); p1=p2; p2=p3; if n==4 p4=p2; else end end % End Calcn of Maxwell Polynomials npt=length(pt); zero=0*ones(1,nquad); %Boyd1 mapping maxpolyn1=0.5*exp(r').*wtfcn'.*p4.*p4; maxpolyn2=0.5*exp(r').*wtfcn'.*p2.*p2; integ1=sum(wt'.*maxpolyn1); integ2=sum(wt'.*maxpolyn2); err1=log10(abs(1-integ1)); err2=log10(abs(1-integ2)); acc1=[acc1 err1]; acc2=[acc2 err2]; plot(xn,acc1,'-^k',xn,acc2,'-ok','linewidth',1.6,'markerfacecolor','k','markersize',10); hold on axis([10 60 -10 0]) set(gca,'FontSize',32) set(gca,'Ytick',[-10:2:0],'linewidth',1.6) set(gca,'Xtick',[10:10:60],'linewidth',1.6) xlabel('$N$','Interpreter','LaTex','FontSize',32) ylabel('$\log_{10}|1-{\rm I}_{approx}|)$','Interpreter','LaTex','FontSize',32) %Boyd1 mapping str1={'Exponential map'}; str2={'${x=1-2e^{-r}}$'}; text(12,-8,str1,'Interpreter','latex','fontsize',32) text(12,-9,str2,'Interpreter','latex','fontsize',32) set(gcf, 'Units','centimeters','Papersize',[36,36]) set(gcf, 'Units','centimeters','Position',[3 3 24 20]) end % Legendre quadrature points and weights function ptwt = legptwt(n) xn=[0:1:n-1]; xnsq=xn.*xn; b=xnsq./(4*xnsq-1); rtb=sqrt(b); rtb(1)=[]; t=diag(rtb,-1)+diag(rtb,1); [f,lambda]=eig(t); pt=diag(lambda); wt=2*f(1,:).^2; ptwt=[pt,wt'];