有限元之有限元法的实现

一、单元刚度矩阵及单元荷载

二、总刚度矩阵及总荷载的合成

三、边界条件处理

四、算例实现

4.1 C++代码

4.2 计算结果

五、结论

前三节我们介绍了有限元的基本概念、变分理论及有限元空间的构造，本节我们探讨如何实现有限元法。我们继续以二维椭圆型方程的初边值问题（公式1）为研究对象，探讨以三角形一次有限元来数值求解该问题，也就是 $V_{h}\subset H^{1}_{0}(\Omega)$ 为分片连续的一次多项式函数空间，求解离散的变分问题式（公式2），

$\left\{\begin{matrix} -\Delta u=-(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})=f(x,y),(x,y)\in \overset{^{\circ}}{\Omega}\\ u(x,y)=0,\;\;\; (x,y)\in \partial \Omega = \Gamma \end{matrix}\right.\;\;\;(1)$

求 $u_{h}(x,y)\in V_{h}$ ，使得 $a(u_{h},v_{h})=(f,v_{h})\;\;\; \forall v_{h}(x,y)\in V_{h}\;\;\;(2)$

即求 $u_{h}(x,y)\in V_{h}$ ，使得

$\iint_{\Omega}\bigtriangledown u_{h}\cdot \bigtriangledown v_{h}dxdy=\iint_{\Omega}fv_{h}dxdy,\;\;\;\;\forall v_{h}(x,y)\in V_{h}$

也就是

$\underset{e\in \tau_{h}}{\sum }\iint_{e} \bigtriangledown u_{h} \cdot \bigtriangledown v_{h} dxdy=\underset{e\in \tau_{h}}{\sum }\iint_{e} fv_{h}dxdy,\;\;\;\;\forall v_{h}(x,y)\in V_{h}\;\;(3)$

一、单元刚度矩阵及单元荷载

考虑 $P_{i},P_{j},P_{k}$ 为顶点的单元e。因为 $u_{h}|_{e}=u_{i}\lambda _{0}+u_{j}\lambda _{1}+u_{k}\lambda _{2}$ 以及 $\forall v_{h}(x,y)\in V_{h}$ 的任意性，其中 $u_{i},u_{j},u_{k}$ 待定，则有

$\underset{e\in \tau_{h}}{\sum }\iint_{e}(u_{i}\bigtriangledown \lambda_{0}+u_{j}\bigtriangledown \lambda_{1}+u_{k}\bigtriangledown \lambda_{2})\cdot \bigtriangledown \lambda_{1} dxdy=\underset{e\in \tau_{h}}{\sum }\iint_{e}f\lambda_{l}dxdy,\;\;l=0,1,2$

这里的 $\lambda_{0},\lambda_{1},\lambda_{2}$ 都是相应于具体单元e的。这样很容易得到单元e上相应于 $u_{i},u_{j},u_{k}$ 的单元刚度矩阵及单元荷载分别为

$\begin{pmatrix} \iint_{e}\bigtriangledown \lambda_{0}\cdot \bigtriangledown \lambda_{0}dxdy & \iint_{e}\bigtriangledown \lambda_{1}\cdot \bigtriangledown \lambda_{0}dxdy & \iint_{e}\bigtriangledown \lambda_{2}\cdot \bigtriangledown \lambda_{0}dxdy\\ \iint_{e}\bigtriangledown \lambda_{0}\cdot \bigtriangledown \lambda_{1}dxdy & \iint_{e}\bigtriangledown \lambda_{1}\cdot \bigtriangledown \lambda_{1}dxdy &\iint_{e}\bigtriangledown \lambda_{2}\cdot \bigtriangledown \lambda_{1}dxdy \\ \iint_{e}\bigtriangledown \lambda_{0}\cdot \bigtriangledown \lambda_{2}dxdy & \iint_{e}\bigtriangledown \lambda_{1}\cdot \bigtriangledown \lambda_{2}dxdy & \iint_{e}\bigtriangledown \lambda_{2}\cdot \bigtriangledown \lambda_{2}dxdy \end{pmatrix},\begin{pmatrix} \iint_{e}f\lambda_{0}dxdy\\ \iint_{e}f\lambda_{1}dxdy\\ \iint_{e}f\lambda_{2}dxdy \end{pmatrix}$

