给2011年将要学习“计算方法”课程的学生

此文为西安电子科技大学前副校长 陈怀琛 老师推荐的文章,中文部分为陈老师翻译共享给我们年轻教师看的,这里也共享给我的同行和学生

-------------------------------------------------------------------------------------------------

译者说明:

哈尔滨工程大学的廖振鹏院士向我推荐了这篇文章,他在来信中说附件是我最近看到的一篇短文。我觉得写得好,对科学计算和线性代数的地位作了简洁、清楚和准确的概述,虽然我对其中个别观点有所保留。”我看了以后,很赞同廖先生的意见,觉得这篇从数学家的角度对数值分析和线性代数的综述文章对我们的项目很有价值,这也表示了廖先生对我们项目的关心与指导,我觉得有必要让参加项目的全体老师都得益,所以发给大家。为了使大多数读者能够快速浏览,我化了几天把它粗粗地翻译了一下,有些文字我也没有琢磨透,仅供参考,看不懂译文就请看原文,并请把您的译文发给我hchchen1934@vip.163.com。

从这篇文章至少可以看出以下几点:

1. 数值分析的全部发展过程和重点都是以面向应用(即面向科学计算)的需求作为主导方向的;

2.  在这些应用中,方程的阶次n都是成千上万的,低于千算是小的,大多数要考虑十万以上;高阶系统的计算速度和精度一直是工程师和数学家们最求的目标;

3.  所有的计算都是用计算机来实现的,算法的效率要跟上硬件的效率,这种不断的跟进是促进线性代数和数值分析近几十年来发展的主要动力;

4.  线性代数近几十年来的发展是围绕着计算的稳定性。从主元交换、LU分解、QR分解、特征值分解、奇异值分解、… 的不断出现,其实都是为了计算机实现的可靠的。

5.  线性代数算法的进步也不断促进了软件包的发展:EISPACK, LAPACK, LINPACK再到MATLAB, Maple和Mathematica, 直到1993年把LINPACK作为世界高速计算机排行榜TOP500的测试标准;

6.  数值分析课程必须要与计算程序的实现结合起来讲,才能辨别算法的价值,所以必须要结合一门语言,结合计算机。不能都从最基本的语法起步,要利用成熟软件的子程序,才能避开细节而突出算法。工程师和数学家算题都要用商用软件作为工具的。

7.  常微分方程和偏微分方程的解要以矩阵求解作为基础;涉及偏微分方程(组)的三维多物理场的计算和仿真已经有了商用软件,COMSOL就是很优秀的一个,你只要指明问题的物理基础(即方程性质),输入对象形状的图纸和材料,规定其三维的边界条件,根本无需知道计算的细节,就可以得出解。

8.  数值最优化的方向将导致程序的自适应化,把程序和算法相结合是未来发展的方向。

由此可知,即使按照数学家的观点,线性代数和数值分析也是随科学计算的需要而发展起来的,特别是近五十年的飞速进步都是紧跟计算机才得到的。脱离了计算机,这两门数学就失去了前进的生命力。我们把两者结合起来,不仅是工程与数学结合的需要,也是恢复了它们在数学发展上应有的历史面目。

这篇文章主要对象是数学老师,特别是线性代数和数值计算课的教师。对于学生,特别是非数学系的学生,这篇文章偏深一些,可以引用它的某些历史事实,更多从工科的需求来讲这个问题。我们很希望得到数学界和工程界权威和老师们提供这样的文章。

                                                            西安电子科技大学   陈怀琛 2011-1-25

 

----------------------------------------------------------------------------------------------

以下是翻译的部分,有希望看原文的同学请向我索取

 

Numerical Analysis

Lloyd N. Trefethen

Oxford University

March 2006

致谢:

这一个文章将会在不久就要出版的普林斯顿数学伴侣中刊登, 由Timothy Gowers和June Barrow-Green编辑,普林斯顿大学出版社出版。

在准备这一篇短文时,我得益于来自许多同事的忠告,他们也帮助改正了若干事实和论述的错误。我也不总是听从他们的忠告,因此我将对书中的错误和省略负完全责任。(以下是致谢的名单,略去)

 

1 对数值计算的需要

每个人知道当科学家和工程师需要数学问题的数值答案的时候,他们就转向计算机。然而,存在着很大的关于这一个过程的误解。

数值力量是特别巨大的。 人们时常认为,当伽利略等定下了一项原则,对每件事物一定要进行测量,从那时起科技革命就开始起动了。数值测量导致了用数学方法表达物理定律。并且,在全部周期中,其成果都围绕着我们,精细的测量导致精确的定律, 进而导致较好的技术和更精确的测量。 那种离开数值数学而获得某种物理科学的进步,或获得某种意义重大的工程产品发展的时代,早已过去了。

在这一个故事中计算机确实扮演一个角色,不过,它们的角色是什么却存在一个误会。许多人想像,由科学家和数学家产生公式,然后,藉着将数值插入进这些公式之内,计算机就制造出必需的结果。实际情况完全不是这样。 真的进行的是执行运算法则的一个更为有趣程序。在大部份的情形下,照着公式做这件工作甚至无法完成,大多数的数学问题不能够靠一个有限步操作的序列来解决。相反的是快速算法则很快地收敛于精密到三或十位,甚至一百位的数值近似答案。对于科学或工程应用来说, 这样一个答案可能和精确答案一样有用。

可以举例说明正确和近似解的复杂性。假如我们有一个四次多项式,

p(z)= c0+ c1z+ c2z2+ c3z3+c4z4;

而另外有一个五次多项式,

q(z)= d0+ d1z+ d2z2+ d3z3+d4z4+d5z5:

