偏微分方程的有限元法求解
- 1、下载文档前请自行甄别文档内容的完整性,平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益,请联系客服反馈,我们会尽快为您处理(人工客服工作时间:9:00-18:30)。
%对于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;%dζ1/dξ
dN2dxi=1/2;%dζ2/dξ
dN1dx=dN1dxi*dxidx;%-1/(x j-x j-1)
dN2dx=dN2dxi*dxidx;%1/(x j-x j-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=K\b;%v为Kx=b的解
plot(x,v,'*-');%画图并比较
hold on;
plot(xx,uex);
hold off;
xlabel('x');
ylabel('u');