偏微分方程的有限元法求解_精品文档.doc

偏微分方程的有限元法求解_精品文档.doc

《偏微分方程的有限元法求解_精品文档.doc》由会员分享，可在线阅读，更多相关《偏微分方程的有限元法求解_精品文档.doc（22页珍藏版）》请在冰豆网上搜索。

%对于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ξ/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;

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=K\b;%v为Kx=b的解

plot（x,v,'*-'）;%画图并比较

holdon;

plot（xx,uex）;

holdoff;

xlabel（'x'）;

ylabel（'u'）;