加星标，才能不错过每日推送！方法见文末动图
函数迭代是数学中一个非常重要和有趣的主题，它在不同的领域有着不同的应用和着眼点。在动力系统中，函数迭代可以揭示复杂系统的演化规律和混沌现象；在计算数学中，函数迭代可以求解各种非线性方程和优化问题。关于动力系统和混沌领域里的迭代函数，《返朴》之前已刊发数篇相关文章，本文则是从计算数学的角度，重新向读者介绍函数迭代的基本内容。
我在最近撰写的《这么说迭代，你一定能懂》以及几篇后续文章中，介绍了函数迭代的基本概念，展示了倍周期分叉的美妙图形，历经了周期三点引爆出的“天下大乱”，并运用高中代数、大学初等微积分加上几何图形的直观，帮助读者理解从有序到混沌沿途各种现象之“所以然”，从中领略数学推理在发现自然规律的过程中的奇妙作用。
然而，上述几篇数学科普文章的立足点位于“动力系统”的地盘，或更具体地说是在它的子领域“离散动力系统”。离散动力系统这门学科的主要任务就是研究“函数迭代”。自然，这里的“函数”是在最广泛的意义上理解的，即它的定义域可以是任何“空间”——带有某个数学结构的一类集合。当然，最简单的空间就是实数轴上的单位区间或坐标平面上的单位圆周，即便对于这样“简单”的空间所对应的“一维离散动力系统”，其主题却可以写成一本有五百页的大书，需要的数学语言来自分析、拓扑、几何、代数等几乎所有现代数学的分支。
在离散动力系统的范围谈论函数迭代，讨论的重点是不动点和周期点的稳定性问题，探索的目光聚焦在迭代点轨道的最终性态，而混沌学家则更感兴趣于最终走向不可预测的不规则行为，并由此进一步用概率统计的工具研究所有这些行踪诡异的轨道的整体性质。这些都是我在前几篇文章中试图向读者介绍的基本概念和知识。
然而，在出发点和目的性完全不一样的另一个极其有用的数学领域——计算数学里，“函数迭代”同样扮演着不可或缺的角色。“计算数学”，顾名思义，就是以计算为特征的数学，它设计可行、快速、有效、稳定的算法，通过电子计算机的程序实现，数值求解科学技术中建立起来的各类方程：常微分方程、偏微分方程、积分方程或一般的算子方程。离散化这些“连续性方程”就得到线性或非线性的代数方程组，求解这些方程组的一个主要方法就是“迭代法”。
与离散动力系统研究函数迭代的多重目的不太一样的是，计算数学研究迭代的目标很一致：寻求迭代轨道的收敛性及其速率。为此，计算数学领域里的学者谈到迭代时对周期点几乎是“冷若冰霜”的，因为它导致迭代法不收敛，这不是计算数学家和工程学家所乐见的。所以当我们从计算数学的观点看函数迭代时，我们的注意力只需放在迭代轨道收敛性这点上。
四十一年前，我在南京大学数学系修读计算数学最优化理论与方法硕士研究生的第一门专业基础课是“非线性方程组的迭代解”，选用的教材是马里兰大学数学系James M. Ortega和Werner C. Rheinboldt 于1970年出版的Iterative Solution of Nonlinear Equations in Several Variables（多变元非线性方程的迭代解），分两学期学完。第一学期由老师讲授课程，第二学期由学生分章报告，我们的收获很大，我甚至做了书中全部习题。五年后，在密歇根州立大学读博士学位的课程时，由于修了李天岩教授一门“[0, 1]上的遍历理论”学年课程，我对离散动力系统中的“函数迭代”产生了兴趣，第二年写的博士论文居然与计算遍历理论挂上了钩。可以说，在我求学成长的岁月中，我先后戴着计算数学和动力系统的面具，学过两次“函数的迭代”。
学过的知识有用，就没有白学。这两种“函数迭代”的理论都对我的研究论文有过贡献，所以算是做到了“学以致用”。但是现在的我多了一个想法：将自己年轻时代学来的有用知识服务于公众，传播更多的数学真理。在这篇文章里，我将向对科学计算感兴趣的人介绍迭代法背后的基本思想。与往常一样，我们主要考虑区间函数的简单情形，方便大家理解。
假设G是定义在一个区间I上的连续函数，并将I映到自身，这样在I中取定一个初始点x0，由此出发，就可以周而复始地迭代下去得到迭代点数列
x0, x1, x2, …, xn, …，
记为{xn}，其中第n个迭代点xn是G在第n-1个迭代点xn-1的值，即
xn = G(xn-1), n = 1, 2, …; 给定x0。
如果将n个G复合起来的复合函数记为Gn，则xn = Gn(x0)。由于函数G连续并将I映入I，函数Gn也连续且将I映入I。
计算数学家对函数迭代的基本要求是：迭代点数列{xn}必须收敛。在讨论这个数列何时收敛的问题前，我们先问另一个问题：假设{xn}确实收敛到I中的一点，这个极限点会有什么性质？或言之，数列{xn}的极限x*与迭代的函数G有什么关系？
注意我们已经假设G在定义域区间I上连续，这个信息加上初等微分学中的连续概念就提供了回答上述问话的线索。连续函数有一个基于数列的等价说法，即函数f在点a连续当且仅当任给f定义域内的数列{xn}，只要它收敛于a，则对应的函数值数列{f(xn)}收敛到f(a)。由迭代程式xn = G(xn-1)可知，两个数列{xn}和{G(xn-1)}是恒等的数列，既然{xn}收敛到x*，那么{G(xn-1)}也收敛到同一个极限x*。另一方面，由于数列{xn-1}和{xn}一样收敛到x*，故G在x*的连续性就推出{G(xn-1)}收敛到G(x*)。这样x*和G(x*)都是{xn}的极限，由收敛数列极限的唯一性，必有x* = G(x*)。这说明，迭代点数列{Gn(x0)}如果收敛到连续函数G的定义域中的一点，则它一定是G的不动点。这个命题的逆否命题可以帮助我们确认何时一个迭代点数列不可能收敛，那就是只要迭代函数G没有不动点。例如，如果一个人迭代自然指数函数G(x) = ex，则无论他从哪一个初始点出发，迭代点数列都不会收敛，因为不动点方程x = ex在实数集合里无解。
上述简单结论的一个逻辑后果是，计算数学中迭代法的目的是数值求解不动点方程x = G(x)。如果这个方程根本无解，那就不必构造什么迭代法求解它了，所以在研究迭代法前，要先假设G至少有一个不动点。但是，如果闭着眼睛随便选取初始点x0开始迭代，能保证{Gn(x0)}收敛吗？
我们再看一个例子：G(x) = x – x2。这时G依然只有唯一的不动点0。运用《这样讲迭代，你一定能懂》中的“图像迭代法”，或通过初等微分学中的“单调收敛定理”可以严格证明：当0 < x0 ≤ 1时，{Gn(x0)}收敛到不动点0，而当x0 < 0或x0 > 1时，{Gn(x0)}发散到负无穷大。虽然G在不动点0处导数的绝对值等于1，但对不动点右侧临近的初始点，迭代点数列收敛到该不动点。读者可以用几何方法画出一个连续函数y = G(x)的图像曲线，它与xy-直角坐标系的对角线y = x相切于一点(x*, x*)，在切点的左侧，曲线位于对角线上方，向上弯曲；在切点的右侧，曲线位于对角线下方，向下弯曲。然后你就会发现，尽管G在x*处的导数绝对值等于1，但从x*左右两边附近的任一初始点x0出发的迭代点数列{Gn(x0)}收敛到x*。
这些例子说明，在迭代函数G可求导数的前提下，当G在其不动点x*的导数绝对值等于1时，在其附近出发的迭代点数列可能收敛到x*，也可能不收敛到x*；而当G在不动点x*的导数绝对值大于1时，则根本不能指望迭代点数列会收敛到x*。因此，要想迭代法xn = G(xn-1)生成的数列收敛到G的一个不动点，我们一般需要假设G在这个不动点的导数绝对值严格小于1。虽然这个假设不是迭代法收敛的必要条件，但确实是保证收敛的一个充分条件。它的证明不难，也容易理解，但需要对导数的极限定义有透彻的理解。
许多人号称学过微积分，但可能还需要温故一下微积分学中最基础性的概念——极限，所以我们先来回忆一下函数在一点的极限的精确定义：函数f当自变量x趋向于数a时有极限L，是指对于任意给定的正数ε，存在正数δ，使得只要位于f定义域内的x减去a的绝对值大于0但又小于δ，其对应的函数值f(x)减去L的绝对值就会小于ε。用数学的符号，就是要找到δ，使得若不等式0 < |x – a| < δ成立，则不等式|f(x) – L| < ε也成立。
由于函数的导数是通过函数的平均变化率的极限来定义的，函数f在其定义域内一点a的导数f’(a) = lim x → a (f(x) – f(a))/(x-a)按照上面的 “ε-δ”语言来说就是：任给正数ε，存在正数δ，使得对f定义域内的所有x，不等式0 < |x – a| < δ推出不等式
|(f(x)–f(a))/(x-a) – f’(a)| < ε。
好了，我们用上面的导数定义来证明关于迭代法收敛性的一个基本定理：设想迭代函数G在其定义域区间内部有一个不动点x*并且在该点有导数。如果|G’(x*)| < 1，则存在以x*为中点的一个小区间(x* - δ, x* + δ)，使得任取此区间中的一点x0作为初始点，迭代点数列{Gn(x0)}中的所有点都在同一区间内并且lim n → ∞ Gn(x0) = x*。
在读懂下面的解析证明之前，不知道这个命题的读者可用图像迭代法看到迭代点的收敛性。我们现在严格证明之。首先因为x*位于G的定义域区间的内部，存在一个正数δ0使得区间(x*-δ0, x*+δ0)包含在G的定义域中。由假设条件，可取ε = (1 - |G’(x*)|)/2 > 0。对于这个正数，上述的导数定义告诉我们，存在正数δ ≤ δ0，使得满足不等式0 < |x – x*| < δ的所有点x也满足不等式
|(G(x)–G(x*))/(x-x*) – G’(x*)| < (1 - |G’(x*)|)/2。
由于G(x*) = x*，上式可以改写成
|(G(x)–x*)/(x-x*) – G’(x*)| < (1 - |G’(x*)|)/2。
运用代数不等式 |a| - |b| ≤ |a – b|，就可推出 |G(x*)-x*|/|x-x*| - |G’(x*)| < (1 - |G’(x*)|)/2。故
|G(x)-x*|/|x-x*| < |G’(x*)| + (1 - |G’(x*)|)/2 = (1 + |G’(x*)|)/2。
记r = (1 + |G’(x*)|)/2，则正数r < 1。上面的不等式即为
|G(x)-x*| < r |x-x*|， 任给(x*-δ, x*+δ)中的点x。
因为r小于1，此不等式说明，当初始点x0 属于(x*-δ, x*+δ)时，第一个迭代点x1 = G(x0)也在(x*-δ, x*+δ)之中，故同一不等式保证第二个迭代点x2 = G(x1) = G2(x0)继续落在(x*-δ, x*+δ)当中，且有|x2 – x*| < r |x1 – x*| < r2 |x0 – x*|。推而广之，我们知道所有迭代点xn 都在区间(x*-δ, x*+δ)之内，并满足不等式
|xn – x*| < rn |x0 – x*|，n = 1, 2, 3, …。
因为lim n → ∞ rn = 0，结论就是：迭代点数列{xn} 趋向于x*。
那么，有名闻遐迩且至少具有超线性收敛的迭代法吗？当然有，而且只需要初等微分学中的泰勒公式就可以构造出许多更高收敛阶的迭代法。然而一般而言，如同俗语所述，“一分耕耘一分收获”，要获得高阶的迭代法，成本也是可观的，计算数学中的“成本”是用“计算工作量”来衡量的。不过，我们几乎都听过一个大名鼎鼎的迭代法，甚至和它常打交道，这就是名字挂的是“牛顿”但实际上应该是“辛普森”的“牛顿迭代法”。为何它应该易名为“辛普森迭代法”，可见我的师兄弟曾钟刚教授和我合写的科普文章《牛顿迭代法传奇（上）：张冠李戴的命名》。
牛顿迭代法用于数值求解非线性方程，它的基本思想用几何表达最直观易懂。解方程f(x) = 0的几何意义就是找到函数y = f(x)的图像与x-轴的交点。在函数f可求导且导数值都不为零的条件下，设x*是方程的一个解，但我们却不知道它是几，于是我们先取一个初始点x0作为x*的一个猜测。自然f(x0)一般不恰好为零。在图像上的点(x0, f(x0))处作曲线的切线，根据导数的几何解释，这根切线的斜率是f’(x0)。由直线的点斜式方程可知，切线的方程是y – f(x0) = f’(x0)(x – x0)，在此方程中令y = 0而解出未知数x的值 x1 = x0 - f(x0)/f’(x0)，它就是切线与x-轴的交点的x-坐标，当x0距离精确解x*不远时应该比x0更靠近x*。对x1重复做上面的事，得到比x1离x*更近的点x2。由此类推，有了下列求解方程f(x) = 0的牛顿迭代法：
xn = xn-1 – f(xn-1)/f’(xn-1)，n = 1, 2, 3, …。
如果将原方程f(x) = 0改写成与之等价的不动点方程x = x – f(x)/f’(x)，那么上述的牛顿法就是用于后一个方程的直接迭代法 xn = G(xn-1), n = 1, 2, …, 其中G(x) = x – f(x)/f’(x)。为了研究牛顿法的收敛速率，我们需要假定f的二阶导数在x*存在，这样就保证了迭代函数G在x*处一阶导数的存在性。求导数的商法则给出
G’(x*) = 1 – [f’(x*)2 – f(x*)f’’(x*)]/f’(x*)2 = f(x*)f’’(x*)/f’(x*)2 = 0，
理由是f(x*) = 0。根据上面证明出的迭代收敛基本定理，存在以x*为中心的某个开区间(x*-δ, x*+δ)，使得只要初始点x0取之于它，牛顿法产生的迭代点数列{xn}不仅收敛到解x*，而且收敛速率是超线性的。事实上，在关于二阶导数稍微更强一点的条件下，牛顿迭代是二阶收敛的。
到目前为止讲到的迭代收敛都是局部性的，那怎么扩大保证收敛的初始点范围呢？当然，有不少途径，其中一个十分简单，但对导数的要求却颇为苛刻。既然简单，便适合在这里提及。应用这个方法需满足的条件是：将定义域区间[a, b]映到自身的迭代函数G有一个不动点x*，且导数处处存在并满足不等式|G’(x)| ≤ r，其中r为小于1的一个正数。
若条件成立，则对[a, b]中的任一初始点x0，由微分学里的拉格朗日中值定理，
G(x0) – x* = G(x0) – G(x*) = G’(c)(x0 – x*)，
其中c是位于x0和x*之间的一个数。这样，
|G(x0) – x*| = | G’(c)(x0 – x*)| = |G’(c)| |x0 – x*| ≤ r |x0 – x*|，
即G(x0)不仅位于[a, b]，而且比x0更靠近x*。用归纳法，得到
|Gn(x0) – x*| ≤ rn |x0 – x*| → 0。
因此，迭代点数列{Gn(x0)}收敛到x*。
如果我们再看一次证明就会发现，施加于导数的要求事实上是为了实现不等式
|G(x) – G(x*)| ≤ r |x – x*|，任给定义域[a, b]中的点x。
只要这个不等式满足，该迭代法对任意初始点都会收敛。然而该不等式并不要求G必须可导，甚至也不需要G在x ≠ x*处连续。当然，用于迭代的函数一般都是至少连续的，所以可以将上述不等式中的特殊点x*改为与x有同等地位的y，即
|G(x) – G(y)| ≤ r |x – y|，任给闭区间[a, b]中的点x和y，
其中G将[a, b]映到自身，且0 < r < 1。满足这些条件的G称为在区间[a, b]上是压缩的。区间上的压缩函数是处处连续的，但不一定处处可导。压缩函数一定存在不动点，因为这是“压缩”性质的一个直接推论，证明的基本思想来自公比绝对值小于1的等比级数的收敛性和实数的完备性这两个事实，细说如下：
在区间[a, b]中任取一点x0迭代G，得迭代点数列{xn}，其中xn = G(xn-1) = Gn(x0)。则对任意非负整数n有
|xn+1 – xn| = |G(xn) – G(xn-1)| ≤ r |xn – xn-1|。
累次使用上述过程，就有 |xn+1 – xn| ≤ rn |x1 – x0|。从而对任意自然数p，由代数中的三角不等式及已得的不等式，
|xn+p – xn| ≤ |xn+1 – xn| + |xn+2 – xn+1| + … + |xn+p – xn+p-1|
≤ (1 + r + … + rp) |xn+1 – xn| ≤ rn |x1 – x0|/(1-r)。
既然lim n → ∞ rn = 0，上式右端数列的极限为0，说明{xn}是一个柯西数列，即双下标数列{|xm – xn|}当m和n都趋于无穷大时的极限为0。因为在实数集内，任何柯西数列都是收敛的，故极限lim n → ∞ xn = x*存在。由于对所有n都有a ≤ xn ≤ b，极限x*属于[a, b]；由于G处处连续，x*是G的不动点。如果在几行之上的不等式中取p → ∞时的极限，我们还可以得出迭代到第n步后的“误差估计”：
|x* – xn| ≤ rn |x1 – x0|/(1-r)，
这是计算数学家和工程师最爱看到的结果。
更进一步，压缩函数不仅有不动点，而且仅有一个不动点，因为如果有相异不动点x*和y*的话，则 |x*- y*| = |G(x*) – G(y*)| ≤ r |x* - y*|，因为x* - y* ≠ 0及r < 1，这个不等式是不可能的。这个唯一性是许多其他著名的不动点定理如“布劳威尔不动点定理”所缺乏的，一维的布劳威尔不动点定理本质上就是微分学里的介值定理，它只要求函数连续，所以少了一点限制条件，不动点的个数就有可能大于一。这和家长对孩子读书的框框条条效果类似，限制越多，自由越少，子女今后的成就也就可能越少。
在数学上，“唯一性”是很重要的一条性质，因此压缩函数的概念被抽象成泛函分析中所谓“距离空间”内的“压缩映射”。这门学科的集大成者、波兰利沃夫数学学派的领袖巴拿赫 (Stefan Banach，1892-1945)，在他进入而立之年时，提出了一个“压缩映射定理”，现在该定理又以它的创立者命名为：巴拿赫不动点定理。其专业性陈述是这样的：
令(X, d)为一个完备的距离空间。若G: X → X为一压缩映射，即存在小于1的一个数r，使得
d(G(x), G(y)) ≤ r d(x, y)
对X中所有的元素x和y都成立，则G有并只有一个不动点x*，且从X中任一初始点x0出发的迭代点列{Gn(x0)}收敛到x*，其第n个迭代点的逼近上界为
d(x*, xn) ≤ rn d(x1, x0)/(1-r)。
这几个抽象术语自然需要解释，但只要发现上面定理的内容，包括条件结论中的两个不等式，与之前实变量的压缩函数几乎是一个模式，就会知道巴拿赫不动点定理理解起来一点也不难。这也充分说明任何抽象理论都源自简单的具体模型。
(3)d(x, y) ≤ d(x, z) + d(z, y)。
其次，什么叫“完备”的距离空间？在微积分中用两数x和y之间的通常距离|x – y|可以定义数列的收敛性：数列{xn}收敛到数x，假如距离数列{|xn – x|}当n趋于无穷大时趋向于0。用完全同样的方式可以定义距离空间(X, d)中元素序列的收敛性：我们说X中的序列{xn}收敛到X的一个元素x，如果距离数列{d(xn, x)}当n趋于无穷大时趋向于0。类似地，实数集合中的柯西数列也可直接推广到距离空间中的序列：X中的序列{xn}被称为是柯西序列，若双指标数列{d(xm, xn)} 当m和n都趋于无穷大时趋向于0。我们知道在实数中任何柯西数列都有极限，但并非每一个距离空间都有这一性质，就像有理数组成的柯西数列不一定收敛到有理数那样。如果一个距离空间(X, d)中的每一个柯西序列都有极限，那么(X, d)称为是完备的。
既然闭区间是我们大学时代就熟悉的一个代表性完备距离空间，上面关于压缩函数的不动点定理的全部证明，只要将绝对值符号改为距离符号，完全就可以移植为巴拿赫不动点定理的证明，读者可以自行写出这一证明。巴拿赫不动点定理的一个标准应用是常微分方程初值问题解的存在唯一性定理证明，这时初值问题的解被表达成将连续函数映成连续函数的某个积分算子的不动点，这里定义在某个闭区间[a, b]上的所有连续函数全体构成一个距离空间，它的距离函数定义为d(f, g) = max{|f(x) – g(x)|: a ≤ x ≤ b}。
再回到牛顿迭代法，我们已知它总是局部收敛的，但如果初始点取得距离方程f(x) = 0的解x*太远了，那会有什么结果？牛顿法非局部的收敛性问题是一门大学问，苏联数学家康特诺维奇 (Leonid Kantorovich，1912-1988) 于上世纪中叶发展了半局部收敛理论，在James M. Ortega和Werner C. Rheinboldt的那本名著中有详细的讨论。该理论的特点是，即便初始点没能选在保证局部收敛的局部收敛域内，只要某些关键不等式被满足，就能保证迭代点的序列收敛到方程的解。自然，许多初始点却满足不了这些条件，容易构造一个例子，在许多点出发的牛顿迭代数列发散，如下图所示。
事实上，让牛顿法不收敛的那些初始点可以形成复杂无比的点集。上世纪80年代的某个学期，康奈尔大学数学教授哈伯德 (John Hubbard，1945-) 在法国一所大学度学术假，同时给一年级新生讲授一门微积分课程。有次他灵光一闪，将自己从事的一部分事业——复离散动力系统——和计算数学中古老的牛顿法联系在一起。复离散动力系统，顾名思义，就是迭代复变函数看看最终走向如何。
中学生就已经知道像 2 + 3 i这样的复数是什么，这是几百年前为了迎合求解二次方程 x2 + 1 = 0 的需要不得不创造出来的“虚无缥缈”之数。这“由实数到虚数”的一个大跳，解决了许多数学大问题，例如代数基本定理——在复数范围内，每个代数方程都有根。复动力系统和分形几何这两门当代学科紧密相连。
教授了太多次微积分课本里介绍的牛顿法，有点厌倦标准教法的哈伯德想换一换花样。他把目光转向在复数平面上用牛顿法解最简单的三次方程z3 – 1 = 0，即算出1的三个立方根，分别是实数1和两个互相共轭的复数(-1 + i√3)/2 和(-1 - i √3)/2。这三个根在复平面上形成一个等边三角形。任取一个复数作为初始点，他让学生们看看牛顿法将引向三个根中的哪一个。这个创造性的习题给学生铺下了一条发现之路！
自然，根据牛顿法的局部收敛性，当初始点取得靠近三根之一时，迭代点复数列将收敛到那个根。但是，哈伯德让学生用计算机编程找出复平面上哪些初始点将走到第一个根，哪些点趋向于第二个根，哪些点会导致第三个根。这三类初始点分别用三种不同的颜色区别开来。在粗糙的选点下，牛顿法迭代的最终性态果然如他所猜把平面分成三个扇形，但随着选点越来越精细，他和学生们发现这三个区域的分界线越来越不清楚：三种颜色互相缠绕，只要两种颜色靠近一些，第三种颜色便乘虚而入，挤进来夹在它们中间，就好像红黄蓝三股绳绞在一起彼此纠缠不清。这又引起一连串新的自相似的涌入，似乎没有哪个点可以分开任意两种颜色。就这样，美国数学教授哈伯德和修其课程的法国大学生们意外地发现了已经广泛使用了几百年的牛顿法所制造的分形图。
出品：科普中国
联系客服
