本文将介绍矩形区域上Poisson方程 − Δ u = f , Δ = ∂ ∂ x 2 + ∂ ∂ y 2 -\Delta u=f,\Delta= \frac{\partial}{\partial x^2}+\frac{\partial}{\partial y^2} −Δu=f,Δ=∂x2∂​+∂y2∂​的有限元程序编写，包括主要包括三角网格剖分，基函数构造以及函数在任意三角区域二重积分的计算。进一步构造基函数的相关偏导数，荷载向量，刚度矩阵以及替换边界值。

为了使过程清晰，每一部分的详细单独写了分析过程。

对于在区域 Ω \Omega Ω上Poisson问题： − Δ u = f , Δ = ∂ ∂ x 2 + ∂ ∂ y 2 -\Delta u=f,\Delta= \frac{\partial}{\partial x^2}+\frac{\partial}{\partial y^2} −Δu=f,Δ=∂x2∂​+∂y2∂​其中 u u u满足如下边值条件： u ∣ ∂ Ω = α ( x , y ) u|_{\partial \Omega}=\alpha(x,y) u∣∂Ω​=α(x,y)，它等价的变分形式为：求 u ∈ H 1 ( Ω ) , u ∣ ∂ Ω = α ( x , y ) u\in H^1(\Omega),u|_{\partial \Omega}=\alpha(x,y) u∈H1(Ω),u∣∂Ω​=α(x,y)，使 a ( u , v ) = ( f , v ) a(u,v)=(f,v) a(u,v)=(f,v)，对任意 v ∈ H 0 1 ( Ω ) v\in H_0^1(\Omega) v∈H01​(Ω)成立，其中 a ( u , v ) = ∬ Ω ( ∂ u ∂ x ∂ v ∂ x + ∂ u ∂ y ∂ v ∂ y ) d x d y a(u,v) = \iint_{\Omega}(\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial y})dxdy a(u,v)=∬Ω​(∂x∂u​∂x∂v​+∂y∂u​∂y∂v​)dxdy， ( f , v ) = ∬ Ω f v d x d y (f,v) = \iint_{\Omega}fvdxdy (f,v)=∬Ω​fvdxdy。

过程：将区域 Ω \Omega Ω剖分为 N E NE NE个小三角形，即节点个数为 N P NP NP，则有

a ( u , v ) = ∬ Ω ( ∂ u ∂ x ∂ v ∂ x + ∂ u ∂ y ∂ v ∂ y ) d x d y             = ∑ l = 1 N E ∬ e l ( ∂ u ∂ x ∂ v ∂ x + ∂ u ∂ y ∂ v ∂ y ) d x d y a(u,v) = \iint_{\Omega}(\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial y})dxdy\\\ \ \ \ \ \ \ \ \ \ \ =\sum_{l=1}^{NE}\iint_{e_l}(\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial y})dxdy a(u,v)=∬Ω​(∂x∂u​∂x∂v​+∂y∂u​∂y∂v​)dxdy           =∑l=1NE​∬el​​(∂x∂u​∂x∂v​+∂y∂u​∂y∂v​)dxdy

( f , v ) = ∬ Ω f v d x d y (f,v) = \iint_{\Omega}fvdxdy (f,v)=∬Ω​fvdxdy

单元刚度矩阵：

在单元 e l e_l el​上，有 u ( x , y ) ≈ u i K i ( x , y ) + u j K j ( x , y ) + u m K m ( x , y ) u(x,y) \approx u_iK_i(x,y)+u_jK_j(x,y)+u_mK_m(x,y) u(x,y)≈ui​Ki​(x,y)+uj​Kj​(x,y)+um​Km​(x,y)，同样，可以使用如上基函数对 v ( x , y ) v(x,y) v(x,y)逼近，即

v ( x , y ) ≈ v i K i ( x , y ) + v j K j ( x , y ) + v m K m ( x , y ) v(x,y) \approx v_iK_i(x,y)+v_jK_j(x,y)+v_mK_m(x,y) v(x,y)≈vi​Ki​(x,y)+vj​Kj​(x,y)+vm​Km​(x,y)

其中， u i , v i u_i,v_i ui​,vi​分别是函数 u ( ⋅ , ⋅ ) u(\cdot,\cdot) u(⋅,⋅)和 v ( ⋅ , ⋅ ) v(\cdot,\cdot) v(⋅,⋅)在 x i , y i x_i,y_i xi​,yi​处的取值， u j , u m u_j,u_m uj​,um​类似。需要注意到点 ( x i , y i ) , ( x j , y j ) (x_i,y_i),(x_j,y_j) (xi​,yi​),(xj​,yj​)和 ( x m , y m ) (x_m,y_m) (xm​,ym​)是三角单元 e l e_l el​的顶点。

因此

