《偏微分方程的有限元法求解_高等教育-微积分.pdf》由会员分享，可在线阅读，更多相关《偏微分方程的有限元法求解_高等教育-微积分.pdf（22页珍藏版）》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载 精品资料 欢迎下载%对于 d2u/dx2=f 的 FEM 解算器，其中 f=x*(1-x)%边界条件 u(0)=0,u(1)=0.%精确解用以比对 xx=linspace(0,1,101);%产生 0
2、-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
3、(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,
4、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);