偏微分方程的有限元法求解_高等教育-微积分.pdf
精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载%对于 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;%d1/d dN2dxi=1/2;%d2/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=Kb;%v 为 Kx=b 的解 plot(x,v,*-);%画图并比较 hold on;plot(xx,uex);hold off;xlabel(x);ylabel(u);