∬ e l ( ∂ u ∂ x ∂ v ∂ x + ∂ u ∂ y ∂ v ∂ y ) d x d y = ∬ e l [ v i v j v m ] [ k i i k i j k i m k j i k j j k j m k m i k m j k m m ] [ u i u j u m ] d x d y \iint_{e_l}(\frac{\partial u}{\partial x}\frac{\partial v}{\partial x}+\frac{\partial u}{\partial y}\frac{\partial v}{\partial y})dxdy\\=\iint_{e_l}\begin{bmatrix}v_i&v_j&v_m \end{bmatrix}\begin{bmatrix}k_{ii}&k_{ij}&k_{im} \\k_{ji}&k_{jj}&k_{jm} \\k_{mi}&k_{mj}&k_{mm} \end{bmatrix}\begin{bmatrix}u_i\\u_j\\u_m \end{bmatrix}dxdy ∬el​​(∂x∂u​∂x∂v​+∂y∂u​∂y∂v​)dxdy=∬el​​[vi​​vj​​vm​​] ​kii​kji​kmi​​kij​kjj​kmj​​kim​kjm​kmm​​ ​ ​ui​uj​um​​ ​dxdy

其中 k s t = ( ∂ K s ∂ x ∂ K s ∂ x + ∂ K t ∂ y ∂ K t ∂ y ) ,        s , t = i , j , m k_{st}=(\frac{\partial K_s}{\partial x}\frac{\partial K_s}{\partial x}+\frac{\partial K_t}{\partial y}\frac{\partial K_t}{\partial y}),\ \ \ \ \ \ s,t=i,j,m kst​=(∂x∂Ks​​∂x∂Ks​​+∂y∂Kt​​∂y∂Kt​​),      s,t=i,j,m

记矩阵 [ K e ] = ∬ e [ k i i k i j k i m k j i k j j k j m k m i k m j k m m ] d x d y [K_e]=\iint_{e}\begin{bmatrix}k_{ii}&k_{ij}&k_{im} \\k_{ji}&k_{jj}&k_{jm} \\k_{mi}&k_{mj}&k_{mm} \end{bmatrix}dxdy [Ke​]=∬e​ ​kii​kji​kmi​​kij​kjj​kmj​​kim​kjm​kmm​​ ​dxdy为单元刚度矩阵。

单位荷载向量：

同理可得 ( f , v ) = ∬ e l f v d x d y             = ∬ e l [ v i v j v m ] [ K i K j K m ] f ( x , y ) d x d y (f,v) = \iint_{e_l}fvdxdy\\\ \ \ \ \ \ \ \ \ \ \ =\iint_{e_l}\begin{bmatrix}v_i&v_j&v_m \end{bmatrix}\begin{bmatrix}K_{i}\\K_{j}\\K_{m}\end{bmatrix}f(x,y)dxdy (f,v)=∬el​​fvdxdy           =∬el​​[vi​​vj​​vm​​] ​Ki​Kj​Km​​ ​f(x,y)dxdy，

记向量 { F e } = ∬ e [ K i K j K m ] f ( x , y ) d x d y \{F_e\}=\iint_{e}\begin{bmatrix}K_{i}\\K_{j}\\K_{m}\end{bmatrix}f(x,y)dxdy {Fe​}=∬e​ ​Ki​Kj​Km​​ ​f(x,y)dxdy

对于不同的单元 e l , l = 1 , 2 , 3... N E e_l,l=1,2,3...NE el​,l=1,2,3...NE来说，其对应的 i , j , m i,j,m i,j,m是不一样的，因此，为使在不同单元上的单元刚度矩阵和单元荷载能够求和，将单位刚度矩阵和单位荷载向量分别扩充为 N P × N P NP\times NP NP×NP的矩阵和 N P × 1 NP\times 1 NP×1的向量，如下：

[ K e ] ˉ = ∬ e [ ⋮ ⋮ ⋮ ⋯ k i i ⋯ k i j ⋯ k i m ⋯ ⋮ ⋮ ⋮ ⋯ k j i ⋯ k j j ⋯ k j m ⋯ ⋮ ⋮ ⋮ ⋯ k m i ⋯ k m j ⋯ k m m ⋯ ⋮ ⋮ ⋮ ] d x d y \bar{[K_e]}=\iint_{e}\begin{bmatrix} & \vdots& &\vdots & &\vdots & \\ \cdots& k_{ii}&\cdots&k_{ij}&\cdots&k_{im}&\cdots \\ & \vdots& &\vdots & &\vdots & \\ \cdots&k_{ji}&\cdots&k_{jj}&\cdots&k_{jm}&\cdots \\ & \vdots& &\vdots & &\vdots & \\ \cdots&k_{mi}&\cdots&k_{mj}&\cdots&k_{mm}&\cdots \\ & \vdots& &\vdots & &\vdots &\end{bmatrix}dxdy [Ke​]ˉ​=∬e​ ​⋯⋯⋯​⋮kii​⋮kji​⋮kmi​⋮​⋯⋯⋯​⋮kij​⋮kjj​⋮kmj​⋮​⋯⋯⋯​⋮kim​⋮kjm​⋮kmm​⋮​⋯⋯⋯​ ​dxdy

