有限元之有限元法的實現

一、單元剛度矩陣及單元荷載

二、總剛度矩陣及總荷載的合成

三、邊界條件處理

四、算例實現

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的直徑。