在实际单元刚度矩阵的计算中，由于 $\lambda_{0},\lambda_{1},\lambda_{2}$ 都是一次多项式，所以 $\bigtriangledown \lambda_{l}(l=0,1,2)$ 均为常向量，可以得到

$\bigtriangledown \lambda_{0}=\frac{1}{2S_{e}}\begin{pmatrix} y_{j}-y_{k}\\ x_{k}-x_{j} \end{pmatrix},\bigtriangledown \lambda_{1}=\frac{1}{2S_{e}}\begin{pmatrix} y_{k}-y_{i}\\ x_{i}-x_{k} \end{pmatrix},\bigtriangledown \lambda_{2}=\frac{1}{2S_{e}}\begin{pmatrix} y_{i}-y_{j}\\ x_{j}-x_{i} \end{pmatrix}$

从而单位刚度矩阵可以写作

$\frac{1}{4S_{e}}\begin{pmatrix} (y_{j}-y_{k})^{2} & (y_{k}-y_{i})(y_{j}-y_{k}) & (y_{i}-y_{j})(y_{j}-y_{k})\\ (y_{j}-y_{k})(y_{k}-y_{i}) & (y_{k}-y_{i})^{2} & (y_{i}-y_{j})(y_{k}-y_{i})\\ (y_{j}-y_{k})(y_{i}-y_{j}) & (y_{k}-y_{i})(y_{i}-y_{j}) & (y_{i}-y_{j})^{2} \end{pmatrix}+$

$\frac{1}{4S_{e}}\begin{pmatrix} (x_{k}-x_{j})^{2} & (x_{i}-x_{k})(x_{k}-x_{j}) & (x_{j}-x_{i})(x_{k}-x_{j})\\ (x_{k}-x_{j})(x_{i}-x_{k}) & (x_{i}-x_{k})^{2} & (x_{j}-x_{i})(x_{j}-x_{k})\\ (x_{k}-x_{j})(x_{j}-x_{i}) & (x_{i}-x_{k})(x_{j}-x_{i}) & (x_{j}-x_{i})^{2} \end{pmatrix}$

至于单元荷载的计算，通常情况下需要借助数值积分。常用的二维Hammer数值积分公式（在标准单元 $\widehat{e}$ 上）为

$\iint_{\widehat{e}}g(\lambda_{0},\lambda_{1},\lambda_{2})d\lambda_{0}d\lambda_{1}=\int_{0}^{1}d\lambda_{0}\int_{0}^{1-\lambda_{0}}g(\lambda_{0},\lambda_{1},1-\lambda_{0}-\lambda_{1})d\lambda_{1}=\frac{1}{2}g(\frac{1}{3},\frac{1}{3},\frac{1}{3})+O(h^{2})\;\;\;(4)$

或

$\iint_{\widehat{e}}g(\lambda_{0},\lambda_{1},\lambda_{2})d\lambda_{0}d\lambda_{1}=\frac{1}{6}(g(0,\frac{1}{2},\frac{1}{2})+g(\frac{1}{2},0,\frac{1}{2})+g(\frac{1}{2},\frac{1}{2},0))+O(h^{3})\;\;(4)$

因此在计算单元荷载时，先要将一般单元e上的积分通过坐标变换式（公式5）变到标准单元 $\widehat{e}$ 上，

