偏微分方程地有限元法求解Word文档下载推荐.docx

偏微分方程地有限元法求解Word文档下载推荐.docx

《偏微分方程地有限元法求解Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《偏微分方程地有限元法求解Word文档下载推荐.docx（22页珍藏版）》请在冰豆网上搜索。

%对力项设置高斯点的数目

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ξ/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;

%用高斯积分估计力项的积分

fornn=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;

%在x=0处设置Dirichlet条件

kn1=1;

K（kn1,:

）=zeros（size（1,Ne+1））;

K（kn1,kn1）=1;

b（kn1）=0;

%在x=1处设置Dirichlet条件

%求解方程

v=K\b;

%v为Kx=b的解

plot（x,v,'

*-'

）;

%画图并比较

holdon;

plot（xx,uex）;

holdoff;

xlabel（'

x'

ylabel（'

u'