有限差分法求解偏微分方程

自然科学与工程技术中种种运动发展过程与平衡现象各自遵守一定的规律。这些规律的定量表述一般地呈现为关于含有未知函数及其导数的方程。我们将只含有未知多元函数及其偏导数的方程，称之为偏微分方程。初始条件和边界条件称为定解条件，未附加定解条件的偏微分方程称为泛定方程。对于一个具体的问题，定解条件与泛定方程总是同时提出。定解条件与泛定方程作为一个整体，称为定解问题。

在介绍有限差分法求解偏微分方程的过程中，我们会用到Jacobi迭代法与Gauss-Seidel迭代法的相关内容，如果读者对此不甚了解，可以参阅下文：

Jacobi迭代法与Gauss-Seidel迭代法

欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji，为保证公式、图表得以正确显示，强烈建议你从该地址上查看原版博文。本博客主要关注方向包括：数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。

椭圆方程

由于大多数工程问题都是二维问题，所以得到的微分方程一般都是偏微分方程，对于一维问题得到的是常微分方程，解法与偏微分方程类似，故为了不是一般性，这里只讨论偏微分方程。由于工程中高阶偏微分较少出现，所以本文仅仅给出二阶偏微分方程的一般形式，对于高阶的偏微分，可进行类似地推广。二阶偏微分方程的一般形式如下：
$A Φ x x + B Φ x y + C Φ y y = f (x, y, Φ, Φ x, Φ y)$
其中， Φ 表示一个连续函数。当 A,B,C 都是常数时，上式成为准线性，有三种准线性方程形式：

如果 Δ=B2−4AC<0 ，则称为椭圆型方程；
如果 Δ=B2−4AC=0 ，则称为抛物型方程；
如果 Δ=B2−4AC>0 ，则称为双曲型方程；

椭圆方程是工程技术应用中所涉及的偏微分方程里最为普遍的一种形式。根据椭圆方程的具体形式又可以将其分为以下三种形式：

拉普拉斯（Laplace）方程： ∇2u=0 ；
泊松（Poisson）方程： ∇2u=g(x,y) ；
亥姆霍兹（Helmholtz）方程： ∇2u+λu=0 。

其中 u 是关于 x 和 y 的二元函数。

有限差分法

差分方法又称为有限差分方法或网格法，是求偏微分方程定解问题的数值解中应用最广泛的方法之一。它的基本思想是：先对求解区域作网格剖分，将自变量的连续变化区域用有限离散点（网格点）集代替；将问题中出现的连续变量的函数用定义在网格点上离散变量的函数代替；通过用网格点上函数的差商代替导数，将含连续变量的偏微分方程定解问题化成只含有限个未知数的代数方程组（称为差分格式）。如果差分格式有解，且当网格无限变小时其解收敛于原微分方程定解问题的解，则差分格式的解就作为原问题的近似解（数值解）。因此，用差分方法求偏微分方程定解问题一般需要解决以下问题：

选取网格；
对微分方程及定解条件选择差分近似，列出差分格式；
求解差分格式；
讨论差分格式解对于微分方程解的收敛性及误差估计。

下面我们就以拉普拉斯方程的数值解法为例来演示一下有限差分法的基本思路。我们首先写出完整的拉普拉斯方程如下：
$\partial 2 f \partial x 2 + \partial 2 f \partial y 2 = 0$
现在的问题其实是要求我们在一个给定的二维区域中求解满足方程的每一点 (x,y) 。一些区域中的点将被用来给出边界条件（hold boundary conditions）。

于是我们将整个二维区域离散化成若干个点，如下图所示为其中的五个相邻点：

根据偏导数的定义则有：

\partial f \partial x \approx 1 Δ [f (x + Δ 2, y) - f (x - Δ 2, y)]

\partial 2 f \partial x 2 \approx 1 Δ {1 Δ [f (x + Δ, y) - f (x, y)] - 1 Δ [f (x, y) - f (x - Δ, y)]} = 1 Δ 2 [f (x + Δ, y) - 2 f (x, y) + f (x - Δ, y)]