广为人知的是:p的根可由显式(由Ferrari在大约 1540 年发现)求得, 但是q的根却没有这样的公式(Ruffini和Abel在250年之后; 证明了它的无解性)。因此,哲学上会意义到,p 和 q 的求根问题是完全不同的。然而在实际应用中它们却难于区别。如果一位科学家或一个数学家想要知道一个多项式的根, 他将会转向一部计算机而且在小于几毫秒的时间内得到16位数值精度的答案。 计算机使用了一个显式吗?在q的情况,答案肯定为不,但p怎么样?也许是,也许不。大部份的时间,使用者既不知道也不关心,一百个数学家中也许找不到一个能凭记忆写下求p的根的公式。

这里再举出另外三个例子,就像对于p的求根那样,它们是能用一级数初等运算求解的。

(1) 线性方程组: 解含n个未知数的n个一次方程组。

(2) 线性规划: 在m 个线性约束下,将含n个变量的一次函数减到最小。

(3) 旅行售货员问题: 找到在 n 城市之间的最短旅游路线。

而下面的五个例题, 则像对于q求根那样,通常是不能够用初等运算求解的。

(4) 求 n x n 矩阵的特征值。

(5) 求多变量函数的最小值。

(6) 计算积分。

(7) 解常微分方程(ODE问题) 。

(8) 解偏微分方程(PDE) 。

我们能否得出结论 (1)-(3) 在实际中将会比(4)-(8)容易 ? 完全不是!如果 n 是数百或数千,问题 (3) 通常是非常难解的。 问题 (6)和(7) 通常相当容易,至少如果积分是一维的。 问题 (1)和(4) 几乎完全有相同的难度:当n很小的时候(例如 100)比较容易,而当 n大的时候,(例如1000000)就很难。事实上, 在这些问题中,哲学只能对实践作很差的指导, 在问题 (1)-(3)中, 当n和m很大的时候,人们一般不去求精确的解,而使用近似的(但却是快速的!)解法。

数值分析是研究连续问题运算法则的数学, 这意味着命题包含实变量或复变量。(这个定义包括在实数域定义的线性规划和旅行售货员那样的问题, 但不是它们的离散对应物.) 在本文的其他部分,我们将综述它的主要分支、已有成就和可能的未来趋向。

2 简短的历史

在整个世纪中,领先的数学家已经参与了科学应用, 而且在许多情况下,这已经导致今天的仍然在使用的算法的发现。高斯就是一个杰出的例子。在许多其他的贡献中, 他在最小二乘数据拟合(1795)、线性方程组求解(1809)、和数值积分(1814)方面,都推动了决定性的进步。他在发明快速傅立叶变换 (1805)方面也一样, 虽然后者直到1965年Cooley和Tukey 把它再发现后,才变成广为人知。

大约在1900年左右,在数学家的研究活动中,数值分析开始变得不大活跃。因为技术上的理由,当时数学的进步主要集中在严格性的问题上。举例来说, 二十世纪初期数学家的许多结果要用新的关于无穷大的严格方法来论证,这些命题和数值计算相去甚远。

一代人过去了,在1940年代发明了计算机。从这时刻开始,数值数学爆炸了。但是主要地在专家手中。涌现了很多新的杂志,如Mathematics of Computation (1943)和Numerische Mathematik (1959)。这一革命与硬件交互映辉,但是它包括的却是与硬件没有多大关系的数学和算法。从1950年代起的半世纪中,计算机的速度加快了大约 109,但是某些问题闻名的最好运算法也加快了那么多。两者组合后速度的增加几乎难以置信。

半世纪来,数值分析已经发展成数学中最大的分支之一,数以千计跨越多种科学和工程学科的研究人员在数十个数学杂志和应用杂志中发表了文章。由于过去几十年的中这些人的努力,和由于强有力的计算机, 我们已经达到这个水平,即大部份物理学遇到的古典数学问题能被数值方法以高的准确性解决。使这成为可能的大部份的算法是1950 以后发明的.

数值分析是建立在一个坚强的基础上的: 那就是数学中近似值理论。这个领域包含插值的古典问题,级数展开和与牛顿,傅立叶、高斯有关的调和分析和其它;半经典多项式问题和与 Chebyshev 和伯恩斯坦等的名字相关的有理数的极大极小近似值问题,以及样条函数、径向基函数和小波。 我们没有篇幅来讨论这些主题,但是在几乎数值分析的每个领域中,迟早总要涉及到近似值理论。

3 机器算术和舍入误差

人们都知道,计算机不能够准确地表现实数和复数。比如1/7这个商数在一部计算机上将不能得出精确的结果。(如果我们专门设计了基7的机器,那会是另一回事)。计算机是用浮点算法系统来近似表示实数的,其中每个数值是用一个科学符号的等价物来表现,这样,除非数值超过了上溢值或低于下溢值,否则比例尺的大小不引起精度问题。浮点算术由 Konrad Zuse 在1930年代在柏林发明,在 1950 年代底之前它成为计算机的工业标准。

在1980年代之前,不同的计算机有不同的算术特性。在 1985 年,在数年的讨论之后, 二进制浮点算术被采用为 IEEE(电气和电子学工程师学会) 标准,或IEEE短的算术标准。这一个标准后来已用于几乎全世界许多类型的处理器上。一个IEEE(双精度)实数包括一个64位二进制,53位二进制表示一个带符号的分数,11 位表示一个带符号的指数。因为 ,IEEE 数可以表现实数的相对精度到16位小数。因为 ,这个系统可以表示 范围内的数值。

计算机不仅仅表示数,当然; 他们要对数值进行运算,例如加、减、乘、除,并从这些初等运算序列中获得更多复杂的结果。

