偏微分方程的有限元法求解

%对于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');