Ritz-Galerkin法求解边值问题

Ritz-Galerkin法求解边值问题2010-05-3013:37functionvarargout=saxplliu(varargin)%function%-u''+ru=f(x)0x2*pi%u(0)=0,u(2*pi)=0n=8;r=1;%A1=ones(n,n);%A2=ones(n,n);%b=ones(n,1);[x1,w1]=lglnodes(n)fori=1:nA1=i^2*pi/4;A2=pi;s=0;forj=1:n+1s=s+f(pi*(x1(j)+1))*fG(pi*(x1(j)+1),i)*w1(j);endb=pi*s;c(i)=b/(A1+r*A2);end%A%bh=0.2;x=0:h:2*pi;n1=length(x);y1=x;y2=x;fori=1:n1y1(i)=u(x(i));endfori=1:n1y20=0;forj=1:ny20=y20+c(j)*fG(x(i),j);endy2(i)=y20;endfigure(1)plot(x,y1,x,y2,'ro');legend('精确解','数值解')figure(2)plot(x,abs(y1-y2));title('误差图')function[x,w]=lglnodes(N)N1=N+1;x=cos(pi*(0:N)/N)';P=zeros(N1,N1);xold=2;whilemax(abs(x-xold))epsxold=x;P(:,1)=1;P(:,2)=x;fork=2:NP(:,k+1)=((2*k-1)*x.*P(:,k)-(k-1)*P(:,k-1))/k;endx=xold-(x.*P(:,N1)-P(:,N))./(N1*P(:,N1));endw=2./(N*N1*P(:,N1).^2);functiony=fG(x,I)y=sin(I*x/2)functiony=f(x)y=(1/pi^2+1)*2*sin(x/pi)/sin(2)-x/pi;functiony=u(x)y=2*sin(x/pi)/sin(2)-x/pi