在浮点算术中,每个初等运算计算的结果可以在下列的意义下改正: 如果“*"是这四种运算的理想形式之一,而 则是在计算机上做的相同的运算,则对于任何浮点数值 x 和 y, 假定没有上溢或下溢,可以写成:

 

这里ε是非常小的量,它的绝对值小于被称为machine epsilon的最小数值,写成εmach。这个数值也度量了计算机的精度,在IEEE系统中, 。

因此,在一部计算机上,区间[1,2],被分为约1016个数值。有趣的是比较一下这个数值离散化和物理学的离散化的精细程度,取一把固体或液体或一个气球大的气体,将它的的原子或分子排成一列,从一根线上的一端到另一端,其原子数大约为 108(Avogadro数的立方根). 一个这样的连续性系统足以作为衡量我们对实际量的定义, 如密度、压力、应力、应变和温度等。然而,计算机算术却比它超过了百万倍。和物理学的另外一个比较是基本的常数的已知精度, 如重力G为4位数字, 蒲朗克常数h和基本电荷e为7位数字,光速 c有8位数字,电子磁矩和包尔磁子之比 为12位。目前,在物理学方面几乎找不到准确性的超过 12位的数字。因此 IEEE 数字是比任何科学数字更精确的量级。(当然,纯粹的数学量π是另一回事)

因此,在两个意义上,浮点算术远比物理量更靠近它的理想值。然而却有一种奇怪的现象,很多人竟把浮点算术而不是物理法则看作为一丑陋的和危险的折中。数值分析家他们自己部分是为这个感觉而提出责备。在1950和1960年代中,这个领域的先行者们发现不正确的算术可能是一个危险的来源, 导致本来应该正确的结果却是错误的。这个问题的来源是数值不稳定:就是舍入误差在特定模态的计算中被放大,从显观的变到巨观的比例。包括von Neumann, Wilkinson, Forsythe, 和 Henrici这些人,费尽心力宣传在盲目相信机器算术会导致的危险。这些危险确实非常真实,但是信息传播的太成功了, 导致了广泛的印象:以为数值分析的主要工作是控制舍入误差。 事实上, 数值分析的主要工作是设计收敛得很快的运算法则;舍入误差的分析是时常讨论的一部份,但很少成为中心议题。如果舍入误差消失,90% 的数值分析会保留下来。

4 数值线性代数

线性代数变成大学本科数学课程的标准的主题开始于1950和 1960 年代, 而且保持至今。这有多个理由,但是我认为最根本的一个理由是: 线性代数的重要性是从计算机抵达以后发生爆炸的。

这个主题以高斯消元法为出发点,这是一个能在n个一次方程组中,使用用n3数量级算术运算次数,解出n个未知数的程序。相等地,它可解决形如Ax=b的等式,其中A是一个nxn的矩阵而x和b都是n维的列矢量。在全球的计算机上每解一次线性方程组,就要调用一次高斯消元法。即使 n大到1000, 在一个典型的 2006年制造的台式机上需要的时间不到一秒。

消元法的思想最初产生于约2000年以前的中国学者,较近的贡献者包括拉格朗日,高斯和 Jacobi。直到1944年Dwyer才用现代方法叙述这种算法。假如,用α乘以第一行并从第二行被减去。 这一个运算可以解释为把A阵的左边乘一个下三角矩阵M1,它包含一个单位矩阵以及另外一个非零元素m21 =-α.进一步类似的行运算对应于进一步左乘下三角矩阵 Mj。如果经过k步后,A转换成上三角矩阵矩阵U,则有MA= U, 其中M= Mk…M2M1.如果,设L=M-1,则有

A =LU

这里L是单位下三角矩阵,也就是说,所有对角元素均为一的下三角矩阵. 因为 U代表目标结构,而L则代表实现U的运算编码,我们可以说高斯消元法是一个分解为下三角乘上三角的过程。

许多其他数值线性代数的运算法则也都是基于将一个矩阵写成几个有特殊特性的矩阵的乘积。从生物学中借一个词汇,我们可以说在这个领域有一则中心的教条:

运算法则↔矩阵分解:

按这一个结构,我们能很快地描述出需要的下一个运算法则。不是每个矩阵都有LU 因子分解的; 举一个2x2的矩阵作为反例:

 

在应用计算机之后很快就观察观察到,即使确实有 LU 因子分解的矩阵,它对于纯粹形式的高斯消元仍可以是不稳定的, 其舍入误差被可能扩大得很厉害。在消元期间交换行,以便把最大的元素到对角线上(主元交换法)就可以实现稳定。因为主元交换是行运算,它又可用左乘一个其他矩阵来实现。具有主元交换功能的高斯消元法的分解公式就是

PA=LU;

其中U 是上三角矩阵, L 是单位下角矩阵,而P是交换矩阵,也就是产生行交换的单位矩阵。如果选择的交换是把列 k 的对角线下最大的元素放到(k,k)的位置上(在第k步消元之前),那L将具有附加的特性,对所有i,j,元素

主元交换法发现得很早,但是它的理论上的分析证明却令人惊异地难。 在实践中,主元交换法几乎使高斯消元法达到了完全的稳定,并且成为几乎所有的电脑解线性方程组时都采用的程序。它的证明是在大约 1960 年由威尔金森和其他人完成的,不过对于某些特定的例外矩阵,高斯消元法(甚至主元交换后)仍然是不稳定的。对这个差别的缺乏解释表现了在数值分析的核心还存在令人困窘的缝隙。实验建议,由高斯消元法放大的舍入误差的倍数大于 的因子应该在一定意义下指数式地随 减小。其中n为维数(例如,具有独立正态分布函数作为元素的随机矩阵),但是这个定理从未被证明过。