{ F e } ˉ = ∬ e [ ⋮ K i ⋮ K j ⋮ K m ⋮ ] f ( x , y ) d x d y \bar{\{F_e\}}=\iint_{e}\begin{bmatrix}\vdots\\K_{i}\\\vdots\\K_{j}\\\vdots\\K_{m}\\\vdots\end{bmatrix}f(x,y)dxdy {Fe​}ˉ​=∬e​ ​⋮Ki​⋮Kj​⋮Km​⋮​ ​f(x,y)dxdy

上面两个矩阵中，其中虚点表达0元素，也就是在矩阵 [ K e ] ˉ \bar{[K_e]} [Ke​]ˉ​中至多有九个非零元，非零元的位置如上；在 { F e } ˉ \bar{\{F_e\}} {Fe​}ˉ​中至多有三个非零元，非零元的位置如上。记 U = [ u 1 u 2 ⋮ u N P ] U = \begin{bmatrix} u_1\\u_2\\\vdots\\u_{NP}\end{bmatrix} U= ​u1​u2​⋮uNP​​ ​, V = [ v 1 v 2 ⋮ v N P ] V = \begin{bmatrix} v_1\\v_2\\\vdots\\v_{NP}\end{bmatrix} V= ​v1​v2​⋮vNP​​ ​

记 K ˉ = ∑ l N E [ K e l ] ˉ \bar K=\sum_l^{NE}\bar{[K_{e_l}]} Kˉ=∑lNE​[Kel​​]ˉ​， F ˉ = ∑ l N E { F e l } ˉ \bar F=\sum_l^{NE}\bar{\{F_{e_l}\}} Fˉ=∑lNE​{Fel​​}ˉ​

则 a ( u , v ) = ( f , v ) a(u,v)=(f,v) a(u,v)=(f,v)近似于 V T K ˉ U = V T F ˉ V^T\bar{K}U=V^T\bar{F} VTKˉU=VTFˉ，即 V T ( K ˉ U − K ˉ ) = 0 V^T(\bar{K}U-\bar{K})=0 VT(KˉU−Kˉ)=0,由于 v ( x , y ) v(x,y) v(x,y)的任意性，因此对于任意的 V V V， V T ( K ˉ U − F ˉ ) = 0 V^T(\bar{K}U-\bar{F})=0 VT(KˉU−Fˉ)=0总是成立的，故 K ˉ U − F ˉ = 0 \bar{K}U-\bar{F}=0 KˉU−Fˉ=0。

理论可能写的不够详细，感觉有点乱，见谅。

在有限元-区域划分博客中，已经着重介绍了矩形区域上三角网格的剖分，包括保存节点位置的矩阵 P o s Pos Pos和保存三角形顶点位置的连接矩阵 T T T，将其封装为函数，如下

在线性插值函数的基函数构造博客中，已经着重介绍了如何构造过不在同一直线上的三个点的平面，以及整理后得到的线性插值的基函数，这里将其封装，通过输入三个不在同一直线上的节点，返回三个插值基函数，具体函数如下：

在函数在任意三角区域二重积分的计算博客中，已经着重介绍了任意三角区域二重积分的计算，在下面这个函数中，你只需要输入三角形三个顶点对应坐标以及函数 f 1 ( x , y ) f1(x,y) f1(x,y)，就可以得到函数 f 1 ( x , y ) f1(x,y) f1(x,y)在这三个点确定的三角区域上的二重积分值。

根据理论知识可以构造荷载向量，代码如下：

刚度矩阵的构造同理：

最后处理边界点，首先找到边界点的对应的位置，将刚度矩阵中对应行替换为向量 ( 0 , 0 , . . . 0 , 1 , 0 , . . . 0 ) (0,0,...0,1,0,...0) (0,0,...0,1,0,...0)，向量中1的位置为边界点对应在向量U中的位置，由于边界值是已知的，进一步替换荷载向量对应位置为所取边界值。

这个实例是博客有限元方法入门：有限元方法简单的二维算例（三角形剖分）中使用的例子， − Δ u = f , u ∣ ∂ Ω = 0 -\Delta u=f,u|_{\partial \Omega}=0 −Δu=f,u∣∂Ω​=0，其中 Ω = ( 0 , 1 ) 2 \Omega = ( 0 , 1 ) ^2 Ω=(0,1)2且 f = 2 π 2 sin ⁡ ( π x ) sin ⁡ ( π y ) f=2\pi^2\sin(\pi x)\sin(\pi y) f=2π2sin(πx)sin(πy)，真解为 u ( x , y ) = sin ⁡ ( π x ) sin ⁡ ( π y ) u(x,y)=\sin(\pi x)\sin(\pi y) u(x,y)=sin(πx)sin(πy)。

结果如下：

实例代码

这个例子在《数值方法(MATLAB版)(第四版)》（John H.Mathews与Kurtis D.Fink所著）第402页例10.7，如下：

所得数值解

实例代码如下：

这段代码和上面的代码几乎一样，只是右端函数不同，以及边界处的赋值代码不一样：

注意：使用到的函数实例只是验证我们的代码可以求解问题。