同理可得：

\partial 2 f \partial y 2 \approx 1 Δ 2 [f (x, y + Δ) - 2 f (x, y) + f (x, y - Δ)]

将上述两个结果带入拉普拉斯方程可得：

1 Δ 2 [f (x + Δ, y) + f (x - Δ, y) + f (x, y + Δ) + f (x, y - Δ) - 4 f (x, y)] = 0

方程组求解

回想雅各比迭代法（可以参考本文最开始给出的文章链接），假设我们有一个由 n 个线性方程组成的系统（也就是线性方程组）：
$A x = b$
那么Jacobi迭代可以描述为：
$x k i = 1 a i , i [b i - \sum j \neq i a i, j x k - 1 j]$
其中上标 k 表示第 k 轮迭代。

注意我们在上一小节最后得出的拉普拉斯方程离散化形式给出了（离散化后）区域上众多点中的一个点的求解方程，所有点的求解方程合在一起就构成了一个大的方程组。我们把求解某点 (x,y) 的方程重写成Jacobi迭代的形式，则有：
$f k (x, y) = 1 4 [f k - 1 (x - Δ, y) + f k - 1 (x, y - Δ) + f k - 1 (x + Δ, y) + f k - 1 (x, y + Δ)]$
重复应用上述迭代式，最后方程就会收敛到解的附近。
本来连续的一个区域经过离散化处理之后就变成了一个网格结构，假设网格的大小是 n2 ，它们的标签为 x1,x2,⋯,xn2 ，如下图所示

上面这种自然排列的点序可以得出不超过

n2 个五元线性方程：

x i - n + x i - 1 - 4 x i + x i + 1 + x i + n = 0

注意我们不对处于边界上的点（例如：

x1,x2,⋯,x4,x5,x8,x9,⋯ ）应用上述方程。最后我们将得到如下所示的一个大型的稀疏线性方程组。

有限差分法求解偏微分方程

除了Jacobi迭代之外，（如果你看了本文最开始推荐的文章也应该知道）我们还可以采用Guass-Siedel迭代来加速方程组解的收敛速度。Guass-Siedel迭代的一般形式为：

x k i = 1 a i , i [b i - \sum j = 1 i - 1 a i, j x k j - \sum j = i + 1 n a i, j x k - 1 j]

此时拉普拉斯方程需用下式进行求解（我们不再做详细讨论）：

f k (x, y) = 1 4 [f k (x - Δ, y) + f k (x, y - Δ) + f k - 1 (x + Δ, y) + f k - 1 (x, y + Δ)]

一个讨论：这些内容有什么用？

我们前面提到偏微分方程在工程中有重要应用。但是在信息技术中有没有离我们比较近的应用实例呢？事实上，偏微分方程在图像处理中就有重要应用！本文主要是以拉普拉斯方程的数值解为例来讨论的，而本文前面我们也提到过椭圆方程中除了拉普拉斯方程之外，还有一类叫做泊松方程。图像处理中基于泊松方程的算法构成了一大类的具有广泛应用的算法，可以用于图像融合、图像去雾、图像拼接等等。例如下图就是基于解泊松方程的方法实现的图像泊松编辑的效果图：

更多关于泊松融合、泊松编辑方面的内容还可以参考如下链接：
http://blog.csdn.net/baimafujinji/article/details/50583497

更多关于数学在图像处理方面的应用，或者图像处理中所需的数学知识和数学原理，你还可以参考我的新书《图像处理中的数学修炼》

参考文献与推荐阅读材料

【1】张文生. 科学计算中的偏微分方程有限差分法，高等教育出版社，2006.
【2】左飞. 图像处理中的数学修炼，清华大学出版社，2016.
【3】百度文库：有限差分法求解偏微分方程，最后浏览时间为2016年11月.

秒客网