$\left\{\begin{matrix} x=x_{i}\lambda_{0}+x_{j}\lambda_{1}+x_{k}\lambda_{2}\\ y=y_{i}\lambda_{0}+y_{j}\lambda_{1}+y_{k}\lambda_{2} \end{matrix}\right.\;\;\;\; or\;\;\;\left\{\begin{matrix} x=(x_{i}-x_{k})\lambda_{0}+(x_{j}-x_{k})\lambda_{1}+x_{k}\\ y=(y_{i}-y_{k})\lambda_{0}+(y_{j}-y_{k})\lambda_{1}+y_{k} \end{matrix}\right.\;\;\;(5)$

从而有

$\iint_{e}p(x,y)dxdy=2S_{e}\iint_{\widehat{e}}p(x_{i}\lambda_{0}+x_{j}\lambda_{1}+x_{k}\lambda_{2},y_{i}\lambda_{0}+y_{j}\lambda_{1}+y_{k}\lambda_{2})d\lambda_{0}d\lambda_{1}$

再利用Hammer积分式（公式4）并忽略高阶小项可得

$\iint_{e}p(x,y)dxdy\approx S_{e}\cdot p(G),$ 误差为 $O(h^{2})$ ，其中G为单元e的重心 $(6)$

或者

$\iint_{e}p(x,y)dxdy\approx \frac{S_{e}}{3}p((P_{jk})+p(P_{ki})+p(P_{ij}))$ ，误差为 $O(h^{3})\;\;\;(6)$

其中， $P_{jk},P_{ki},P_{ij}$ 为单元e的3条边的中点。更高精度的Hammer积分需要借助更多的积分点。这样，由公式（6）单元荷载可近似为

$\frac{S_{e}f(G)}{3}\begin{pmatrix} 1\\ 1\\ 1 \end{pmatrix} \;\;or\;\;\frac{S_{e}}{6}\begin{pmatrix} f(P_{ki})+f(P_{ij})\\ f(P_{ij})+f(P_{jk})\\ f(P_{jk})+f(P_{ki}) \end{pmatrix}\;\;\;(7)$

如果 $V_{h}$ 为分片连续二次多项式函数空间，则单元刚度矩阵将会是6x6阶的矩阵，计算会更复杂。

二、总刚度矩阵及总荷载的合成

为简单起见，以图1的剖分为例进行单元刚度矩阵的合成。可知，总刚度矩阵A是一个 $30((m+1)(n+1))$ 阶的方阵。在进行总刚度矩阵A的合成之前，需要初始化A，让其所有元素均为0。接下来我们考察第⑩号单元的单元刚度矩阵在总刚度矩阵A中的位置。由于第⑩号单元的3个顶点的整体编号为6,7,12，局部编号为0,1,2，可见其相应于 $u_{6},u_{7},u_{12}$ 单元刚度矩阵为

$\frac{1}{4S_{e}}\begin{pmatrix} (y_{7}-y_{12})^{2} & (y_{12}-y_{6})(y_{7}-y_{12}) & (y_{6}-y_{7})(y_{7}-y_{12})\\ (y_{7}-y_{12})(y_{12}-y_{6}) & (y_{12}-y_{6})^{2} & (y_{6}-y_{7})(y_{12}-y_{6})\\ (y_{7}-y_{12})(y_{6}-y_{7}) & (y_{12}-y_{6})(y_{6}-y_{7}) & (y_{6}-y_{7})^{2} \end{pmatrix}+$

$\frac{1}{4S_{e}}\begin{pmatrix} (x_{12}-x_{7})^{2} & (x_{6}-x_{12})(x_{12}-x_{7}) & (x_{7}-x_{6})(x_{12}-x_{7})\\ (x_{12}-x_{7})(x_{6}-x_{12}) & (x_{6}-x_{12})^{2} & (x_{7}-x_{6})(x_{7}-x_{12})\\ (x_{12}-x_{7})(x_{7}-x_{6}) & (x_{6}-x_{12})(x_{7}-x_{6}) & (x_{7}-x_{6})^{2} \end{pmatrix}$

其中

$S_{e}=\frac{(x_{b}-x_{a})(y_{d}-y_{c})}{2mn}$

在实际编程中，上述三阶单元刚度矩阵的元素可以存储在一个二维数组ea中。如，令

$ea[1][0]=\frac{(y_{7}-y_{12})(y_{12}-y_{6})+(x_{12}-x_{7})(x_{6}-x_{12})}{4S_{e}}$

$ea[2][2]=\frac{(y_{6}-y_{7})^{2}+(x_{7}-x_{6})^{2}}{4S_{e}}$

然后再将数组ea合成至总刚度矩阵A中，其中ea[1][0]在A中的位置是A[7][6]，ea[2][2]在A中的位置是A[12][12]，ea[1][2]在A中的位置是A[7][12]等等。这里用到了局部编号与整体编号之间的对应关系，这样将ea中的9个元素都存到A中相应的位置。同样处理第⑪号单元，但由于第⑪号单元的顶点中也有整体编号为7和12分别对应局部编号2和1的节点，因此这个单元上的单元刚度矩阵中ea[2][1]在A中的位置也是A[7][12]，它需要与已有的A[7][12]的值相加，从而实现总刚度矩阵的合成。最后，对所有单元都进行上述过程，则可最终得到总刚度矩阵A。

类似的，进行总荷载的合成。可知，总荷载矩阵rhs（表示右端项right-hand-side）是一个 $30((m+1)(n+1))$ 维的列向量。在进行总荷载矩阵rhs的合成之前，也要初始化，让其全部元素为0。同时，考察第⑩号单元的单元荷载在总单元荷载中的位置。可以取单元荷载为公式（7）中精度较高的那个，即相应于 $u_{6},u_{7},u_{12}$ 的单元荷载为

$\frac{S_{e}}{6}\begin{pmatrix} f(P_{12,6})+f(P_{6,7})\\ f(P_{6,7})+f(P_{7,12})\\ f(P_{7,12})+f(P_{12,6}) \end{pmatrix}=\frac{S_{e}}{6}\begin{pmatrix} f(\frac{x_{12}+x_{6}}{2},\frac{y_{12}+y_{6}}{2})+f(\frac{x_{6}+x_{7}}{2},\frac{y_{6}+y_{7}}{2}) \\ f(\frac{x_{6}+x_{7}}{2},\frac{y_{6}+y_{7}}{2})+f(\frac{x_{7}+x_{12}}{2},\frac{y_{7}+y_{12}}{2}) \\ f(\frac{x_{7}+x_{12}}{2},\frac{y_{7}+y_{12}}{2})+f(\frac{x_{12}+x_{6}}{2},\frac{y_{12}+y_{6}}{2}) \end{pmatrix}$

在实际编程中，上述列向量可以存储在一个一维数组g中，这样，g[0]在总荷载rhs中的位置是rhs[6]，g[1]的位置是rhs[7]，g[2]的位置是rhs[12]。再处理第⑪号单元，由于第⑪号单元的顶点中也有整体编号为7和12分别对应局部编号为2和1的节点，因此这个单元上的单元荷载中g[2]在A中的位置也是rhs[7]，它需要与已有的rhs[7]的值相加，g[1]在A中的位置也是rhs[12]，它需要与已有的rhs[12]的值相加，从而实现总荷载的合成。最后，对所有单元都进行上述过程，则可最终得到总荷载矩阵rhs。

综上，可得原离散问题的矩阵表示为

$Au=rhs,\;\;u=(u_{0},u_{1},\cdots,u_{28},u_{29})^{T}\;\;(8)$

三、边界条件处理

在进行线性方程组（公式8）求解前，还需要对边界条件进行处理。由于原问题的边界条件是 $u(x,y)|_{\partial \Omega }=0$ ，且有限元空间要求 $V_{h}\subset H^{1}_{0}(\Omega)$ ，这样就限制了 $u_{h}(P)=0$ ，其中P为边界节点，也就是要求 $u_{l}=0$ （ $l\in \Lambda =$ {所有边界节点的整体编号}）。这样，只要修改系数矩阵A，令其第 $l$ 行和第 $l$ 列元素均为0，但 $A[l][l]=1,l\in\Lambda$ ，得到新的系数矩阵记作 $\tilde{A}$ 。再令右端项 $rhs[l]=0,l\in\Lambda$ ，得到新的右端项记作b。完成以上修改后再去求解 $\tilde{A}u=b$ ，即可得到 $u_{0},u_{1},\cdots,u_{28},u_{29}$ ，从而得到数值解 $u_{h}$ 。

事实上，连续双线性型 $a(\cdot,\cdot)$ 的对称性就确定了有限元离散后的总刚度矩阵A是对称的， $a(\cdot,\cdot)$ 的V-椭圆性就保证了这个刚度矩阵还是正定的，即使做了边界条件处理后仍然是对称正定的，从而离散后的线性系统是唯一可解的，可以利用高斯消元法求解。

四、算例实现

利用三角形一次有限元求解椭圆型方程边值问题：

$\left\{\begin{matrix} -(\frac{\partial^{2}u(x,y)}{\partial x^{2}}+\frac{\partial^{2}u(x,y)}{\partial y^{2}})=2(x+y-x^{2}-y^{2}),(x,y)\in \overset{\circ }{\Omega}\\ u|_{\partial \Omega}=0 \end{matrix}\right.$

其中， $\Omega =[0,1]\times[0,1]$ 。已知此问题的精确解为 $u(x,y)=xy(1-x)(1-y)$ 。分别取第一种剖分m=n=16和第二种剖分m=n=32，输出在节点 $(x,y)=(0.125i,0.5),1\leqslant i\leqslant 7$ 处的数值解，并给出误差。

4.1 C++代码

#include <cmath>
#include <stdlib.h>
#include <stdio.h>

int m,n;
int node,elem,limnode;
int **lnd;
double area,*xco,*yco,*u;

int main(int argc, char *argv[])
{
        int i,j,k,t,e,row,col;
        int *limnd;
        double x,y,q,sum,summ;
        double a,b,c,d,dx,dy;
        double **A,*rhs,*w,*graduh;
        double alpha[3],beta[3],ea[3][3],g[3];
        double v(double x, double y);
        double vx(double x, double y);
        double vy(double x, double y);
        double f(double x, double y);
        double *fun_lambda(int e, double x, double y);
        double *d_lambda(int e);
        double uh(int e, double x, double y);
        double *d_uh(int e);
        double *GaussElimination(double **a, double *b, int N);

        a=0.0;  //x方向的求解区域[a,b]
        b=1.0;
        c=0.0;  //y方向的求解区域[c,d]
        d=1.0;
        m=16;  //x方向的剖分数
        n=16;  //y方向的剖分数
        printf("m=%d,n=%d.\n",m,n);

        dx=(b-a)/m;  //x方向的步长
        dy=(d-c)/n;  //y方向的步长
        node=(m+1)*(n+1);  //总节点数
        elem=2*m*n;  //总单元数（三角形剖分）
        limnode=2*(m+n);  //边界节点数（受限节点数）
        area=(b-a)*(d-c)/elem;  //单元面积

        limnd=(int *)malloc(sizeof(int)*(limnode));
        for(i=0;i<=m;i++)   //边界节点的编号
        {
                limnd[i]=i;  //下底边节点编号
                limnd[limnode-i-1]=node-1-i;  //上顶边节点编号
        }   //此时i=m+1
        for(j=1;j<n;j++)
        {
                limnd[i]=j*(m+1);  //左侧边上的节点编号
                limnd[i+1]=limnd[i]+m;  //右侧边上的节点编号
                i=i+2;
        }

        //各单元局部节点编号与整体编号之间的关系
        lnd=(int **)malloc(sizeof(int*)*elem);
        for(i=0;i<elem;i++)
                lnd[i]=(int*)malloc(sizeof(int)*3);
        lnd[0][0]=0;  //0号单元的三个节点局部编号是0,1,2,整体编号是0,1,m+1
        lnd[0][1]=1;
        lnd[0][2]=m+1;
        lnd[1][0]=m+2;  //1号单元的三个节点局部编号是0,1,2,整体编号是m+2,m+1,1
        lnd[1][1]=m+1;
        lnd[1][2]=1;

        for(e=2;e<2*m;e++)  //2~2m-1号单元的节点编号情况
        {
                for(i=0;i<3;i++)
                        lnd[e][i]=lnd[e-2][i]+1;
        }
        for(e=2*m;e<elem;e++)  //2m-1~elem-1号单元的节点编号情况
        {
                for(i=0;i<3;i++)
                        lnd[e][i]=lnd[e-2*m][i]+m+1;
        }

        //各节点的坐标
        xco=(double*)malloc(sizeof(double)*node);
        yco=(double*)malloc(sizeof(double)*node);
        for(i=0;i<=m;i++)
        {
                xco[i]=a+i*dx;
                yco[i]=c;
                for(j=1;j<=n;j++)
                {
                        t=j*(m+1)+i;
                        xco[t]=xco[i];
                        yco[t]=c+j*dy;
                }
        }

        //初始化刚度矩阵A及右端项rhs
        A=(double**)malloc(sizeof(double*)*node);
        for(i=0;i<node;i++)
                A[i]=(double*)malloc(sizeof(double)*node);
        rhs=(double*)malloc(sizeof(double)*node);
        u=(double*)malloc(sizeof(double)*node);
        for(i=0;i<node;i++)
        {
                rhs[i]=0.0;
                for(j=0;j<node;j++)
                        A[i][j]=0.0;
        }
        
        for(e=0;e<elem;e++)
        {
                i=lnd[e][0];
                j=lnd[e][1];
                k=lnd[e][2];

                //单元荷载
                g[0]=area*(f((xco[i]+xco[k])/2,(yco[i]+yco[k])/2)/2+f((xco[i]+xco[j])/2,(yco[i]+yco[j])/2)/2)/3;
                g[1]=area*(f((xco[j]+xco[k])/2,(yco[j]+yco[k])/2)/2+f((xco[i]+xco[j])/2,(yco[i]+yco[j])/2)/2)/3;
                g[2]=area*(f((xco[i]+xco[k])/2,(yco[i]+yco[k])/2)/2+f((xco[k]+xco[j])/2,(yco[k]+yco[j])/2)/2)/3;

                w=d_lambda(e);
                //计算单元刚度矩阵
                for(i=0;i<3;i++)
                        for(j=0;j<3;j++)
                                ea[i][j]=(w[i]*w[j]+w[i+3]*w[j+3])*area;

                //合成总刚度矩阵
                for(i=0;i<3;i++)
                {
                        for(j=0;j<3;j++)
                        {
                                row=lnd[e][i];
                                col=lnd[e][j];
                                A[row][col]=A[row][col]+ea[i][j];
                        }
                }

                //合成总荷载
                for(i=0;i<3;i++)
                {
                        k=lnd[e][i];  //确定合成总荷载时所在的行
                        rhs[k]=g[i]+rhs[k];
                }
        }

        //处理边界条件
        for(i=0;i<limnode;i++)
        {
                k=limnd[i];
                for(t=0;t<node;t++)
                {
                        A[t][k]=0;
                        A[k][t]=0;
                }
                A[k][k]=1;
                rhs[k]=0;
        }

        //高斯消元法解方程组Au=rhs
        u=GaussElimination(A,rhs,node);

        //输出
        for(i=1;i<m;i++)
        {
                t=(node-1-m)/2+i;   //y=0.5时对应的节点编号
                if(t%(m/8)==0)
                        printf("u[%d]=%f, x=%.3f, err=%.4e.\n",t,u[t],xco[t],fabs(u[t]-v(xco[t],0.5)));
        }

        //误差估计
        sum=0;
        summ=0;
        for(e=0;e<elem;e++)
        {
                i=lnd[e][0];
                j=lnd[e][1];
                k=lnd[e][2];
                alpha[0]=(xco[j]+xco[k])/2;
                alpha[1]=(xco[k]+xco[i])/2;
                alpha[2]=(xco[i]+xco[j])/2;
                beta[0]=(yco[j]+yco[k])/2;
                beta[1]=(yco[k]+yco[i])/2;
                beta[2]=(yco[i]+yco[j])/2;
                graduh=d_uh(e);

                for(i=0;i<3;i++)
                {
                        q=v(alpha[i],beta[i])-uh(e,alpha[i],beta[i]); //计算L2范数
                        sum=sum+q*q;
                        x=vx(alpha[i],beta[i])-graduh[0];
                        y=vy(alpha[i],beta[i])-graduh[1];  //计算H1范数
                        summ=summ+x*x+y*y;
                }
        }

        sum=sum*area/3;
        summ=summ*area/3;
        printf("0 norm error is:%.8f.\n",sqrt(sum));
        printf("1 norm error is:%.8f.\n",sqrt(sum+summ));

        for(e=0;e<elem;e++)
                free(lnd[e]);
        for(i=0;i<node;i++)
                free(A[i]);
        free(lnd);free(A);free(rhs);free(xco);free(yco);
        free(limnd);free(u);free(w);free(graduh);

        return 0;
}


//精确解
double v(double x, double y)
{
        double z;
        z=x*y*(1-x)*(1-y);
        return z;
}
//精确解关于x的偏导数
double vx(double x, double y)
{
        double z;
        z=(1-2*x)*(y-y*y);
        return z;
}
//精确解关于y的偏导数
double vy(double x, double y)
{
        double z;
        z=(1-2*y)*(x-x*x);
        return z;
}
//方程的右端项
double f(double x, double y)
{
        double z;
        z=2*(x+y-x*x-y*y);
        return z;
}
//单元e上的基函数λ0,λ1,λ2
double *fun_lambda(int e, double x, double y)
{
        int i,j,k;
        double *lambda;

        lambda=(double*)malloc(sizeof(double)*3);
        i=lnd[e][0];
        j=lnd[e][1];
        k=lnd[e][2];
        lambda[0]=((xco[j]-x)*(yco[k]-y)-(yco[j]-y)*(xco[k]-x))/(2*area);
        lambda[1]=((xco[k]-x)*(yco[i]-y)-(yco[k]-y)*(xco[i]-x))/(2*area);
        lambda[2]=((xco[i]-x)*(yco[j]-y)-(yco[i]-y)*(xco[j]-x))/(2*area);

        return lambda;
}
//单元e上的基函数lambda的关于x,y的偏导数
double *d_lambda(int e)
{
        int i,j,k;
        double *w;

        w=(double*)malloc(sizeof(double)*6);
        i=lnd[e][0];
        j=lnd[e][1];
        k=lnd[e][2];
        w[0]=(yco[j]-yco[k])/(2*area);
        w[1]=(yco[k]-yco[i])/(2*area);
        w[2]=(yco[i]-yco[j])/(2*area);
        w[3]=(xco[k]-xco[j])/(2*area);
        w[4]=(xco[i]-xco[k])/(2*area);
        w[5]=(xco[j]-xco[i])/(2*area);

        return w;
}
//单元e上的数值解uh
double uh(int e, double x, double y)
{
        int i,k;
        double z, *lambda;

        z=0;
        lambda=fun_lambda(e,x,y);
        for(i=0;i<3;i++)
        {
                k=lnd[e][i];
                z=z+u[k]*lambda[i];
        }

        return z;
}
//单元e上uh关于x,y的偏导数
double *d_uh(int e)
{
        int i,j,k;
        double *z,*w;

        i=lnd[e][0];
        j=lnd[e][1];
        k=lnd[e][2];
        w=d_lambda(e);
        z=(double*)malloc(sizeof(double)*2);
        z[0]=z[1]=0.0;

        for(i=0;i<3;i++)
        {
                k=lnd[e][i];
                z[0]=z[0]+u[k]*w[i];
                z[1]=z[1]+u[k]*w[i+3];
        }

        return z;
}
//Gauss消去法子程序解N阶线性方程组Ax=b
double *GaussElimination(double **a, double *b, int N)
{
        int i,j,k;
        double *x,sum;

        x=(double*)malloc(sizeof(double)*N);
        for(k=0;k<N;k++)
        {
                if(fabs(a[k][k])<1e-8)
                {
                        printf("Error!\n");
                        exit;
                }
                b[k]=b[k]/a[k][k];
                for(j=k+1;j<N;j++)
                        a[k][j]=a[k][j]/a[k][k];
                for(i=k+1;i<N;i++)
                {
                        b[i]=b[i]-a[i][k]*b[k];
                        for(j=k+1;j<N;j++)
                                a[i][j]=a[i][j]-a[i][k]*a[k][j];
                }
        }

        x[N-1]=b[N-1];
        for(i=N-2;i>=0;i--)
        {
                sum=0;
                for(j=i+1;j<N;j++)
                        sum=sum+a[i][j]*x[j];
                x[i]=b[i]-sum;
        }

        return x;
}

4.2 计算结果

当m=n=16时，计算结果为：

m=16,n=16.
u[138]=0.027253, x=0.125, err=9.0703e-05.
u[140]=0.046726, x=0.250, err=1.4885e-04.
u[142]=0.058413, x=0.375, err=1.8101e-04.
u[144]=0.062309, x=0.500, err=1.9127e-04.
u[146]=0.058413, x=0.625, err=1.8101e-04.
u[148]=0.046726, x=0.750, err=1.4885e-04.
u[150]=0.027253, x=0.875, err=9.0703e-05.
0 norm error is:0.00039064.
1 norm error is:0.01518861.

当m=n=32时，计算结果为：

m=32,n=32.
u[532]=0.027321, x=0.125, err=2.2726e-05.
u[536]=0.046838, x=0.250, err=3.7299e-05.
u[540]=0.058548, x=0.375, err=4.5356e-05.
u[544]=0.062452, x=0.500, err=4.7926e-05.
u[548]=0.058548, x=0.625, err=4.5356e-05.
u[552]=0.046838, x=0.750, err=3.7299e-05.
u[556]=0.027321, x=0.875, err=2.2726e-05.
0 norm error is:0.00009800.
1 norm error is:0.00760401.

五、结论

从计算结果可知，数值结果与下面的定理一致

定理：

设 $u,u_{h}$ 分别是连续问题“求 $u(x,y)\in H^{1}_{0}(\Omega)$ ，使得 $a(u,v)=(f,v), \forall v(x,y)\in H^{1}_{0}(\Omega)$ ”和离散问题“求 $u_{h}(x,y)\in V_{h}$ ，使得 $a(u_{h},v_{h})=(f,v_{h}), \forall v_{h}(x,y)\in V_{h}$ ”的解，则对于剖分细度 $h=\underset{e\in \tau_{h}}{max}(diam(e))$ ， $\left \| u-u_{h} \right \|_{1,\Omega},\left \| u-u_{h} \right \|$ 分别是一阶、二阶收敛的。其中， $diam(e)$ 表示剖分 $\tau_{h}$ 中单元e的直径。