在1950年代后期开始, 数值线性代数在另外一个方向得到发展: 利用基于正交或酉矩阵为基础的算法。实酉矩阵满足 ,复酉矩阵满足 ,其中 指的是共轭转置。这一发展的起点是QR分解的概念。如果A是一个 的矩阵,A的QR分解就是如下的乘积:

A=QR

其中Q由正交的列组成,而R是上三角矩阵。 可以把这个公式解释为Schmidt正交化的概念,其中Q的各个列向量q1,q2,…是逐个依次产生的。这些列运算对应于在A的右边乘以单位上三角矩阵。可以说Gram-Schmidt算法的目的是求Q,副产品是R,所以它是一个三角正交化过程。当Householder在1958年给出了具有正交三角化的双重策略的算法对许多目的更为有效的时候,这成了一件大事。在他的方法中,藉由连续应用多个初等矩阵运算,每个初等矩阵运算反映一个超平面Rm, 人们通过正交运算把矩阵A简化为上三角矩阵:目标是R,而Q成了副产品。Householder在数值方面更为稳定,因为正交运算能保持范数不变,不会扩大在每个步骤中引入的舍入误差。

在 1960 年代中,从QR分解引出了许多线性代数运算法则。QR 因子分解能独自解决最小二乘问题以及构造标准正交基。更显著的是把它作为一个步骤使用在其他的算法中。尤其,数值线性代数的中心问题之一是求方阵A的特征值和特征向量,如果A有一组完整的特征向量,然后用这些特征向量列构成矩阵X,用对应的特征值构成对角矩阵D,我们获得

AX=XD;

因为 X 是非奇异的, ,称特征值分解。特别当A是 hermitian时, 总是存在一组完全的标准正交特征向量,给出

 

其中Q是酉矩阵。 计算这些因子分解的标准算法是在1960年代早期由Francis, Kublanovskaya 和Wilkinson开发的。因为不存在五次以上多项式的求根公式,所以特征值一般地不能用闭合形式求出,因此QR算法必然地是迭代的,它应该包含一系列的(原则上是无穷的)QR分解。 然而,它的收敛特别迅速。在对称矩阵中,对于一个典型的矩阵A,QR算法按立方规则收敛, 其意义是每一次迭代后,特征值和特征向量的精确数字差不多增加三位。

