如果您无法下载资料，请参考说明：
1、部分资料下载需要金币，请确保您的账户上有足够的金币
2、已购买过的文档，再次下载不重复扣费
3、资料包下载后请先用软件解压，在使用对应软件打开
%对于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=[11];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);%对所有单元循环计算刚性和残差forii=1:Ne,kn1=ii;kn2=ii+1;x1=x(kn1);x2=x(kn2);dx=x2-x1;%每一个单元的长度dxidx=2/dx;%dξ/dxdxdxi=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;%用高斯积分估计力项的积分fornn=1:NGf%NGf=2xiG=xiGf(nn);%得到高斯点的ξN1=0.5*(1-xiG);%求N1和N2(即在xiG的权重/插值)形状函数在ξ的值N2=0.5*(1+xiG);%ζ值fG=xiG*(1-xiG);%对ξ点求fgG1=N1*fG*dxdxi;%在节点处估计权函数在高斯点的被积函数gG2=N2*fG*dxdxi;%估计是个积分值b(kn1)=b(kn1)+aGf(nn)*gG1;%aGf为1b(kn2)=b(kn1)+aGf(nn)*gG2;endend%在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,'*-');%画图并比较holdon;plot(xx,uex);holdoff;xlabel('x');ylabel('u');
立即下载
