有限元之初探门径

蒲公英的孩子

已于 2024-05-28 15:07:32 修改

阅读量684

点赞数 11

分类专栏：有限元文章标签：算法 c++ 数学建模

于 2024-05-22 11:07:07 首次发布

本文链接：https://blog.csdn.net/L_peanut/article/details/139091636

版权

有限元专栏收录该内容

8 篇文章 3 订阅 ¥15.90 ¥99.00

订阅专栏

超级会员免费看

自然科学与工程领域中很多问题的数学模型都是以偏微分方程或积分方程的形式给出的，这些方程是很难通过解析的方式获得其精确解的。通过数值方法，我们将连续情形下的原问题（通常是无限维的）离散为一个有限维的问题，从而得到只有有限个未知量的离散系统。它可以是线性系统，也可以是非线性系统，然后再通过合适的算法来求解。

前面四个专栏我们介绍了针对不同类型偏微分方程的差分格式，差分法在求解常微分方程、抛物型方程和双曲型方程的定解问题上都比较简单，但在处理椭圆型方程的边值问题时，无论是格式设计、理论误差分析还是编程实践上都比较复杂。而解决椭圆型偏微分方程边值问题最有效的方法是有限元法（Finite Element Method，FEM）。

常微分方程_蒲公英的孩子的博客-CSDN博客https://blog.csdn.net/l_peanut/category_12628498.html 偏微分方程（抛物型）_蒲公英的孩子的博客-CSDN博客https://blog.csdn.net/l_peanut/category_12635343.html 偏微分方程（双曲型）_蒲公英的孩子的博客-CSDN博客https://blog.csdn.net/l_peanut/category_12646078.html 偏微分方程（椭圆型）_蒲公英的孩子的博客-CSDN博客https://blog.csdn.net/l_peanut/category_12654271.html

有限元法是20世纪50年代由工程师们在求解结构力学问题（如梁、板问题等）中提出的，基本做法是将复杂的结构体分成有限小块，俗称“单元”。先在每个结构性质相对简单的小单元上作分析，最后进行总的合成。在20世纪60年代中期，以中科院冯康先生为代表的中国学者和西方学者独立并行的研究了有限元法的数学理论，发现其理论基础主要是20世纪初数学理论中的变分法和Sobolev空间理论。与常规的有限差分相比，有限元法从数学物理问题的变分原理出发，从整体来描述原问题，而差分法则是局部描述原问题；有限元法对求解区域的剖分离散可以是多样的，如对二维平面区域可以使用三角形剖分，也可以使用矩形剖分，剖分后的小单元甚至可以是曲边的，而差分法则局限于处理矩形区域，而且只进行矩形剖分；有限元法最后得到的数值结果是分片的，数值解在剖分区域内的没一点都有定义，而差分法最后得到的数值解只在孤立的节点有定义。

接下来，我们在本专栏对有限元法进行简单介绍，需要深入学习的小伙伴请参考专业的书籍。

一、常微分方程两点边值问题的等价形式

考虑以下的模型问题——常微分方程两点Dirichlet边值问题：

