有限元之有限元法的实现

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

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

三、边界条件处理

四、算例实现

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+1for(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+1lnd[0][1]=1;lnd[0][2]=m+1;lnd[1][0]=m+2;  //1号单元的三个节点局部编号是0,1,2,整体编号是m+2,m+1,1lnd[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及右端项rhsA=(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=rhsu=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的直径。