偏微分方程的有限元法求解.doc

%对于d2u/dx2=f的FEM解算器，其中f=x*(1-x) %%边界条件u(0)=0, u(1)=0. %精确解用以比对 xx=linspace(0,1,101);%产生0-1之间的均分指令，101为元素个数 uex=(1/6).*xx.^3-(1/12).*xx.^4-(1/12).*xx; %对力项设置高斯点的数目 NGf=2; if (NGf==2) xiGf=[-1/sqrt(3);1/sqrt(3)];%ξ1、ξ2的值 aGf=[1 1]; else, NGf=1; xiGf=[0.0]; aGf=[2.0]; end %单元数目 Ne=5; %建立网格节点 x=linspace(0,1,Ne+1); %零刚性矩阵 K=zeros(Ne+1,Ne+1); b=zeros(Ne+1,1); %对所有单元循环计算刚性和残差 for ii=1:Ne, kn1=ii; kn2=ii+1; x1=x(kn1); x2=x(kn2); dx=x2-x1;%每一个单元的长度 dxidx=2/dx;%dξ/dx dxdxi=1/dxidx;%dx/dξ dN1dxi=-1/2;%dζ1/dξ dN2dxi=1/2;%dζ2/dξ dN1dx=dN1dxi*dxidx;%-1/(xj-xj-1) dN2dx=dN2dxi*dxidx;%1/(xj-xj-1) K(kn1,kn1)=K(kn1,kn1)-2*dN1dx*dN1dx*dxdxi;%Rj的第二项 K(kn1,kn2)=K(kn1,kn2)-2*dN1dx*dN2dx*dxdxi; K(kn2,kn1)=K(kn2,kn1)-2*dN2dx*dN1dx*dxdxi; K(kn2,kn2)=K(kn2,kn2)-2*dN2dx*dN2dx*dxdxi; %用高斯积分估计力项的积分 for nn=1:NGf%NGf=2 xiG=xiGf(nn);%得到高斯点的ξ N1=0.5*(1-xiG);%求N1和N2(即在xiG的权重/插值) 形状函数在ξ的值 N2=0.5*(1+xiG);%ζ值 fG=xiG*(1-xiG);%对ξ点求f gG1=N1*fG*dxdxi;%在节点处估计权函数在高斯点的被积函数 gG2=N2*fG*dxdxi;%估计是个积分值 b(kn1)=b(kn1)+aGf(nn)*gG1;% aGf为1 b(kn2)=b(kn1)+aGf(nn)*gG2; end end %在x=0处设置Dirichlet条件 kn1=1; K(kn1,:)=zeros(size(1,Ne+1)); K(kn1,kn1)=1; b(kn1)=0; %在x=1处设置Dirichlet条件 kn1=1; K(kn1,:)=zeros(size(1,Ne+1)); K(kn1,kn1)=1; b(kn1)=0; %求解方程 v=K\b;%v为Kx=b的解 plot(x,v,*-);%画图并比较 hold on; plot(xx,uex); hold off; xlabel(x); ylabel(u);