4.3.1 有限差分法的原理及应用
从以上几种数值解法的介绍上看，FDM的出现早于FVM和FEM，相对而言更简单易学，FDM中网格化和离散化的思想也是FVM和FEM的重要特征；从实际应用上看，在工程设计和科学研究中，FDM目前依然是常用的一种方法。因此，学习FDM不仅有利于我们学习和理解FVM和FEM，还具有一定的实用性。需要说明的是，本书主要介绍微流控芯片技术的建模分析方法，偏微分方程和数值解法是相应基础知识。在后续内容中，主要以介绍COMSOL Multiphysics软件的应用为主，该软件主要采用FEM来求解偏微分方程。至于FVM，本书不做具体介绍。接下来，我们从FDM开始，了解如何求解偏微分方程。
针对一个二维问题，我们采用简单的正方形网格来作为FDM的计算网格，如图4.1所示。将方格作为计算网格的好处是节点位置序号和其空间分布天然契合，无须单独考虑网格节点的位置信息。
图4.1 FDM计算网格
图4.1中，矩形区域中划分出若干个方格，设方格的边长为h，坐标(i，j)表示中心点的位置坐标。进一步地，周围点的坐标也可以相应地表示出来。
假设采用FDM考虑温度的稳态扩散问题，u为温度变量。首先考虑x方向，在i+1、i-1点分别对(i，j)点的温度变量进行泰勒展开，可以得到如下形式：
式(4.8)与式(4.9)的区别在于右侧奇数项符号相同，偶数项符号相反。从右侧第三项开始舍弃高阶项，可得：
基于式(4.10)、式(4.11)，通过移项，可得在x方向的前、后差分公式：
同理，在y方向，也可得到相应一阶导数的前、后差分公式：
以上差分格式只具有一阶精度，而中心差分法可以提高精度至二阶。将式(4.8)与式(4.9)左右两端分别作差，并舍去高阶项，即可得一阶导数的中心差分格式：
若将式(4.8)与式(4.9)左右两端分别作和，则可以得到二阶偏导数的差分格式：
类似地，在y方向为
利用以上二阶偏导数的差分格式，可以得到拉普拉斯方程的“五点差分格式”为
即
如此，拉普拉斯方程便可以转化为一个线性方程。这样的代数方程就可以通过计算机编程来求解了。下面举一个具体的例子，看看如何用“五点差分格式”来编程求解拉普拉斯方程。
例4.2 在区域Ω={(x，y)|0≤x≤4，0≤y≤4}内，边界条件为
求解偏微分方程
图4.2 网格划分与边界条件设置
解：将求解区域划分成如图4.2所示的5×5方格网格，为方便列写方程组，在边界点处标注相应的数值。U1～U9是在相应网格节点上需要求解的变量。
下面，根据“五点差分格式”列方程组。易知，U1～U9都满足拉普拉斯方程。
由式(4.18)，对于U1可得U4+U2+20+80-4U1=0，合并移项可得4U1-U2-U4=100。同理，也可得U2～U9对应的其他8个线性方程。因此，9个未知数对应一个具有9个独立方程的线性方程组，在理论上可以求出U1～U9。将这9个方程转化为系数矩阵形式：
求解以上线性方程组可得：
U=[U1 U2 U3 U4 U5 U6 U7 U8 U9]
=[55.7 43.2 27.1 79.6 70.0 45.4 112.9 111.8 84.3]
用图示法表示解的分布(角点采用相邻节点进行线性插值)，如图4.3所示。
图4.3 求解域内的解分布(5×5网格)(书后附彩插)
当网格划分得比较粗时，采用差商代替导数会引起比较大的数值误差，因此在实际操作中，所采用的网格划分需要满足精度要求。所谓的精度要求，就是求解的结果与网格划分的尺度接近无关。如何证明求解的结果与网格划分的尺度接近无关呢?具体来说，就是对同一个求解域由粗到细进行网格划分，然后比较求解结果的变化。若当前求解的结果与进一步细分网格获得的结果之间的差值小于设定值，则认为当前的网格划分尺度是合适的，而不再采用进一步细分的网格尺度。这是因为，尽管通过进一步细分网格可以获得更高的精度，但是会导致更大的计算量，这相对于所获得的精度提升，计算量的增加所导致的低效率求解过程可谓得不偿失。例如，将求解域划分为5×5、11×11、21×21、41×41四种尺度，相应网格的边长分别为
，那么所需求解的线性方程组个数分别为(5-2)×(5-2)=9、(11-2)×(11-2)=81、(21-2)×(21-2)=361、(41-2)×(41-2)=1 521。可见，随着网格变细，所需求解的方程将越来越多，因此合适的网格规模是必要的。
对于例4.2，若划分成9×9的网格，就需要求解(9-2)×(9-2)=49个方程，求解结果如图4.4所示。
图4.4 求解域内的解分布(9×9网格)(书后附彩插)
对照图4.3中所示结果，易知图4.4中的解分布更加平滑，这说明加密网格有利于求解精度的提升。
例4.2采用“五点差分格式”求解拉普拉斯方程时，存在一个问题：当网格数量增大，就会存在非常大的稀疏系数矩阵(结合式(4.20)，可以想象具有成千上万个方程的情况)，这样会占用较多内存。因此，迭代法成为解决稀疏矩阵方程组求解的高效方法。
D.M.Young于20世纪70年代提出了逐次超松弛(Successive Over Relaxation)迭代法，简称SOR方法，这是一种经典的迭代算法。它是为了求解大规模系统的线性方程组提出来的、在Gauss-Seidel迭代法基础上为提高收敛速度并采用加权平均而得到的新算法。由于超松弛迭代法的公式简单，容易编制程序求解，因此很多工程学、计算数学中都会应用超松弛迭代法。对于逐次超松弛迭代法，松弛因子的选取对算法的收敛速度有较大影响。通常，对于方程组A x=B，若A为正定矩阵，则当0＜ω＜2时，逐次超松弛迭代恒收敛[3]。
相应的迭代格式如下：
采用SOR方法，设置31×31和61×61网格规模求解例4.2的MATLAB代码如下：
图4.5 求解域内的解分布(书后附彩插)
(a)网格规模为31×31；(b)网格规模为61×61
从图4.5可见，定性比较下两种网格下的解分布几乎相同。但在数值模拟方法中，对于结果不能以定性比较下结论。在此选择x方向中点对应剖面上的求解结果进行比较，得到的结果如图4.6所示。
图4.6 x方向中点对应剖面上的求解结果比较
由图4.6可见，在网格规模为31×31、61×61下的数值解基本重合，说明两种网格规模下的解均具有几乎相同的精度。从选择网格规模的角度，由于31×31网格规模的计算量更小，因此可将其作为求解例4.2的网格。
以上简单介绍和示例了FDM在求解偏微分方程中的应用。总结起来，可知采用FDM求解偏微分方程的主要步骤是：首先，将微分形式离散为差分格式；然后，构建线性方程组；最后，采用求解线性方程组的方式得到偏微分方程的解。FDM的原理为FEM的学习提供了一些参考。首先，对于相同的数值问题，不论是采用FDM还是采用FEM，得到的解结果应相同，因为问题本身的答案是唯一和确定的，与采用何种方法来解决这个问题无关；其次，FDM中的网格划分、局部离散、求解线性方程组等思路，在FEM中也有类似的应用。但是，在线性方程组的构造方法上，FDM和FEM大有不同，接下来将重点介绍FEM求解偏微分方程的思路。