偏微分方程数值实验报告六实验题目：用 Ritz-Galerkin 方法求边值问题—u+u 二 X2u(0)二 0,u⑴ 二 1的第 n 次近似 u(x)，基函数为©(x)二 sin(inx),i 二 1,2,...n 并用表格列出ni0.25,0.5,0.75 三点处的真解和 n 二 1,2,3,4 时的数值解。实验算法：将上述边值问题转化为基于虚功方程的变分问题为：uGH1(I)，满足 a(u,v)=—(x,x),VveH1(I)00其中a(u,v)=(Lu,v)=f1(u''v+uv)dx=f1(—u'v'+uv)dx00记 w(x)=sin(兀 x)，引入 H1(I)的 n 维近似子空间 U={©,©，…,©}0n12n©=sin(inx),i=1,2,n 利用 Ritz-Galerkin 公式：i工 a(©,©)=(f,©),j=1,2,…,n，则问题关于 U 下的近似变分问题解 jijnj=1u(x)=工 c©(x)中的系数 c=(c,c,c)TGRn满足工 a(©,©)=(x2,©)ii12njiji=1j=1程序代码：%主函数 pde=modeldata();I=[0,1];%积分精度 quadorder=10;n=[1,2,3,4];fori=1:length(n)uh{i}=Galerkin(pde,I,n(i),quadorder);endshowsolution(uh{1},'-bo');holdonshowsolution(uh{2},'-rx');holdonshowsolution(uh{3},'-.ko');holdonshowsolution(uh{4},'--k*');holdofftitle('thesolutionoftheun');xlabel('x');ylabel('u');legend('n=1','n=2','n=3','n=4');%计算这三点的数值解x=[1/4,1/2,3/4];un=zeros(length(n),3);fori=1:length(n)[v,]=basis(x,n(i));un(i,:)=(v'*uh{i})';endformatshortedisp('un');disp(un);%存储数据functionpde=modeldata()pde=struct('f',@f);functionz=f(x)z=x.*x;endend%Galerkin 方法functionuh=Galerkin(pde,I,n,quadorder)h=I(2)-I(1);[lambda,weight]=quadpts1d(I,quadorder);nQuad=length(weight);A=zeros(n,n);b=zeros(n,1);forq=1:nQuadgx=lambda(q);w=weight(q);x=[0.25;0.5;0.75];[phi,gradPhi]=basis(gx,n);A=A+(-gradPhi*gradPhi'+phi*phi')*w;b=b+pde.f(gx)*phi*w;endA=h*A;b=h*b;uh=A\b;end%构造基函数function[phi,gradPhi]=basis(x,n)m=length(x);w=sin(pi*x);v=ones(n,m);v(2:end,:)=bsxfun(@times,v(2:end,:),x);v=cumprod(v,1);phi=bsxfun(@times,v,w);gw=pi*cos(pi*x);gv=[zeros(1,m);v(1:end-1,:)];gv(3:end,:)=bsxfun(@times,(2:n-1)',gv(3:end,:));gradPhi=bsxfun(@times,v,gw)+bsxfun(@times,gv,w);end%数值解的图像functionshowsolution(u,varargin)x=0:0.01:1;n=length(u);[v,]=basis(x,n);y=v'*u;plot(x,y,varargin{:});%varargin:aninputvariableinafunctionend%...
