偏微分方程数值算法及matlab实验报告(6).docxVIP

偏微分方程数值实验报告六实验题目：用Ritz-Galerkin方法求边值问题的第n次近似，基函数为并用表格列出0.25，0.5,0.75三点处的真解和时的数值解。实验算法：将上述边值问题转化为基于虚功方程的变分问题为：求，满足其中记，引入的n维近似子空间利用Ritz-Galerkin公式：，则问题关于下的近似变分问题解中的系数满足程序代码：%主函数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(the solution of the un);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);function z=f(x) z=x.*x;endend%Galerkin方法function uh=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);for q=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: an input variable in a functionend%系数矩阵Afunction [lambda,weight] = quadpts1d(I,quadorder)numPts = ceil((quadorder+1)/2); ifnumPts 10numPts = 10; endswitchnumPtscase 1 A = [0 2.0000000000000000000000000];case 2 A = [0.5773502691896257645091488 1.0000000000000000000000000 -0.5773502691896257645091488 1.0000000000000000000000000];case 3 A = [ 0 0.8888888888888888888888889 0.7745966692414833770358531 0.5555555555555555555555556 -0.774596669241