Numerical Analysis.docx
- 文档编号:28812379
- 上传时间:2023-07-19
- 格式:DOCX
- 页数:24
- 大小:122.26KB
Numerical Analysis.docx
《Numerical Analysis.docx》由会员分享,可在线阅读,更多相关《Numerical Analysis.docx(24页珍藏版)》请在冰豆网上搜索。
NumericalAnalysis
NumericalAnalysis
LloydN.Trefethen
OxfordUniversity
March2006
致谢:
这一个文章将会在不久就要出版的普林斯顿数学伴侣中刊登,由TimothyGowers和JuneBarrow-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)求nxn矩阵的特征值。
(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年代发明了计算机。
从这时刻开始,数值数学爆炸了。
但是主要地在专家手中。
涌现了很多新的杂志,如MathematicsofComputation(1943)和NumerischeMathematik(1959)。
这一革命与硬件交互映辉,但是它包括的却是与硬件没有多大关系的数学和算法。
从1950年代起的半世纪中,计算机的速度加快了大约109,但是某些问题闻名的最好运算法也加快了那么多。
两者组合后速度的增加几乎难以置信。
半世纪来,数值分析已经发展成数学中最大的分支之一,数以千计跨越多种科学和工程学科的研究人员在数十个数学杂志和应用杂志中发表了文章。
由于过去几十年的中这些人的努力,和由于强有力的计算机,我们已经达到这个水平,即大部份物理学遇到的古典数学问题能被数值方法以高的准确性解决。
使这成为可能的大部份的算法是1950以后发明的.
数值分析是建立在一个坚强的基础上的:
那就是数学中近似值理论。
这个领域包含插值的古典问题,级数展开和与牛顿,傅立叶、高斯有关的调和分析和其它;半经典多项式问题和与Chebyshev和伯恩斯坦等的名字相关的有理数的极大极小近似值问题,以及样条函数、径向基函数和小波。
我们没有篇幅来讨论这些主题,但是在几乎数值分析的每个领域中,迟早总要涉及到近似值理论。
3机器算术和舍入误差
人们都知道,计算机不能够准确地表现实数和复数。
比如1/7这个商数在一部计算机上将不能得出精确的结果。
(如果我们专门设计了基7的机器,那会是另一回事)。
计算机是用浮点算法系统来近似表示实数的,其中每个数值是用一个科学符号的等价物来表现,这样,除非数值超过了上溢值或低于下溢值,否则比例尺的大小不引起精度问题。
浮点算术由KonradZuse在1930年代在柏林发明,在1950年代底之前它成为计算机的工业标准。
在1980年代之前,不同的计算机有不同的算术特性。
在1985年,在数年的讨论之后,二进制浮点算术被采用为IEEE(电气和电子学工程师学会)标准,或IEEE短的算术标准。
这一个标准后来已用于几乎全世界许多类型的处理器上。
一个IEEE(双精度)实数包括一个64位二进制,53位二进制表示一个带符号的分数,11位表示一个带符号的指数。
因为
IEEE数可以表现实数的相对精度到16位小数。
因为
这个系统可以表示
范围内的数值。
计算机不仅仅表示数,当然;他们要对数值进行运算,例如加、减、乘、除,并从这些初等运算序列中获得更多复杂的结果。
在浮点算术中,每个初等运算计算的结果可以在下列的意义下改正:
如果“*"是这四种运算的理想形式之一,而
则是在计算机上做的相同的运算,则对于任何浮点数值x和y,假定没有上溢或下溢,可以写成:
这里ε是非常小的量,它的绝对值小于被称为machineepsilon的最小数值,写成εmach。
这个数值也度量了计算机的精度,在IEEE系统中,
。
因此,在一部计算机上,区间[1,2],被分为约1016个数值。
有趣的是比较一下这个数值离散化和物理学的离散化的精细程度,取一把固体或液体或一个气球大的气体,将它的的原子或分子排成一列,从一根线上的一端到另一端,其原子数大约为108(Avogadro数的立方根).一个这样的连续性系统足以作为衡量我们对实际量的定义,如密度、压力、应力、应变和温度等。
然而,计算机算术却比它超过了百万倍。
和物理学的另外一个比较是基本的常数的已知精度,如重力G为4位数字,蒲朗克常数h和基本电荷e为7位数字,光速c有8位数字,电子磁矩和包尔磁子之比
为12位。
目前,在物理学方面几乎找不到准确性的超过12位的数字。
因此IEEE数字是比任何科学数字更精确的量级。
(当然,纯粹的数学量π是另一回事)
因此,在两个意义上,浮点算术远比物理量更靠近它的理想值。
然而却有一种奇怪的现象,很多人竟把浮点算术而不是物理法则看作为一丑陋的和危险的折中。
数值分析家他们自己部分是为这个感觉而提出责备。
在1950和1960年代中,这个领域的先行者们发现不正确的算术可能是一个危险的来源,导致本来应该正确的结果却是错误的。
这个问题的来源是数值不稳定:
就是舍入误差在特定模态的计算中被放大,从显观的变到巨观的比例。
包括vonNeumann,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,andSylvester在十九世纪末发现的,并由Golub和其他的数值分析家在大约1965年把它做出名的。
如果A是一个
的矩阵,则A的SVD是一个因子分解
其中U是具有正交列的mxn矩阵,V是nxn具有标准正交的列的酉矩阵,而Σ是对角元素
的对角矩阵。
人们可以把SVD联系到AA*或A*A的特征值问题来求解,但这会出现数值不稳定的问题;比较好的方法是用不把A变成方阵的QR分解算法的一个变种。
当A是方阵并且非奇异时,计算SVD是求范数
以及A的逆的范数
或它们的乘积(条件数)的标准路线:
它也是进一步计算一些特殊问题中的一个步骤,包括欠秩系统的最小二乘问题,范围和零空间的计算,求秩、“全最小二乘“、低秩近似、求子空间之间的夹角等。
以上所有的讨论都涉及产生1950-75时期的“古典的”数值线性代数,此后的四分之一世纪引进一组全新的工具:
基于Krylov子空间迭代为基础的解大规模的问题方法。
这些迭代的思路如下。
假如给定一个含有大维数矩阵的线性代数问题,例如n>>1000。
它的解的特点可能是一个满足一定变化特性(例如使
最小或者)的向量
,也可能是
的一个静止点(当A为对称时Ax-λx的解)。
现在如果Kk是Rn的k维子空间k< 对Kk的奇妙选择是Krylov子空间 对于起始向量q。 因为与近似值理论有密切联系的理由,如果A的特征值分布是适当的话,随着k的增加,这些子空间的解时常非常快速地收敛到在Rn中的正确的解。 例如,时常可用来解决包括105个未知数的矩阵问题,只用数百次迭代,就可得到10位数字精度。 与古典的运算法则相比,其加速可能提高了数千倍。 Krylov子空间迭代开始于在1952年发表的共轭梯度和Lanczos迭代,但是在那时候,计算机还不够强大,不足以在大维数问题上表现出这个方法的竞争力。 1970年代,随着Reid和Paige,尤其是vanderVorst和Meijerink(他们提出了著名的预梳理概念)的工作,这个方法得到起飞。 对一个系统Ax=b进行预梳理,就是把它用一个数学上等价的系统来代替 MAx=Mb 对于一些非奇异的矩阵M.如果恰当地选择M,新的包括矩阵MA的问题可能具有比较合理分布的特征向量,因此,Krylov子空间迭代可以较快地得到它的解。 从1970年代以后,预梳理矩阵迭代已经成为计算科学不可缺少的工具。 作为其突起的一个指标,我们可以注意到,在2001年,汤姆森ISI宣布,在1990年代中引证最多的数学文章是1989年vanderVorst介绍的Bi-CGStab,那是一篇关于非对称的矩阵共轭梯度通用化的文章。 最后,我们一定要提到数值分析中最大的被未解决问题。 任意的nxn矩阵A一定能用O(nα)次运算(α>2)求逆吗? (这个问题和解方程Ax=b或计算矩阵乘积AB的问题是等价的)。 高斯消元法的α=3,1990年Coppersmith和Winograd发表的某些迭代算法(虽然不实用的)说明指数α可能缩减到2.376.是不是还有一种“快速矩阵求逆”的方法没被发现呢? 5微分方程的数值解 数学家发展解决分析的问题数值方法比线性代数要早得多。 数值积分或求面积的问题可以回溯到高斯和牛顿,甚至阿基米德。 最经典的求积公式起源于数据插值的概念,即用n次多项式对n+1个数据点进行插值,然后准确地积分这个多项式。 用等间隔插值点可得到Newton-Cotes公式,它对小的多项式是有用的,但是当n→∞时会按2n的速率发散: Runge现象。 如果这些点被最佳地选择,则结果是高斯求积,积分会快速地收敛并且数值稳定。 因为这些最佳的点是勒让德多项式的根,它们聚集在端点的附近。 对于大多数的目的,同样好的是Clenshaw-Curtis求积,其中的插值点变成cos(jπ/n),0≤j≤n。 这个求积方法也是稳定和快速收敛的,它不像高斯求积那样可按快速的傅立叶变换的方法用O(nlogn)次运算完成。 要说明为什么有效的求积必需要聚集点的规则,这就涉及到位势理论的主题。 大约在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年代以后已经属于JohnButcher。 对于精度p=6,7,8,最小的级数是s=7,9,11,当p>8时,精确的极小级值是不知道的。 幸亏这些高阶要求在实际问题中是很少遇到的。 在二次世界大战后,当开始用计算机来解微分方程的时候,再一次出现了一种最具有重要实际意义的现象: 数值不稳定。 和以前一样,这一个术语指的是局部误差被一个计算的程序无限放大。 ,但是现在占优势的局部误差通常是离散化误差而并非舍入误差。 这种不稳定典型地表现为在计算结果中一个振动的误差,当采取更多的计算步骤,这种误差以指数规律飞速增加。 关心这个现象的数学家是GermundDahlquist,他认识到这个现象可能用很强的力量和普遍性进行分析。 一些人视为他1956年的发表文章标记了现代数值分析的诞生。 这一划时代的文章介绍了可能被称为数值分析的基本定理: 一致性+稳定性=收
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Numerical Analysis