QR 算法是数值分析的伟大成就之一, 它对广泛使用过的软件产品产生的冲击是巨大的。这算法和以它为基础的分析在 1960 年代中引到计算机的Algol和FORTTRAN中,稍后又用到软件库EISPACK(特征系统软件包" 及其后裔LAPACK中。这个方法也已经被吸收在通用的数值计算库如NAG,IMSL 和数值程序集中,也应用在如 MATLAB,MAPLE 和 Mathematica 的求解环境中。这些发展已经是如此的成功以致于矩阵特征值的计算已经事实上变成每位科学家的运算‘黑箱工具’,没有几个专家真正知道它的实现细节。有一个有趣的相关故事, EISPACK的亲戚,用来解线性方程组的LINPACK有了一个料想不到的功能: 它成为计算机制造业者测试他们的计算机的速度标准的出发点。如果一个超级计算机能幸运地进入前500名(TOP500)目录(该名单自从 1993 以后一年更新两次),那是因为它在解决维数从100到数百万的特定 矩阵方程Ax=b方面是的表现超凡。

所有的数学家都熟悉特征值分解, 但是数值线性代数也已经把它的较年轻的堂兄弟带到舞台上:那就是奇异值分解(SVD).SVD是由Beltrami, Jordan, and Sylvester 在十九世纪末发现的,并由Golub 和其他的数值分析家在大约 1965 年把它做出名的。

如果A是一个 的矩阵,则A的SVD是一个因子分解

 

其中U是具有正交列的mxn矩阵,V是 nxn具有标准正交的列的酉矩阵, 而Σ是对角元素 的对角矩阵。人们可以把SVD联系到AA*或A*A的特征值问题来求解,但这会出现数值不稳定的问题;比较好的方法是用不把A变成方阵的QR分解算法的一个变种。当A是方阵并且非奇异时,计算SVD是求范数 以及A的逆的范数 或它们的乘积(条件数)的标准路线:

 

它也是进一步计算一些特殊问题中的一个步骤,包括欠秩系统的最小二乘问题,范围和零空间的计算,求秩、“全最小二乘“、低秩近似、求子空间之间的夹角等。

以上所有的讨论都涉及产生1950-75时期的“古典的”数值线性代数, 此后的四分之一世纪引进一组全新的工具: 基于Krylov子空间迭代为基础的解大规模的问题方法。这些迭代的思路如下。假如给定一个含有大维数矩阵的线性代数问题,例如n>>1000。它的解的特点可能是一个满足一定变化特性(例如使 最小或者)的向量 ,也可能是 的一个静止点(当A为对称时Axx的解)。现在如果Kk 是Rn的k维子空间k<<n,它有可能在那个子空间中加快解决相同的变化问题。对Kk的奇妙选择是 Krylov 子空间

 

对于起始向量q。 因为与近似值理论有密切联系的理由, 如果A的特征值分布是适当的话,随着k的增加,这些子空间的解时常非常快速地收敛到在Rn中的正确的解。例如, 时常可用来解决包括 105个未知数的矩阵问题,只用数百次迭代,就可得到10位数字精度。与古典的运算法则相比,其加速可能提高了数千倍。

Krylov 子空间迭代开始于在1952年发表的共轭梯度和Lanczos迭代, 但是在那时候,计算机还不够强大,不足以在大维数问题上表现出这个方法的竞争力。1970年代,随着Reid和Paige,尤其是van der Vorst和Meijerink(他们提出了著名的预梳理概念)的工作,这个方法得到起飞。对一个系统Ax=b进行预梳理,就是把它用一个数学上等价的系统来代替

MAx=Mb

对于一些非奇异的矩阵M.如果恰当地选择M, 新的包括矩阵MA的问题可能具有比较合理分布的特征向量,因此,Krylov 子空间迭代可以较快地得到它的解。

从 1970 年代以后, 预梳理矩阵迭代已经成为计算科学不可缺少的工具。作为其突起的一个指标,我们可以注意到,在 2001年,汤姆森 ISI 宣布,在1990 年代中引证最多的数学文章是 1989年 van der Vorst 介绍的Bi-CGStab, 那是一篇关于非对称的矩阵共轭梯度通用化的文章。

最后,我们一定要提到数值分析中最大的被未解决问题。任意的nxn 矩阵A一定能用O(nα)次运算(α>2)求逆吗?(这个问题和解方程Ax= b或计算矩阵乘积AB的问题是等价的)。 高斯消元法的α=3, 1990年Coppersmith和 Winograd发表的某些迭代算法(虽然不实用的)说明指数α可能缩减到2.376.是不是还有一种“快速矩阵求逆”的方法没被发现呢?

5 微分方程的数值解

数学家发展解决分析的问题数值方法比线性代数要早得多。 数值积分或求面积的问题可以回溯到高斯和牛顿, 甚至阿基米德。

最经典的求积公式起源于数据插值的概念,即用n次多项式对 n+1个数据点进行插值, 然后准确地积分这个多项式。用等间隔插值点可得到Newton-Cotes公式,它对小的多项式是有用的,但是当n→∞时会按 2 n 的速率发散:Runge 现象。如果这些点被最佳地选择,则结果是高斯求积, 积分会快速地收敛并且数值稳定。因为这些最佳的点是勒让德多项式的根,它们聚集在端点的附近。对于大多数的目的,同样好的是Clenshaw-Curtis求积,其中的插值点变成 cos(jπ/n),0≤j≤n。这个求积方法也是稳定和快速收敛的,它不像高斯求积那样可按快速的傅立叶变换的方法用O(nlog n)次运算完成。要说明为什么有效的求积必需要聚集点的规则,这就涉及到位势理论的主题。

大约在1850年另一个解析问题开始得到注意: ODE(常微分方程)问题的解。Adams公式是以等间隔点插值为基础的,在实践中一般少于10点。这是现在称为多步法的ODE问题的数值解法。这里的概念是对于自变量为t的初值问题 ,我们取一个小的时间步长△t>0 并且考虑一个时间值的有限集

 

然后用一个我们能够计算的代数近似值的问题代替ODE,计算出一序列的近似值

 

(此处上标只代表上标,不是指数)。这类最简单的近似公式(用Euler法)是

 

或者使用缩写

 

ODE问题本身和它的数值近似值两者都可能包括一个或多个方程,在多个方程时, u(t,x) 和 vn 成为适当维数的向量。Adams公式是 Euler's 的公式在高阶条件下的一般化,它能更为有效地产生正确的解。 举例来说,四阶Adams-Bashforth 公式是

 

术语“四阶”反映了处理数值问题的一种新的元素,说明在△t→0时问题收敛的外貌. 上述公式为四阶,在某种意义上表明它将会按O((△t)4)的速率收敛。在实践中应用的阶数多半在3~6之间,它能使各种计算都得到卓越的精度,通常有3-10位有效数字。当需要更多准确性时,偶或也用更高阶的公式。

最不幸地是,数值分析文献习惯上所说的不是这些优秀高效方法的收敛性, 而是它们的误差,更准确地说,它们是与舍入误差不同的离散化误差或截断误差。到处存在的误差分析的语言,语调阴郁, 但是似乎无法根除。

在二十世纪的转折点,出现了第二类优秀的ODE问题算法,人们熟知的Runge- Kutta法或一步积分法被 Runge,Heun和Kutta发明了。举例来说,下面是著名的四阶Runge-Kutta法公式,它借助于对函数f的四次估值,把一个数值解(标量或向量)从时刻tn推进到tn+1。

 

Runge-Kutta法比较容易实现,但是分析起来有时比多步公式难。举例来说,对于任何的s,推导s步Adams-Bashforth公式的系数是一个很小的事,它应该具有p=s的精度等级。但相比之下,Runge-Kutta法在级的数目(每一步需估值的函数数)和所能达到的精度等级之间没有简单的关系。s=1,2,3,4的经典方法早已由Kutta在 1901年得到为p=s,但是直到 1963 年才证明了s=6级所能达成的精度为p=5。分析这个问题涉及图论和其他领域的美丽数学,其中一个关键领域从 1960 年代以后已经属于John Butcher。 对于精度p=6,7,8,最小的级数是 s=7, 9,11,当p>8时,精确的极小级值是不知道的。幸亏这些高阶要求在实际问题中是很少遇到的。

在二次世界大战后,当开始用计算机来解微分方程的时候,再一次出现了一种最具有重要实际意义的现象:数值不稳定。和以前一样,这一个术语指的是局部误差被一个计算的程序无限放大。,但是现在占优势的局部误差通常是离散化误差而并非舍入误差。这种不稳定典型地表现为在计算结果中一个振动的误差,当采取更多的计算步骤,这种误差以指数规律飞速增加。关心这个现象的数学家是Germund Dahlquist,他认识到这个现象可能用很强的力量和普遍性进行分析。一些人视为他1956年的发表文章标记了现代数值分析的诞生。

这一划时代的文章介绍了可能被称为数值分析的基本定理:

一致性 + 稳定性 = 收敛性:

这个理论是基于对这三个概念的精确定义。 一致性是离散化公式的局部性误差满足正的精度等级因而正确地为ODE问题建了模。稳定是指的在某个时间步中引入的误差不会在以后无限增长。收敛性是指当△t→0,并且没有舍入误差时, 数值解收敛于正确的结果。在Dahlquist的报告之前, 关于稳定性和收敛性的等价性的认识几乎是空的,也许在某种意义上,实际工作者知道,如果一个数值方案不是不稳定的, 它或许会得到一个接近正确的答案的近似值。他的理论把这个概念以严密的形式用到了很广阔的各类数值方法中。

当ODE问题的计算机方法正在发展时, 同样的研究对于更大的PDE(偏微分方程)的主题也正在进行。大约1910年,Richardson发明了应用于应力分析和气象学的解PDE的离散数值方法,并更进一步由 Southwell发展;1928年,由Courant,Friedrich 和 Lewy发表了一篇关于有限差分法的理论文章。虽然 Courant-Friedrich- Lewy 的工作稍后变得出名,但在计算机出现前的这些概念所产生的影响向前被约束有限的。此后这个主题发展得很快。 在早期特别有影响的是在洛萨拉摩斯实验室的von Neumann周围的研究者群体, 包括年轻的Peter Lax。

如同ODE问题一样, Neumann 和他的同事发现,一些解PDE的方法也受制于悲惨的数值不稳定。举例来说,要用数值方法解决线性波动方程 ut= ux,我们可能按通常的网格提取空间和时间步长△x 和△t

  

并用代数公式代替PDE计算一序列的近似值:

 

众所周知,为这一个目的所做的离散化是Lax-Wendroff公式:

 

其中 它也能被推广到一维服从双曲线保守律的非线性系统。对于 ut=ux, 如果λ被固定在低于或等于1,当△x 和△t→0时(不考虑舍入误差),这个方法将会收敛到正确的解,另一方面,如果λ大于1,它就会爆炸。 Von Neumann 和其它人认识到,可以测试出来有没有这种不稳定,至少对线性常系数问题,可以用对x 的离散傅立叶分析:“Von Neumann 分析”来解决。经验指出,在实际的情况下,一个方法如果不是不稳定的,它就会成功。而且很快出现一个给这观察以严格证明的理论:Lax等价定理是由Lax和 Richtmyer在1956年发表的,和Dahlquist的报告在同一年。许多细节是不同的—这理论被约束为线性方程组然而为ODE问题的 Dahlquist理论也适用于非线性方程—但是广义地说,新的结果符合使收敛性等于一致性加稳定性的相同模式。数学上,关键点是一致和无界的原理。

自从Von Neumann逝世半世纪以来,Lax- Wendroff公式和它的亲戚已经发展成一个慑人心魄有力的主题—计算流体动力学。早期的在一维空间处理的线性和非线性方程式很快转移到二维,最后发展到三维空间。现在在三个方向各有几百个网格点,总计有几百万个变量的问题已经是常规的工作了。方程组可以是线性的或非线性的; 网格是均匀的或不均匀的, 通常可以自适应地加密以便对边界层和其他变更剧烈的区域以特别的注意;几乎应用在各个地方。 数值方法首先被用来为螺旋桨建模,然后是机翼, 最后是整个的飞机。工程师设计时仍然使用风洞,但是他们更多地信赖计算的结果。

这些成就中的大部份是被1960年代内出现的另一项解PDE的数值技术推进的,那就是有限元法,它来自工程学和数学的不同根基:它不是用一个差商来近似一个微分算子,有限元法用能够瓦解为简单小块的函数f来逼近的是解本身。举例来说,人们可以把f的定义域划分初等形状小块的集, 如三角形或四边形,而且坚持f对每块的约束是低阶的多项式。在对应的有限维子空间解形式变化了的PDE来获得所需的解,而且时常可以发现算得的解在那个子空间里面是最佳的。有限元方法已经利用函数分析的工具发展得非常成熟了。这些方法在运算复杂的几何学方面以灵活性闻名, 特别是他们的在结构力学和土木工程中占有完全的优势。关于有限元法的已经被出版的书和文章数已经超过了10000。

在巨大的和成熟PDE的数值解的的领域中,现在状态中的什么方面会使Richardson 或 Courant,Friedrich和Lewy 吃惊呢? 我认为那是全世界的对线性代数的运算法则的依赖。立体的大规模的 PDE 问题的解可能需要在每个时间步长中解一百万个方程的系统。利用有限差分预梳理器,和依赖于另一个多网格预梳理器的双CSGtab迭代,再用GMRES(广义最小余量法)矩阵迭代,这可能被达成。把工具作如此堆积,当然是早期的计算机开拓者们无法想像的。对它的需要最后归结到数值不稳定, 如Crank和Nicilson在 1947年首先注意的那样,与不稳定奋斗的决定性工具是使用隐含的公式, 它耦合了新的时间步骤 tn+1 处的未知数, 而在实现这一个耦合时,又要用到方程组的解。

这里是一些例子,说明今天的科学和工程学对PDE数值解的成功与信赖: 化学,(Schrodinger 方程式),结构力学 (弹性方程式); 天气预报(地球风方程式);涡轮设计 (navier-Stokes方程式);声学 (Helmholtz 方程式);无线通信(Maxwell方程式); 宇宙学 (爱因斯坦方程式); 油矿探寻(迁移方程式); 地下水再调度 (Darcy定律); 集成电路设计 (杂质扩散方程式); 海啸模型 (局部浅水方程式);光纤;(非线性波动方程式);图像加强 (Perona-Malik方程式);冶金(Cahn-Hilliard方程式);期权定价问题(Black-Scholes 方程式).

6 数值最优化

数值分析的第三个重要分支是最优化,那就是把多变量函数的值减到最小限度,与其相接近的问题是求非线性方程组的解。最优化的发展已经略微独立于数值分析的其余领域,部分是由一些与运筹学和经济学接近的学者群体推进的。

由微积分我们知道一个平滑的函数可能在导数为零的点或在边界点达到极值。这两种可能性也表示最优化领域的两大方面。在一端是求不受约束非线性的导数为零内点和极小值方法,涉及多变量微积分问题。在另一端是线性规划的问题,其中被减到最少的函数是线性的,因此容易的懂得,所有的挑战是边界约束条件。

不受约束的非线性最优化是一个旧的主题。 牛顿引入了我们现在称之为泰勒级数的前几项来逼近函数的概念; 的确,Arnol'd甚至认为泰勒级数是牛顿的“主要数学发现。人人都知道,要求出实变量x 的函数 F 的零点 x*,牛顿法的思路是:在 第k步,给出一个估计x(k)=x*, 用导数F0(x(k))来定义一个线性函数从而得到一个较好的估计 x(k+1):

 

牛顿(1669)和Raphson(1690)把这一个思路应用到多项式,而辛普森(1740)把它推广到其他的函数 F 和两个方程的系统。用今天的语言, 对于含有n未知数的n阶方程组,我们把 F视为一个n维矢量,它在点 x(k)2 Rn 的导数是 nxn Jacobian 矩阵,它的元素为:

 

这矩阵定义了F(x)的一个线性的近似值,但在x≈x(k)点是精确值。牛顿法然后取它的矩阵形式

 

这在实践中意味着若想从x(k)求得x(k+1),我们需要解一个线性方程组:

 

只要 J 是 Lipschitz 连续和非奇异的,而且对 x的开始猜测足够好,这一迭代的收敛是二次的:

      

学生时常认为把公式的指数提高到3或4可能是一个好思路,然而,这是一个幻想。每次做两步二次收敛的算法产生一个四次收敛算法——因此在二次和四次算法之间在效率方面的差别,充其量是一个常数因子。如果指数 2,3,或 4被任何其他大于1的数值代替,结果是一样的.在所有这些收敛算法之间真实的区别是超线性的,, 牛顿法是其原型, 那些按线性或几何规则收敛的,其指数是只有1.

从多变量微积分学的观点,从解一个方程组以使一个标量函数f 2 Rn 减到最少只是一小步,想得到一个(局部)最小量,我们寻求梯度 g(x)的零点,它是一个n维矢量。g的导数是Jacobian 矩阵,称为f 的 Hessian,其元素为:

 

人们可在一个牛顿迭代之前利用它来找到g(x)的一个零点, Hessian 具有的新特征是它总是对称的。

虽然对于最小化和求零的牛顿公式是早就形成的,计算机的到来也创建了一个新的数值最优化的领域。它很快地遇到的障碍之一是如果开始的猜测不好,牛顿法时常会失败。对这一个问题已经从实际上和理论上进行了广泛的研究,采用的算法技术是线搜索和可信赖区域。

对于多变量的问题,人们很快地发现在每个步骤评估 Jacobians 或 Hessians可能是过度的。需要利用不太精确但却较快速的方法算 Jacobians 或 Hessians,来解出伴随的线性方程组的不精确解, 最后仍然达成高线性的收敛。这一类型的最早突破性的发展是1960年由Broyden,Davidon,Flecher和Powell发现的准牛顿法,其中部分的数据用来产生稳定地改良对真实的 Jacobian 或 Hessian 或它的矩阵因素的估计。能说明这个主题在那时紧急性的一个例证事实是在1970年,秩二对称正定准牛顿更新公式被独立地被不少于四位作家同时发表, 即 Broyden,Flecher, Goldfarb 和 Shanno;由此他们的发现已经从此被称为著名的BFGS公式。在后来的数年中,处理问题的规模指数地增加,新的思路变成很重要, 包括自动微分,一种能使被计算的函数是导数被自动地决定的技术: 计算机程序本身是“可求导的",以便在产生数值输出的同时,它也产生了它们的导数。“自动求导数”是一个旧的梦,但是由于种种理由,部分是由于稀疏线性代数的进步, 直到1990年代Bischof, Carle和 Griewank的工作出现之前,它没有能够实现。

无约束的最优化问题相对比较容易,但是它们不是典型的; 已经被发展的能处理约束的方法才显示出这个领域的真正深度。假如函数 f: Rn! R 是被减到最小的函数,它受到等式约束 cj(x)=0 和不等式约束 dj(x)_0,其中 fcjg 和 fdjg 也是从 Rn 到 R的函数. 对这样的问题求解办法即使在局部最优性情况也是不简单的,它包括拉格朗日乘子和在主动的和非主动的约束之间建立区别。 这一个问题现在被所谓的 KKT 条件解决了。它是在 1951 年被Kuhn和Tucker以及, 后来知道, 它是12 年之前由Karush引入的。有约束的非线性最优化算法的开发今天继续是一个活跃的研究主题。

有约束问题把我们带到了另一个数值最优化领域—线性规划。 这一个主题在 1930 年代和 1940 年代由在苏联的Dantzig和美国的Kantorovich诞生。如一个自然的发展

作为在战争中对于的美国空军工作的副产品,Dantzig在1947年发明了出名的解决线性的程序算法—单纯形法。线性规划是将一个受到 m个线性的等式及[或]不等式约束的n个变数的线性函数减到最小的问题。这如何能称作挑战? 答案是 m 和 n 可能是很大的。大规模的问题可能因连续问题离散化而出现,也可能由他们自身产生。一个出名的早例子是Leontiev's 的理论经济学的投入产出模型, 在1973年使他嬴得了诺贝尔奖。甚至在 1970 年代中,苏联用了包括数以千计变量的计算机输入-输出模型作为经济计划的一个工具。

单纯形算法使得中规模和大规模的线性规划问题易于计算。这样一个问题是由它的目标函数定义的,被减到最小的函数f(x),和它的可实行区域,即所有能满足约束条件的矢量集 x 2 使 Rn。

对于一个线性规划,可实行的区域是一个多面体,一个由超平面包围的闭域,而且 f 的最佳值必定在其一个顶点达到。(如果一个点是定义约束方程组的某个子集的唯一解,就称为一个顶点)。单纯形算法由有系统地从一个顶点到另外的顶点作下山移动直到达到最佳的点为止。全部迭代结构都位于可实行区域的边界上。。

在 1984 年, 在 AT&T 贝尔实验室的 Narendra Karmarkar 在这个领域引起了一个剧变。Karmarkar 证明了人们有时凭借在可实行的区域内部工作,可以做得比单纯形算法好得多。 Karmarker's 的方法很快的向前发展了, 尤其把一个问题的初始和双重串接在一起进行研究的决定性思路。现在被称为初始-双重内点法。 这些算法非常有力,而且他们有不强烈地靠线性的优良特征。因此,除了变更线性规划领域之外,Karmarkar 和他的继承者还带来了把最优化的线性和非线性两个方面健康地相互靠近在一起的效果。今天, 这领域得到了很大的进步,它处理的线性或非线性问题的规模已可包括数百万个变量和约束。

7 未来

数值分析是从数学中跳出来的; 然后它繁衍到计算机科学。1960年代,当大学开始设立计算机科学系的时候,数值分析家时常是领引者。经过两代人之后, 现在他们的大部分却到了数学系。发生什么事? 部份的答案是数值分析家喜欢处理连续的数学问题,而计算机科学家偏爱不连续的问题,这个缝隙显得真是很宽。

然而,数值分析的计算机科学这边是非常重要的,而且我想要以一个强调这一个方面的预测来结束这篇文章。传统地人们可能把一个数值算法当做老套的程序——运行一些类型的循环直到满足一个预定的结束标准为止。 对于某些计算这个现象是正确的。 另一方面,1960 年代由de Boor,Lyness和 Rice的工作领头,开始出现了比较不确定型的数值计算: 自适应算法。 在一个最简单类型的自适应求积程序中,在特定的网格的不同部分上同时计算着两个积分,然后进行比较,得到对局部误差的一个估计。基于这一个估计, 可以把部分的网格改得较密以改善准确性。 这一程序被迭代运行,直到答案符合使用者预先规定的目标容差为止。

大多数的这样计算不能保证结果的准确性,但是令人兴奋的新发展是更复杂的后验误差控制技术,它有时确实提供保证。当和间隔算法技术相结合时,甚至舍入误差和离散化误差的期望准确性也可得到保证。

首先,电脑程序对于求积分变成自适应的; 然后是对ODE问题的程序。对于PDE,向自适应程序的推进是在一个比较长的时段上发生的。最近已经在傅立叶转换,最优化, 和大规模的数值线性代数,以及一些新的适应计算机构造和数学问题的算法有了相应的发展。 在一些著名的算法可以解决每个问题的世界中, 我们逐渐发现,最强健的电脑程序将是一个具有各种能力并能自适应地对这些能力进行安排调配的程序。换句话说,数值计算逐渐地正在埋入在智能控制环路。我相信,这一个过程将会继续, 正如已经发生在许多其他技术领域的那样,把科学家从他们的计算细节中远移,代之以稳定增长的计算能力。我预计2050年电脑的数值程序中99% 将会是只能的“程序包",只有1%是真正的“算法”,如果这样区别还有意义的话。没有人会知道程序是如何工作的,但是这些程序将会特别有力和可靠,并时常给出有保证的准确结果。

这个故事将会有一个数学推论。在数学方面的基本的区别之一在于线性问题能用一步来解决,而非线性问题则通常需要迭代。相关的区别是在前向问题(一步) 和倒转的问题 (迭代)之间。当数值算法正在逐渐地被埋入在智能控制环中,几乎每个问题将会被迭代处理,不管它的哲学状态。代数学的问题将会被分析的方法解决; 而且在线性的和非线性之间,前向的和倒转的问题之间,区别将会模糊起来。

附录: 一些主要的数值算法

表 1 的清单尝试确认一些在数值分析的历史发展上最有意义的算法(当做理论的对立面)发展。在每个情形下,都或多或少按年代引证了一些早期的关键面貌并给出了关键的最早年份,当然,任何这的历史简述一定是过分简化了的。省略名字的痛苦在清单各处都发生——包括超过一半的 EISPACK,LINPACK和LAPACK库的作者。甚至日期也是有问题的:举例来说,快速的傅立叶变换列在1965年,因为那是使它在世界闻名的文章发表的年份,虽然高斯在160年之前有了相同的发现。同样地,线性的规划的内点方法是按 Karmarkar的突破性发展列为1984年,虽然Khachiyan的重要相关的工作早五年就出现了。也不能想像从1991到现在的若干年是空白! 无疑未来我们将找到这时期值得放在这个表中适当位置的成就。

posted @ 2013-04-18 19:22  dutlei  阅读(698)  评论(0编辑  收藏  举报