$\left\{\begin{matrix} -u^{''}(x)=f(x),x\in(0,1),\\ u(0)=u(1)=0 \end{matrix}\right.\;\;(1)$

实际上它是一个一维的椭圆型方程边值问题。这个边值问题与一个变分公式：

求 $u\in V$ ，使得 $(u^{'},v^{'})=(f,v),\forall v\in V \;\;(2)$

以及一个泛函（从一个集合到实数集的映射）的极小问题：

求 $u\in V$ ，使得 $J(u)\leqslant J(v),\forall v \in V,$ 即 $J(u)=\underset{v \in V}{min} J(v)\;\;(3)$

相互等价。其中

V={v|v在[0,1]连续，v’在[0,1]上分片连续、有界，且v(0)=v(1)=0} $(4)$

$(v,w)= \int_{0}^{1}v(x)w(x)dx,\;J(v)=\frac{1}{2}(v^{'},v^{'})-(f,v),\forall v,w\in V\;\; (5)$

二、模型问题的有限元法

连续情形下的模型问题与变分公式等价，从而对公式（1）的分析可以转化为对公式（2）的分析。有限元法就是直接从公式（2）进行理论与数值分析的。具体情况为：

令 $V_{h}$ 为 $V$ 一个有限维子空间，通常是分片连续多项式函数空间，先对变分公式（2）进行离散化，即

求 $u_{h}\in V_{h}$ ，使得 $(u_{h}^{'},v_{h}^{'})=(f,v_{h}),\forall v_{h}\in V_{h}\;\;(6)$

即在 $V$ 的子空间 $V_{h}$ 中寻找一个形式上满足公式（6）的解 $u_{h}$ ，称为原问题式（1）的数值解。于是需要对求解区域进行剖分，为简单起见，采用等距剖分。将[0,1]分成m份，得到m+1个节点 $x_{i}=ih,i=0,1,\cdot\cdot\cdot,m$ 和m个单元 $[x_{i},x_{i+1}],i=0,1,\cdot\cdot\cdot,m-1$ 。其中，h=1/m。然后不妨取 $V_{h}$ 为分片连续一次多项式函数空间，即 $V_{h}$ 中的元素 $v(x)$ 在每个单元即子区间 $[x_{i-1},x_{i}]$ 上都是线性函数、在每个节点处都连续，且 $v(0)=v(1)=0$ 。这样的函数图像显示就是一条始于原点、终于点（1,0）的折线段。要准确地描述出这些分片连续函数，只要确定函数在各节点上的取值即可。为此引入 $V_{h}$ 中的基函数：

$\varphi_{i}(x)=\left\{\begin{matrix} \frac{x-x_{i-1}}{h},x_{i-1}\leqslant x \leqslant x_{i},\\ \frac{x_{i+1}-x}{h},x_{i}<x\leqslant x_{i+1}=\left\{\begin{matrix} 1-\frac{|x-x_{i}|}{h},x_{i-1}\leqslant x \leqslant x_{i+1},\\ 0,\;else \end{matrix}\right. 1 \leqslant i \leqslant m-1\\ 0,\;else \end{matrix}\right.$

且 $\varphi_{0}(x)=\left\{\begin{matrix} \frac{x_{1}-x}{h},x\in[x_{0},x_{1}]\\ 0,\;else \end{matrix}\right.,\varphi_{m}(x)=\left\{\begin{matrix} \frac{x-x_{m-1}}{h},x\in[x_{m-1},x_{m}]\\ 0,\;else \end{matrix}\right.$

其中， $\varphi_{i}(x)\in V_{h}$ 均定义在[0,1]上，故有

$\varphi_{i}(x_{j})=\left\{\begin{matrix} 1,i=j\\ 0,i\neq j \end{matrix}\right.$

也就是说基函数 $\varphi_{i}(x)$ 在节点 $x_{i}$ 处取值为1，在其他节点处取值为0且 $\varphi_{i}(x)(0\leqslant i \leqslant m)$ 张成一个m+1维空间 $V_{h}$ 。这样， $V_{h}$ 中的函数 $v_{h}(x)$ 就可以表示为基函数 $\varphi_{i}(x)$ 的一个线性组合，即 $v_{h}(x)=\sum_{i=0}^{m}v_{i}\varphi_{i}(x)$ 。其中，折线段 $v_{h}(x)$ 的确定依赖于系数 $v_{i}$ ，也就是 $v_{h}(x_{i})$ 。以上即为有限元法实际操作的前两步：

第一步，确定变分公式。第二步，构造有限元空间。第三步，将变分公式（6）具体写出来，得到离散系统，进行数值求解。具体的，设待求的 $u_{h}\in V_{h}$ 形式为 $u_{h}(x)=\sum_{i=0}^{m}u_{i}\varphi_{i}(x)$ ，由由于公式（6）对所有的 $v_{h}\in V_{h}$ 成立，只要使公式（6）对所有 $V_{h}$ 中的基函数 $\varphi_{j}(x),0\leqslant j\leqslant m$ 成立即可，因为 $V_{h}$ 由 $\varphi_{j}(x)$ 张成。于是有

$(\sum_{i=0}^{m}u_{i}\varphi^{'}_{i}(x),\varphi^{'}_{j}(x))=(f,\varphi_{j}(x)),0\leqslant j\leqslant m,$

即 $\sum_{i=0}^{m}u_{i}(\varphi^{'}_{i}(x),\varphi^{'}_{j}(x))=(f,\varphi_{j}(x)),0\leqslant j\leqslant m,\;\;(7)$

对 $0\leqslant j\leqslant m$ ，若记

$a_{i,j}=(\varphi^{'}_{i}(x),\varphi^{'}_{j}(x))=\int_{0}^{1}\varphi^{'}_{i}(x)\varphi^{'}_{j}(x)dx,b_{j}=(f,\varphi_{j}(x))=\int_{0}^{1}f(x)\varphi_{j}(x)dx\;\;(8)$

则公式（7）就成为

$\sum_{i=0}^{m}u_{i}a_{i,j}=\sum_{i=0}^{m}u_{i}a_{j,i}=b_{j}(0\leqslant j\leqslant m)$

也就是

$\sum_{j=0}^{m}a_{i,j}u_{j}=b_{i},0\leqslant i\leqslant m \;\;(9)$

公式（9）可以写成离散系统——线性方程组的形式： $A\mathbf{u}=\mathbf{b}$ ，其中m+1阶对称方程A中的元素为 $a_{i,j},0\leqslant i,j\leqslant m$ ，向量 $\mathbf{u}=(u_{0},u_{1},\cdot\cdot\cdot,u_{m})^{T},\mathbf{b}=(b_{0},b_{1},\cdot\cdot\cdot,b_{m})^{T}$ 。可以证明线性方程组（9）的系数矩阵A是正定的。这是因为一方面有 $\mathbf{u}^{T}A\mathbf{u}=(u^{'}_{h},u^{'}_{h})\geqslant 0$ ，即A是半正定的；另一方面当 $\mathbf{u}\neq \mathbf{0}$ 时，必有 $\mathbf{u}^{T}A\mathbf{u} >0$ ，否则， $0=\mathbf{u}^{T}A\mathbf{u} =(u^{'}_{h},u^{'}_{h})$ 就有 $u^{'}_{h}\equiv 0$ 从而 $u_{h}$ 恒为常数，再由边界条件知 $u_{h}$ 恒为0，从而 $\mathbf{u}\equiv 0$ 得到矛盾。此外，A还是三对角矩阵，因为当节点 $x_{i}$ 和 $x_{j}$ 不相邻时，即 $|i-j|>1$ 时，在以 $x_{k}$ 和 $x_{k+1}$ 为端点的单元内， $\varphi_{i}(x)$ 和 $\varphi_{j}(x)$ 至少一个恒为0，从而

$a_{i,j}=\sum_{k=0}^{m-1}\int_{x_{k}}^{x_{k+1}}\varphi^{'}_{i}(x)\varphi^{'}_{j}(x)dx=0 \;\;(|i-j|>1)$

此外，直接计算可知

$a_{i,i}=2/h,a_{i-1,i}=a_{i,i-1}=-1/h(1\leqslant i\leqslant m),a_{0,0}=a_{m,m}=1/h$

一般情况下，公式（9）中右端项 $b_{i}$ 的计算需要用到数值积分，数值积分公式的选取依赖于有限元空间的选取及整个数值格式的精度，通常情况下，用四阶精度的两点高斯公式就足够，即

$\int_{a}^{b}g(x)dx \approx \frac{b-a}{2}(g(\frac{a+b}{2}-\frac{b-a}{2\sqrt{3}})+g(\frac{a+b}{2}+\frac{b-a}{2\sqrt{3}}))\;\;(10)$

三、有限元编程

在编程实现过程中，可按照上述分析直接将系数矩阵中各元素 $a_{i,j}$ 的值输入，并且用数值积分计算出方程组（9）右端向量中各元素 $b_{i}$ 。虽然这种操作方式有效，但没有体现有限元的主要思想，即先化整为零，在每个小单元上考虑，再装配组合，化零为整。

原则上线性方程组（9）的系数矩阵（在结构问题中称为刚度矩阵）A的计算可以利用小的剖分单元上的单元刚度矩阵合成，线性方程组右端的向量（称为总荷载） $\mathbf{b}$ 的计算也是利用晓得剖分单元上的单元荷载合成的。

首先，在每个小单元上单独考虑单元刚度矩阵及单元荷载。具体过程为：

记小单元为 $e_{i}=[x_{i},x_{i+1}],i=0,1,\cdot\cdot\cdot,m-1$ 。其中第i个单元的左右节点 $x_{i},x_{i+1}$ 的整体编号是i和i+1，其局部编号是0和1，然后在每个小单元上考虑这个独立的单元对整个系统式（7）的贡献。不妨考虑在第k个单元 $e_{k}=[x_{k},x_{k+1}]$ 上的情况。首先考虑这个单元上的单元刚度矩阵。注意到公式（7）的左端可以写成

$\int_{0}^{1}(\sum_{i=0}^{m}u_{i}\varphi^{'}_{i}(x))\varphi^{'}_{j}(x)dx=\sum_{k=0}^{m-1}\int_{x_{k}}^{x_{k+1}}(\sum_{i=0}^{m}u_{i}\varphi^{'}_{i}(x))\varphi^{'}_{j}(x)dx$

因此，单元 $e_{k}$ 对它的贡献是

$S_{k}:=\int_{x_{k}}^{x_{k+1}}(\sum_{i=0}^{m}u_{i}\varphi^{'}_{i}(x))\varphi^{'}_{j}(x)dx=\int_{x_{k}}^{x_{k+1}}(u_{k}\varphi^{'}_{k}(x)+u_{k+1}\varphi^{'}_{k+1}(x))\varphi^{'}_{j}(x)dx$

而且只有当j=k和j=k+1时才有实质贡献，其他情况 $S_{k}$ 均为零，从而可以认为对整个公式（7）左端无实质贡献。于是，单元 $e_{k}$ 对整个公式（7）左端的贡献为：

$\begin{pmatrix} \int_{x_{k}}^{x^{k+1}}\varphi^{'}_{k}(x)\varphi^{'}_{k}(x)dx & \int_{x_{k}}^{x^{k+1}}\varphi^{'}_{k+1}(x)\varphi^{'}_{k}(x)dx\\ \int_{x_{k}}^{x^{k+1}}\varphi^{'}_{k}(x)\varphi^{'}_{k+1}(x)dx & \int_{x_{k}}^{x^{k+1}}\varphi^{'}_{k}(x)\varphi^{'}_{k+1}(x)dx \end{pmatrix}\begin{pmatrix} u_{k}\\ u_{k+1} \end{pmatrix}$

实际计算可知

$\int_{x_{k}}^{x_{k+1}}\varphi^{'}_{k}(x)\varphi^{'}_{k+1}(x)dx=-\frac{1}{h},\int_{e_{k}}\varphi^{'2}_{k}(x)dx=\int_{e_{k}}\varphi^{'2}_{k+1}(x)dx=\frac{1}{h}$

从而得到第k个单元上的单元刚度矩阵为

$\frac{1}{h}\begin{pmatrix} 1 & -1\\ -1 & 1 \end{pmatrix}$

然后考虑第k个单元上的单元荷载。公式（7）右端可以写成

$\int_{0}^{1}f(x)\varphi_{j}(x)dx=\sum_{k=0}^{m-1}\int_{x_{k}}^{x_{k+1}}f(x)\varphi_{j}(x)dx$

因此单元 $e_{k}$ 对它的贡献就是

$\int_{x_{k}}^{x_{k+1}}f(x)\varphi_{j}(x)dx$

而且仅当j=k和j=k+1时才有实质贡献。于是，单元 $e_{k}$ 对整个公式（7）右端作出如下贡献：

$\begin{pmatrix} \int_{x_{k}}^{x_{k+1}}f(x)\varphi_{k}(x)dx\\ \int_{x_{k}}^{x_{k+1}}f(x)\varphi_{k+1}(x)dx \end{pmatrix}$

里面的积分一般要借助数值积分公式（10）。这样得到第k个单元上的单元荷载为

$\begin{pmatrix} \int_{x_{0}}^{x_{1}}f(x)\varphi_{0}(x)dx\\ \int_{x_{0}}^{x_{1}}f(x)\varphi_{1}(x)dx \end{pmatrix}$

对其他单元进行相同的分析，可以得到下表：

单元	$e_{0}$	$e_{1}$	$e_{2}$	$\cdots$	$e_{k}$	$\cdots$	$e_{m-1}$
整体编号	0 1	1 2	2 3		k k+1		m-1 m
局部编号	0 1	0 1	0 1		0 1		0 1
*单元刚度矩阵（1/h）**	$\begin{pmatrix} 1 & -1\\ -1 & 1 \end{pmatrix}$	$\begin{pmatrix} 1 & -1\\ -1 & 1 \end{pmatrix}$	$\begin{pmatrix} 1 & -1\\ -1 & 1 \end{pmatrix}$		$\begin{pmatrix} 1 & -1\\ -1 & 1 \end{pmatrix}$		$\begin{pmatrix} 1 & -1\\ -1 & 1 \end{pmatrix}$
单元荷载	$\begin{pmatrix}f_{0}^{0}\\ f_{1}^{0} \end{pmatrix}$	$\begin{pmatrix}f_{0}^{1}\\ f_{1}^{1} \end{pmatrix}$	$\begin{pmatrix} f_{0}^{2}\\ f_{1}^{2} \end{pmatrix}$		$\begin{pmatrix} f_{0}^{k}\\ f_{1}^{k} \end{pmatrix}$		$\begin{pmatrix} f_{0}^{m-1}\\ f_{1}^{m-1} \end{pmatrix}$

其中：

$\begin{pmatrix} f_{0}^{k}\\ f_{1}^{k} \end{pmatrix}=\begin{pmatrix} \int_{x_{k}}^{x_{k+1}}f(x)\varphi_{k}(x)dx\\ \int_{x_{k}}^{x_{k+1}}f(x)\varphi_{k+1}(x)dx \end{pmatrix},0\leqslant k\leqslant m-1$

最后将单元刚度矩阵及单元荷载合成，如下图所示，按照节点的整体编号合成为对应的总刚度矩阵A及总荷载b。

可见，合成的刚度矩阵为

$A=\frac{1}{h}\begin{pmatrix} 1 & -1 & & & & & \\ -1 & 2 & -1 & & & 0 & \\ & -1 & 2 & -1 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -1 & 2 & -1 & \\ & 0 & & & -1 & 2 & -1\\ & & & & & -1 & 1 \end{pmatrix}$

及总荷载为

$\mathbf{b}=\begin{pmatrix} \int_{x_{0}}^{x_{1}}f(x)\varphi_{0}(x)dx\\ \int_{x_{0}}^{x_{1}}f(x)\varphi_{1}(x)dx+\int_{x_{1}}^{x_{2}}f(x)\varphi_{1}(x)dx\\ \int_{x_{1}}^{x_{2}}f(x)\varphi_{2}(x)dx+\int_{x_{2}}^{x_{3}}f(x)\varphi_{2}(x)dx\\ \vdots\\ \int_{x_{m-3}}^{x_{m-2}}f(x)\varphi_{m-2}(x)dx+\int_{x_{m-2}}^{x_{m-1}}f(x)\varphi_{m-2}(x)dx\\ \int_{x_{m-2}}^{x_{m-1}}f(x)\varphi_{m-1}(x)dx+\int_{x_{m-1}}^{x_{m}}f(x)\varphi_{m-1}(x)dx\\ \int_{x_{m-1}}^{x_{m}}f(x)\varphi_{m}(x)dx \end{pmatrix}=\begin{pmatrix} \int_{x_{0}}^{x_{1}}f(x)\varphi_{0}(x)dx\\ \int_{x_{0}}^{x_{2}}f(x)\varphi_{1}(x)dx\\ \int_{x_{1}}^{x_{3}}f(x)\varphi_{2}(x)dx\\ \vdots\\ \int_{x_{m-3}}^{x_{m-1}}f(x)\varphi_{m-2}(x)dx\\ \int_{x_{m-2}}^{x_{m}}f(x)\varphi_{m-1}(x)dx\\ \int_{x_{m-1}}^{x_{m}}f(x)\varphi_{m}(x)dx \end{pmatrix}\;\;\;(11)$

其结果与用公式（8）直接计算的结果完全一致。

在实际计算中，由于 $u_{h}\in V_{h}$ 边界条件已知，所以待求量 $\mathbf{u}=(u_{0},u_{1},\cdots,u_{m})^{T}$ 中已知 $u_{0}=u_{m}=0$ ，只需要求出其它 $u_{i},1\leqslant i\leqslant m-1$ ，为此修改原线性方程组 $A\mathbf{u}=\mathbf{b}$ 为

$\begin{pmatrix} 1 & 0 & & & & & \\ 0 & 2 & -1 & & & 0 & \\ & -1 & 2 & -1 & & & \\ & & \ddots & \ddots & \ddots & & \\ & & & -1 & 2 & -1 & \\ & 0 & & & -1 & 2 & -1\\ & & & & & 0 & 1 \end{pmatrix}\begin{pmatrix} u_{0}\\ u_{1}\\ u_{2}\\ \vdots\\ u_{m-2}\\ u_{m-1}\\ u_{m} \end{pmatrix}=\begin{pmatrix} 0\\ h\int_{x_{0}}^{x_{2}}f(x)\varphi_{1}(x)dx\\ h\int_{x_{1}}^{x_{3}}f(x)\varphi_{2}(x)dx\\ \vdots\\ h\int_{x_{m-3}}^{x_{m-1}}f(x)\varphi_{m-2}(x)dx\\ h\int_{x_{m-2}}^{x_{m}}f(x)\varphi_{m-1}(x)dx\\ 0 \end{pmatrix}$

即第0行和第0列只保留 $a_{0,0}=1$ ，其余元素均为零。第m行和第m列也只保留 $a_{m,m}=1$ ，其余元素也均为零。最后将右端向量的第一个和最后一个分量改为0即可。实质上，除去已知的 $u_{0}=u_{m}=0$ ，公式（11）可以看成是m-1阶的线性方程组。修改上述系数矩阵而不直接去除 $u_{0}$ 和 $u_{m}$ 就是为了编程时下标从0开始算起。这样采用C++编程就不易出错了。

因此，最终通过求解线性方程组（11）就可以得到所有的 $u_{i},0\leqslant i\leqslant m$ ，进而由 $u_{h}(x)=\sum_{i=0}^{m}u_{i}\varphi_{i}(x)$ 得到定义在整个区间[0,1]上的数值解。

四、编程算例

有限元法求解常微分方程两点边值问题：

$\left\{\begin{matrix} -u^{''}=((16\pi^{2}-4)sin(4\pi x)-16\pi cos(4\pi x))e^{2x},0<x<1,\\ u(0)=u(1)=0 \end{matrix}\right.$

已知精确解为 $u(x)=e^{2x}sin(4\pi x)$ 。分别取步长h=1/16和h=1/32，输出在节点 $x=\frac{i}{8},1\leqslant i\leqslant 7$ 处的数值解，并给出误差。

4.1 C++代码

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

int **lnd;  //存放单元的节点编号
double h;  //步长
double *xco,*u;  //xco存放节点的坐标，u存放节点的数值解
double pi=3.14159265359;  //圆周率π
double gauss=0.5773502692;  //数值积分中的高斯点

int main(int argc, char *argv[])
{
        int i,j,k,m,t,e,row,coln,elem;
        double ea[2][2],alpha[2],g[2],sum1,sum2;
        double **a,*b,*a1,*b1,*c1;
        double f(double x);
        double phi(int i, double x);
        double fun1(int i, double x);
        double fun2(int i, double x);
        double fun3(int i, double x);
        double exact(double x);
        double dxexact(double x);
        double integral(double a, double b, int i, double(*fun)(int, double));
        double *chase_algorithm(int i, double *a, double *b, double *c, double *d);

        m=16;  //总剖分数
        elem=m;  //总的单元数
        h=1.0/m;  //网格尺度
        printf("m=%d.\n",m);

        xco=(double*)malloc(sizeof(double)*(m+1));   //为节点坐标动态分配内存，存放在数组xco[]中
        for(i=0;i<=m;i++)
                xco[i]=i*h;   //节点坐标

        lnd=(int**)malloc(sizeof(int*)*elem);
        for(e=0;e<elem;e++)
                lnd[e]=(int*)malloc(sizeof(int)*2);

        //为单元节点编号动态分配内存，存放在二维数组lnd[]中，lnd[e][i]表示第e个单元第i个节点，e是整体单元编号，i是局部节点编号，lnd[e][i]是整体节点编号
        for(e=0;e<elem;e++)
        {
                lnd[e][0]=e;
                lnd[e][1]=e+1;
        }   //编号为e的单元上的第0个节点编号为e，第1个节点编号为e+1

        a=(double**)malloc(sizeof(double*)*(m+1));  //动态分配内存，二维数组a[][]存放总刚度矩阵，即线性方程组系数矩阵
        for(i=0;i<=m;i++)
                a[i]=(double*)malloc(sizeof(double)*(m+1));
        b=(double*)malloc(sizeof(double)*(m+1));  //一维数组b[]存放总荷载，即线性方程组右端项
        
        for(i=0;i<=m;i++)
        {
                for(j=0;j<=m;j++)
                        a[i][j]=0;
                b[i]=0;
        }   //初始化系数矩阵及右端项

        //计算单元刚度矩阵
        alpha[0]=-1.0/h;  //phi0'(x)
        alpha[1]=1.0/h;   //phi1'(x)
        for(i=0;i<2;i++)
                for(j=0;j<2;j++)
                        ea[i][j]=alpha[i]*alpha[j]*h;

        for(e=0;e<elem;e++)
        {
                i=lnd[e][0];
                j=lnd[e][1];
                //利用两点高斯积分公式计算单元荷载
                for(k=0;k<=1;k++)
                        g[k]=integral(xco[i], xco[j], lnd[e][k], fun1);  //放置单元荷载

                //合成整体刚度矩阵
                for(i=0;i<2;i++)
                {
                        for(j=0;j<2;j++)
                        {
                                row=lnd[e][i];  //确定整体行编号以明确合成总刚度矩阵时的位置
                                coln=lnd[e][j];  //确定整体列编号以明确合成总刚度矩阵时的位置
                                a[row][coln]=a[row][coln]+ea[i][j];
                        }
                        k=lnd[e][i];  //确定右端项行编号以明确合成总荷载时的位置
                        b[k]=g[i]+b[k];  //合成总荷载即整体右端项
                }
        }

        //修改边界条件
        for(t=0;t<=m;t++)
        {
                a[t][0]=0;
                a[0][t]=0;
                a[t][m]=0;
                a[m][t]=0;
        }
        a[0][0]=1;
        a[m][m]=1;
        b[0]=0;
        b[m]=0;

        a1=(double*)malloc(sizeof(double)*(m+1));  //存储矩阵a[][]中的下次对角线元素
        b1=(double*)malloc(sizeof(double)*(m+1));  //存储矩阵a[][]中的主对角线元素
        c1=(double*)malloc(sizeof(double)*(m+1));  //存储矩阵a[][]中的上次对角线元素
        a1[0]=0.0;
        c1[m]=0.0;
        for(i=0;i<=m;i++)
        {
                if(i!=0)
                        a1[i]=a[i][i-1];
                if(i!=m)
                        c1[i]=a[i][i+1];
                b1[i]=a[i][i];
        }
        for(i=0;i<=m;i++)
                free(a[i]);
        free(a);

        u=chase_algorithm(m+1, a1, b1, c1, b);  //追赶法求解三对角线性方程组
        free(a1);free(b1);free(c1);free(b);

        k=m/8;
        for(i=k;i<m;i=i+k)
                printf("xco[%d]=%.4f, u[%d]=%f, err=%.4e.\n",i,xco[i],i,u[i],fabs(u[i]-exact(xco[i])));
        //计算数值解与精确解在范数下的误差
        sum1=0.0;
        sum2=0.0;
        for(e=0;e<elem;e++)
        {
                i=lnd[e][0];
                j=lnd[e][1];
                sum1=sum1+integral(xco[i],xco[j],i,fun2);
                sum2=sum2+integral(xco[i],xco[j],i,fun3);
        }
        printf("||u-uh||=%f,||(u-uh)'||=%f.\n",sqrt(sum1),sqrt(sum2));
        for(e=0;e<elem;e++)
                free(lnd[e]);
        free(lnd);free(xco);free(u);

        return 0;
}



double f(double x) //右端项
{
        double z;
        z=((16*pi*pi - 4)*sin(4*pi*x) - 16*pi*cos(4*pi*x))*exp(2*x);
        return z;
}
double phi(int i, double x)
{
        double temp,z;
        temp=fabs(x-xco[i]);
        if(temp<=h)
                z=1.0-temp/h;
        else
                z=0.0;
        return z;
}
double exact(double x)  //精确解
{
        return sin(4*pi*x)*exp(2*x);
}
double dxexact(double x)  //u(x)的导函数
{
        return (2*pi*cos(4*pi*x)+sin(4*pi*x))*exp(2*x)*2;
}
double fun1(int i, double x)  //算单元载荷时的被积函数
{
        return f(x)*phi(i,x);
}
double fun2(int i, double x)  //算||u-uh||^2时的被积函数
{
        double temp,z;
        temp=u[i]*phi(i,x)+u[i+1]*phi(i+1,x);
        z=exact(x)-temp;
        return z*z;
}
double fun3(int i, double x)  //算||(u-uh)'||^2时的被积函数
{
        double temp, z;
        temp=(u[i+1]-u[i])/h;
        z=dxexact(x)-temp;
        return z*z;
}
double integral(double a, double b, int i, double(*fun)(int, double))  //在区间[a,b]上对被积函数fun(i,x)进行数值积分(两点高斯公式)
{
        double mid,w,ans;
        mid=(b+a)/2.0;
        w=(b-a)/2.0;
        ans=w*((*fun)(i,mid+w*gauss)+(*fun)(i,mid-w*gauss));
        return ans;
}
double *chase_algorithm(int n, double *a, double *b, double *c, double *d)
{
        int i;
        double *ans,*g,*w,p;
        ans=(double*)malloc(sizeof(double)*n);
        g=(double*)malloc(sizeof(double)*n);
        w=(double*)malloc(sizeof(double)*n);
        g[0]=d[0]/b[0];
        w[0]=c[0]/b[0];

        for(i=1;i<n;i++)
        {
                p=b[i]-a[i]*w[i-1];
                g[i]=(d[i]-a[i]*g[i-1])/p;
                w[i]=c[i]/p;
        }
        ans[n-1]=g[n-1];
        i=n-2;
        do
        {
                ans[i]=g[i]-w[i]*ans[i+1];
                i=i-1;
        }while(i>=0);

        free(g);free(w);

        return ans;
}

4.2 计算结果

（1）当h=1/16时，结果如下：

m=16.
xco[2]=0.1250, u[2]=1.284629, err=6.0386e-04.
xco[4]=0.2500, u[4]=0.000729, err=7.2911e-04.
xco[6]=0.3750, u[6]=-2.116903, err=9.6762e-05.
xco[8]=0.5000, u[8]=0.000253, err=2.5350e-04.
xco[10]=0.6250, u[10]=3.492002, err=1.6593e-03.
xco[12]=0.7500, u[12]=0.001764, err=1.7641e-03.
xco[14]=0.8750, u[14]=-5.754793, err=1.9041e-04.
||u-uh||=0.139331,||(u-uh)'||=7.796860.

（2）当h=1/32时，结果如下：

m=32.
xco[4]=0.1250, u[4]=1.284062, err=3.6750e-05.
xco[8]=0.2500, u[8]=0.000044, err=4.4014e-05.
xco[12]=0.3750, u[12]=-2.116995, err=5.3512e-06.
xco[16]=0.5000, u[16]=0.000015, err=1.5303e-05.
xco[20]=0.6250, u[20]=3.490444, err=1.0098e-04.
xco[24]=0.7500, u[24]=0.000106, err=1.0650e-04.
xco[28]=0.8750, u[28]=-5.754616, err=1.2827e-05.
||u-uh||=0.035184,||(u-uh)'||=3.909623.