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

-%对于 d u/d*=f 的 FEM 解算器，其中 f=*(1-*)%边界条件 u(0)=0,u(1)=0.%精确解用以比对*=linspace(0,1,101);%产生 0-1 之间的均分指令，101 为元素个数ue*=(1/6).*.3-(1/12).*.4-(1/12).*;%对力项设置高斯点的数目NGf=2;if(NGf=2)*iGf=-1/sqrt(3);1/sqrt(3);%1、2的值aGf=1 1;else,NGf=1;*iGf=0.0;aGf=2.0;end%单元数目Ne=5;%建立网格节点*=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;*1=*(kn1);*2=*(kn2);d*=*2-*1;%每一个单元的长度d*id*=2/d*;%d/d*d*d*i=1/d*id*;%d*/ddN1d*i=-1/2;%d1/ddN2d*i=1/2;%d2/ddN1d*=dN1d*i*d*id*;%-1/(*j-*j-1)dN2d*=dN2d*i*d*id*;%1/(*j-*j-1)K(kn1,kn1)=K(kn1,kn1)-2*dN1d*dN1d*d*d*i;%Rj的第二项K(kn1,kn2)=K(kn1,kn2)-2*dN1d*dN2d*d*d*i;K(kn2,kn1)=K(kn2,kn1)-2*dN2d*dN1d*d*d*i;K(kn2,kn2)=K(kn2,kn2)-2*dN2d*dN2d*d*d*i;%用高斯积分估计力项的积分for nn=1:NGf%NGf=2*iG=*iGf(nn);%得到高斯点的N1=0.5*(1-*iG);%求 N1 和 N2(即在*iG 的权重/插值)形状函数在的值.z.22-N2=0.5*(1+*iG);%值fG=*iG*(1-*iG);%对点求 fgG1=N1*fG*d*d*i;%在节点处估计权函数在高斯点的被积函数gG2=N2*fG*d*d*i;%估计是个积分值b(kn1)=b(kn1)+aGf(nn)*gG1;%aGf为 1b(kn2)=b(kn1)+aGf(nn)*gG2;endend%在*=0 处设置 Dirichlet 条件kn1=1;K(kn1,:)=zeros(size(1,Ne+1);K(kn1,kn1)=1;b(kn1)=0;%在*=1 处设置 Dirichlet 条件kn1=1;K(kn1,:)=zeros(size(1,Ne+1);K(kn1,kn1)=1;b(kn1)=0;%求解方程v=Kb;%v 为 K*=b 的解plot(*,v,*-);%画图并比较hold on;plot(*,ue*);hold off;*label(*);ylabel(u);.z.