有限元编程：基础篇

、前言

相信很多做过有限差分之后又想做做有限元的初学者会有和我一样的困惑，能看懂有限元算法的理论分析，但是真正应用到实际编程当中之前心里发怵，废话不多说，求人不如求己，看懂这篇文章将会让你迅速掌握有限元最基础的编程思想。

二、以经典扩散方程为例(反常扩散方程可类比此例)

考虑如下扩散方程初边值问题

⎧⎩⎨⎪⎪⎪⎪∂u(x,t)∂t=∂2u(x,t)∂x2+f(x,t),u(a,t)=0,u(b,t)=0,u(x,0)=φ(x),a<x<b,0<t≤T,0≤t≤T,a≤x≤b,(1)(1){∂u(x,t)∂t=∂2u(x,t)∂x2+f(x,t),a<x<b,0<t≤T,u(a,t)=0,u(b,t)=0,0≤t≤T,u(x,0)=φ(x),a≤x≤b,

首先我们需要对时间与空间区间进行剖分，其中M,NM,N分别记为空间与时间方向的剖分，h=b−aMh=b−aM, τ=TNτ=TN，h,τh,τ分别为空间与时间方向的步长，则a=x0<x1<⋯<xM=ba=x0<x1<⋯<xM=b, 0=t0<t1<⋯<tN=T0=t0<t1<⋯<tN=T。

空间方向利用有限元逼近， 可得（1）式的变分形式（弱形式）

(∂uh(t)∂t,vh)+a(uh,vh)=(f,vh),vn∈Xh(2)(2)(∂uh(t)∂t,vh)+a(uh,vh)=(f,vh),vn∈Xh

.

其中时间方向利用向后差分离散,可得求解问题（1）的全离散格式

(unh,vh)+τa(unh,vh)=(un−1h,vh)+τ(f(x,tn),vh)(3)(3)(uhn,vh)+τa(uhn,vh)=(uhn−1,vh)+τ(f(x,tn),vh)

其中, 记双线形式a(unh,vh)=∫Ω(unh)′(vh)′dxa(uhn,vh)=∫Ω(uhn)′(vh)′dx为刚度矩阵（stiffness matrix），(unh,vh)=∫Ωunhvhdx(uhn,vh)=∫Ωuhnvhdx为质量矩阵（mass matrix），(f(x,tn),vh)=∫Ωf(x,tn)vhdx(f(x,tn),vh)=∫Ωf(x,tn)vhdx 为载荷向量(load vector)。Ω=[a,b]=∪i=1,2,3,⋯M−1ΩiΩ=[a,b]=∪i=1,2,3,⋯M−1Ωi，Xh={ϕi}Xh={ϕi}为有限元空间，这里我们选用的是线性插值

ϕi=⎧⎩⎨⎪⎪⎪⎪x−xi−1h,xi+1 − xh,0,当x∈Ωi,当x∈Ωi+1,其他.(4)(4)ϕi={x−xi−1h,当x∈Ωi,xi+1 − xh,当x∈Ωi+1,0,其他.

因此，我们需要获得的问题（1）的有限元解unhuhn可表示为

unh=∑j=1M−1unjϕj,(5)(5)uhn=∑j=1M−1ujnϕj,

其中uniuin 为问题（1）在点(xi,tn)(xi,tn)处的数值解,将（5）式代入（3）式，我们可以获得求解问题（1）的有限元格式的代数方程组。.

(τA+B)un=Bun−1+τfˆn.(6)(6)(τA+B)un=Bun−1+τf^n.

其中un=[un1,un2,⋯,unM−1]Tun=[u1n,u2n,⋯,uM−1n]T, A=(aij)M−1×M−1A=(aij)M−1×M−1为刚度矩阵，由（4）式我们可以发现，对任意ii, ϕiϕi只在Ωi=[xi−1,xi]Ωi=[xi−1,xi]与Ωi+1=[xi,xi+1]Ωi+1=[xi,xi+1]两个区间不为零，其余区间即|i−j|≥2|i−j|≥2时aijaij均为零，则通过区间转换计算，我们发现

ai,i=∫Ω(ϕi)′(ϕi)′dx=∑i=1M∫Ωi(ϕi)′(ϕi)′dx=∫xixi−1(ϕi)′(ϕi)′dx+∫xi+1xi(ϕi)′(ϕi)′dx区间转换⟹∫101hdξ+∫101hdξ=2h,i=1,2,⋯,M−1.(7)(7)ai,i=∫Ω(ϕi)′(ϕi)′dx=∑i=1M∫Ωi(ϕi)′(ϕi)′dx=∫xi−1xi(ϕi)′(ϕi)′dx+∫xixi+1(ϕi)′(ϕi)′dx区间转换⟹∫011hdξ+∫011hdξ=2h,i=1,2,⋯,M−1.

ai−1,i=ai,i−1=∫Ω(ϕi)′(ϕi−1)′dx=∑i=1M∫Ωi(ϕi)′(ϕi−1)′dx=∫xi−1xi−1(ϕi)′(ϕi)′dx=∫10−1hdξ=−1h,i=2,3,⋯,M−1.(8)(8)ai−1,i=ai,i−1=∫Ω(ϕi)′(ϕi−1)′dx=∑i=1M∫Ωi(ϕi)′(ϕi−1)′dx=∫xi−1xi−1(ϕi)′(ϕi)′dx=∫01−1hdξ=−1h,i=2,3,⋯,M−1.

可知刚度矩阵AA为一三对角矩阵，同理可知质量矩阵B=(bij)M−1×M−1B=(bij)M−1×M−1也为三对角矩阵，其矩阵元素由如下计算获得

bi,i=∫Ω(ϕi)2dx=∑i=1M∫Ωi(ϕi)2dx=∫xixi−1(ϕi)2dx+∫xi+1xi(ϕi)2dx=h∫10ξ2dξ+h∫10(1−ξ)2dξ=2h3,i=1,2,⋯,M−1.(9)(9)bi,i=∫Ω(ϕi)2dx=∑i=1M∫Ωi(ϕi)2dx=∫xi−1xi(ϕi)2dx+∫xixi+1(ϕi)2dx=h∫01ξ2dξ+h∫01(1−ξ)2dξ=2h3,i=1,2,⋯,M−1.

bi−1,i=bi,i−1=∫Ω(ϕi)(ϕi−1)dx=∑i=1M∫Ωi(ϕi)(ϕi−1)dx=∫xi−1xi−1(ϕi)(ϕi)dx=h∫10ξ(1−ξ)dξ=h6,i=2,3,⋯,M−1.(10)(10)bi−1,i=bi,i−1=∫Ω(ϕi)(ϕi−1)dx=∑i=1M∫Ωi(ϕi)(ϕi−1)dx=∫xi−1xi−1(ϕi)(ϕi)dx=h∫01ξ(1−ξ)dξ=h6,i=2,3,⋯,M−1.

到了这里，我们终于可以松一口气了，主要的东西我们已经准备好了，，，，等等，，，还有载荷向量fˆn=[(f(x,tn),ϕ1),(f(x,tn),ϕ2),⋯,(f(x,tn),ϕM−1)]Tf^n=[(f(x,tn),ϕ1),(f(x,tn),ϕ2),⋯,(f(x,tn),ϕM−1)]T还没计算，接下来给出载荷向量的计算

(f(x,tn),ϕi)=∫Ωf(x,tn)ϕidx=∑i=1M−1∫Ωif(x,tn)ϕidx=∫xixi−1f(x,tn)ϕidx+∫xi+1xif(x,tn)ϕidx=h∫10f(xi−1+hξ,tn)xdx+h∫10f(xi+hξ,tn)(1−ξ)dξ.(11)(11)(f(x,tn),ϕi)=∫Ωf(x,tn)ϕidx=∑i=1M−1∫Ωif(x,tn)ϕidx=∫xi−1xif(x,tn)ϕidx+∫xixi+1f(x,tn)ϕidx=h∫01f(xi−1+hξ,tn)xdx+h∫01f(xi+hξ,tn)(1−ξ)dξ.

为计算以上的一重积分，我在程序中采用的是17点GaussGauss求积公式，绝对符合我们对算法精度的要求！不会GaussGauss积分公式的同学可以查阅文献（黄云清等，数值计算方法，科学出版社，p.116）.具体程序实现如下

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26

function u = single_quad( a,b,f )