function ab=abRys_Jn_with_pts(c,m,nint,npts) %Guass Stieltjes procedure for caluation of recurrence coefficients % alpha and betta for weight function w(x)=x*x*exp(-x*x) format long e xmin=-1;xmax=1; %nint=# of intervals; npts=# of pts per interval [xmin,xmax] pwmd = mdquad(nint,npts,xmin,xmax); ntot=nint*npts; %multi-domain matrix of quad pts p and wts w: %weight function w(x)=x*x*exp(-x*x) cc=c*ones(ntot,1); xm=[1:1:m]; p=pwmd(:,1); w=pwmd(:,2); psq=p.*p; wtfcn=exp(-cc.*psq); srwtfcn=exp(-cc.*psq/2); %wtfcn=(log(p)).^2; wtfcn=exp(-cc./p); % Polynomial Q_0(x) b0=ones(ntot,1); % Calculate recurrence coefficients % First two integrals s1=sum(w.*wtfcn); s2=sum(w.*(p.*wtfcn));mom1=s1; h(1)=s1; alfa(1)=s2/s1; beta(1)=0;k=1; %Norm and alpha_1 and betta_1 myfile1=fopen('ab.dat','wt'); fprintf(myfile1,'%10.8f %10.8f\n',alfa(k),beta(k)); %Polynomial P_1(x) b1=p-alfa(1); %Next two integrals s1=sum(w.*(wtfcn.*(b1.^2))); s2=sum(w.*(p.*(wtfcn.*(b1.^2)))); alfa(2)=s2/s1; h(2)=s1; beta(2)=h(2)/h(1);k=2; %alpha_2, norm and beta_2 fprintf(myfile1,'%10.8f %10.8f\n',alfa(k),beta(k)); % % Recurrence for Q_n(x) - not normalized to unity % for k=3:m pma=p-alfa(k-1); %Recurrence for the next polynomial b2=pma.*b1-beta(k-1)*b0; s1=sum(w.*(wtfcn.*(b2.^2))); s2=sum(w.*(p.*(wtfcn.*(b2.^2)))); % alfa_kl, norm and beta_k alfa(k)=s2/s1; h(k)=s1; beta(k)=h(k)/h(k-1); fprintf(myfile1,'%10.8f %10.8f\n',alfa(k),beta(k)); b0=b1; b1=b2; b3=pma.*b1-beta(k-1)*b0; end fclose(myfile1); zero=0*ones(1,m); plot(p,b3.*srwtfcn,'-k','linewidth',1.6) hold on b=beta(2:m); rtb=sqrt(b); t=diag(rtb,-1)+diag(alfa)+diag(rtb,1); [f,lambda]=eig(t); pt=diag(lambda); plot(pt,zero,'ok','markersize',7,'markerfacecolor','k') hold on axis([-1 1 -2e-6 2e-6]) set(gca,'FontSize',20) set(gca,'Ytick',[-2e-6:0.5e-6:2e-6],'linewidth',1.6) set(gca,'Xtick',[-1:.5:1],'linewidth',1.6) %axis([-3 3 -1 1.2]) xlabel('$x$','Interpreter','LaTex','FontSize',24) ylabel('$\sqrt{e^{-cx^2}}J_{20}(x,c)$','Interpreter','LaTex','FontSize',24) str1 = {'$c = 2$'}; text(0,1.5e-6,str1,'Interpreter','Latex','Fontsize',26) hold off for i=1:m wt(i)=mom1*f(1,i)^2; end ptwt=[pt,wt']; myfile2=fopen('ptswts.dat', 'wt'); for i=1:m fprintf(myfile2,'%22.14e %22.14e\n',pt(i),wt(i)); end fclose(myfile2); function pwmd = mdquad(nint,npts,xmin,xmax) format long e ntot=nint*npts; %Quadrature for mutlidomain with nint intervals and npts per interval dx=(xmax-xmin)/ntot; for i=1:nint a=xmin+(i-1)*npts*dx; b=a+npts*dx; pw=fejer2(a,b,npts); %Fejer quadrature used for each interval if i==1 pwmd=pw; else pwmd=cat(1,pwmd,pw); end end % FEJER Fejer quadrature rule. % The call pw=fejer(n) generates the n-point Fejer quadrature rule. function pw=fejer2(a,b,N) format long e n=N:-1:1; m=1:floor(N/2); th=(2*n-1)*pi./(2*N); p=cos(th'); for k=N:-1:1 s=sum(cos(2*m*th(k))./(4*(m.^2)-1)); w(k)=2*(1-2*s)/N; end r1=(b-a)/2.; r2=(a+b)/2.; ps=r1*p+r2; ws=r1*w; %Map [a,b] onto [-1,1] pw=[ps,ws'];