计算方法绪论与误差等

计算方法

NumericalAnalysis能源与动力工程学院刘火星82316418liuhuoxing@课程介绍

绪论

误差分析

线性方程组的解法

常微分方程的初值问题矩阵的特征值和特征向量

插值和拟合

数值微分和数值积分

非线性方程解法CourseOutline徐翠薇、孙绳武，计算方法引论（第二版），高等教育出版社，2002RichardL.Burden&J.DouglasFaires,NumericalAnalysis(SeventhEdition),高等教育出版社，2001数值分析，李庆祥等编，高等教育出版社，2000主要参考书成绩评定方法平时成绩：20%期末考试:80%先修课程高等数学线性代数程序设计语言绪论Introduction关于计算方法计算数学的过去和未来算法和收敛数学软件目次

计算机解决实际问题的步骤建立数学模型选择数值方法编写程序上机计算实际问题数学问题提供计算方法程序设计上机计算结果分析实际问题工程问题假设实验模型数学模型

校验分析/数值解测量物理定律修改

在计算机上是否根据数学公式编程就能得到正确结果?研究例子:求解线性方程组其准确解为x1=x2=x3=1如把方程组的系数舍入成两位有效数字它的解为x1=-6.222...x2=38.25…x3=-33.65...Numericalanalysisisthestudyofalgorithmsfortheproblemsofcontinuousmathematics----LloydN.Trefethen“计算方法”就是研究在计算机上解决数学问题的理论和数值方法今天的数值计算方法，无论从形式到内容，还是从工具到效果，已远非半世纪前VonNeumann、Lax等先驱们所处的环境和条件了，计算机技术和应用软件的发展，让计算数学展开了双翼。许多迅速发展的其他学科和社会进步给计算数学的发展开拓出更为广阔的新天地计算方法研究的对象研究数值方法的设计、分析和有关理论基础与软件实现。计算方法又称：计算数学、数值方法、数值分析等计算方法的分枝有最优化方法、计算几何、计算概率统计等计算方法的内容

连续系统的离散化

离散性方程的数值求解

串行计算方法与并行计算方法的关系：并行计算方法70年代初随并行计算机的出现而产生，是计算数学中最活跃的新领域。计算数学的发展与科学工程计算是紧密相联的，计算数学的发展历史也就是与其他学科结合，利用计算机不断形成新的理论及数值方法并不断形成新的学科的历史，例如：“计算物理”“计算流体力学”计算数学发展的历史回顾H.Aiken(1900-1973)哈佛大学博士，因做博士论文涉及到空间电荷传导问题的计算，1937年提出方案，1939年得到IBM资助，1944年建成投入使用。这是继电式计算机－MarkI三位计算机设计大师的贡献J.W.Mauchly(1907-1980)宾夕法尼亚物理博士，因从事天气预报需要想设计计算机，1942年提出计算机方案，1945年底竣工，这就是世界上第一台电子计算机－ENIAC机三位计算机设计大师的贡献J.VonNeumann(1903-1957)普林斯顿高级研究所，1945年在普林斯顿研制成MANIAC机，有力地支持美国氢弹研制，称为计算机之父三位计算机设计大师的贡献在研制原子弹和氢弹过程中，许多物

理规律必须通过计算机上的计算摸清

楚。计算物理、理论物理与实验物理

相辅相成相互促进共同发展，形成现

代物理学的三大分支由于核武器研制需要，1950年全球只有15台，到了1962年9月仅美国就有16187台计算机计算数学发展的历史回顾1983年一个由美国著名数学家拉克斯(P.Lax)为首的不同学科的专家委员会向美国政府提出的报告之中，强调“科学计算是关系到国家安全、经济发展和科技进步的关键性环节，是事关国家命脉的大事。”1984年美国政府大幅度地增加对科学计算经费的支持,新建成五个国家级超级计算中心（分别在普林斯顿大学、圣地亚哥、伊里诺大学、康奈尔大学、匹兹堡），配备当时最高性能的计算机，建立NSF-net新网络

80年代中期我国将“大规模科学与工程计算”列入国家资助重大项目1987年起美国NSF把“科学与工程计算”、“生物工程”“全局性科学”作为三大优先资助的领域由于大存储的高速计算机的使用已导致了科学和技术方面的两大突出进展：大量用于设计工作的实验被数学模型的研究逐步取代，如航天飞机设计、反应堆设计、人工心瓣膜设计等能获取和存储大量的数据，并能提取隐秘的信息，如计算机层析X射线摄影，核磁共振等1991年以美国总统倡议的形式提出了“高性能计算与通信计划”。这是为了保持和提高美国在计算和网络的所有先进领域中的领导地位而制定的。计划为期五年（1992－1996），投资的重点是发展先进的软件技术与并行算法，关键技术是可扩展的大规模并行计算要求到1996年高性能计算能力提高14倍，达到每秒万亿次浮点运算速度（1012

Teraops/S）。计算机网络通迅能力提高1百倍，达到每秒109位（Gigabits/S）该计划中列举的“挑战”项目有：磁记录技术、药物设计、催化、燃烧、海洋模拟、臭氧洞、空气污染、高速民用运输机、数字解剖、蛋白质结构设计、金星成像等1993年初美国总统发布“发展信息高速公路”（NII）的总统令1994年4月美国总统发布“建立国家（地球）空间数据基础实施”（NSDI）的总统令数值方法和数值软件过去50年的主要进展Before1940Newton’smethod;Gaussianelimination;Gaussquadrature;leastsquaresfitting;AdamsandRunge-Kuttaformulas;Richardsonextrapolation1940-1970floatingpointarithmetic;Fortran;finitedifferences;finiteelements;FFT;simplexalgorithm;MonteCarlo;orthogonallinearalgebra;splinefunction1970-2000quasi-Newtoniterations;adaptivity;stiffODEsolvers;softwarelibraries;Matlab;multigrid;sparseanditerativelinearalgebra;spectralmethods;interiorpointmethods计算数学未来50年的展望将更多的通过声音，而不是键盘向计算机传递信息，而计算机将更多地以图象而不是数字反映结果数值计算将更具有适应性、迭代性、灵活性。计算能力大得惊人数值计算中更具智能性算法：一系列近似计算步骤的组成，目的是找到问题的近似解算法的特征：收敛性、稳定性计算量：一个算法所需的乘除运算总次数，单位是flop.计算量是衡量一个算法好坏的重要标准算法和收敛矩阵乘积AB的计算量分析a11a12a13…a1na21a22

a23…a2n...

...

…...am1am2

amm-1…amnb11b12b13…b1sb21b22

b23…b2s...

...

…...bn1bn2

bnn-1…

bns=[cij]ms因为cij=aik

bkj

计算量为n所以上面AmnBns的计算量为N=mn

s误差分析ErrorAnalysis误差的来源误差误差限有效数字相对误差和绝对误差误差的传播在近似计算中需要注意的问题目次模型误差观测误差舍入误差截断误差1.1

误差的来源

计算机字长有限，一般实数不能精确存储，于是产生舍入误差。例如：在10位十进制数限制下：1÷3＝0.3333333333本应1÷3＝0.3333333333……1.0000022-1.000004=0本应1.0000022-1.000004

=1.000004000004-1.000004

=0.0000000000041.1.1

舍入误差(Round-offErrors)舍入误差很小，本课程将研究它在运算过程中是否能有效控制。1.1.2截断误差(TruncationError)用近似的值去代替数学上的准确值带来的误差。例如：

泰勒级数•零阶近似:•一阶近似:•二阶近似:完全的泰勒级数:余项(n阶近似)::介于

xiand

xi+1

x

=

xi+1-

xi

余项：Taylor级数表示为:截去的部分•零阶近似:

截断误差

:•

一阶近似

Rn:零阶近似Rn:斜率:1.2

误差误差限有效数字[Def1.1]若用x*表示x准确值的一个近似值。则此近似值x*和准确值x的差称为误差，用e*来表示

e*＝x*－x[Def1.2]若

|e*|＝|x*－x|≤ε*ε*称为近似值x*的误差限。[例1.2]已知x*=π=3.14159…，求近似值x1=3.14，x2=3.142，x3=3.1416的误差限。[解]所以误差限

ε1=0.002，ε2=0.0005，ε3=0.000008有效数字[Def1.3]若用x的近似值x*的误差限是某一位上的半个单位，该位到x*的第一位非零数字共有n位，则称x*有n位有效数字若用x*表示x的近似值，并将x*表示成

x*＝±0.a1a2…an×10m若

|x*－x|≤0.5×10m－n则近似值x*有n位有效数字(1.1)[例1.3]

设x*=0.0270是某数x经“四舍五入”所得，则误差|e(x*)|不超过x*末位的半个单位，即：

|x*－x|≤0.5×10-4

又x*=0.27×10-1,故该不等式又可写为

|x*－x|≤0.5×10-1-3由有效数字定义可知,x*有3位有效数字,分别是2,7，0。[例1.4]

设x＝32.93,x*＝32.89,则

|x*－x|＝0.04＜0.05＝0.5×10-1即

|x*－x|≤0.5×102-3由有效数字定义可知,x*有3位有效数字,分别是3,2,8。由于x*中的数字9不是有效数字，故x*不是有效数。1.3相对误差和绝对误差设x——准确值x*——近似值称为近似值x*的相对误差实用中，常用表示近似值x*的相对误差，称为相对误差限相应的，e*称为绝对误差，ε称为绝对误差限有效数位与误差的关系有效数位n越多，则绝对误差|e*|越小形如(1.1)式的近似数x*具有n位有效数字，则其相对误差限可取为基本算术运算设x*和y*分别是x和y的近似值，把它们的误差近似地看做是相应地微分，即

dx≈x*－x，dy≈y*－y则

d（x±y）＝dx±

dy d（xy）＝xdy±

ydx d（x/y）＝（－xdy＋ydx

）/y21.4

误差传播(1.3)和(1.4)给出了由自变量的误差引起的函数值的误差的近似式（误差传播）。一元函数设y＝f(x)，若x的近似值是x*，用f(x*)去近似f(x)的误差可用Taylor公式估计

(1.3)(1.4)多元函数情形由多元函数的Taylor展开公式类似可得

(1.5)

(1.6)

(1.8)

(1.9)[例1.5]测得某桌面的长a的近似值a*=120cm,宽b的近似值b*=60cm。若已知|e(a*)|≤0.2cm,|e(b*)|≤0.1cm。试求近似面积s*=a*b*

的绝对误差限与相对误差限。解:面积s=ab,在公式（1.5）中,将y=f(x1,x2)

换为s=ab,则相对误差限为1.5

在近似计算中需要注意的问题1.尽量简化计算步骤，减少乘除运算的次数

例如，计算多项式通常运算的乘法次数为若采用递推算法,

则乘法次数仅为n.又如2.防止大数“吃掉”小数当|a|>>|b|时，尽量避免a+b

。例如，假设计算机只能存放10位尾数的十进制数，则

108+0.04=1083.尽量避免相近数相减例如，当x很大时，应当x接近于0时，应4.避免绝对值很小的数做分母当|b|<<|a|时，应尽量避免a/b

5.

选用数值稳定性好的算法，以控制舍入误差高速增长第2章线性代数方程组解法本章研究的对象是n

阶线性代数方程组2.1对象(2.1)线性系统广泛存在于工程、科学以及社会科学、商业和经济问题的定量分析等领域中用矩阵和向量的记法来表示，(2.1)式可写成对象（续）(2.2)其中A＝(aij)是方程组(2.1)的系数aij构成的n×n阶矩阵，称为系数矩阵。B＝{bi},X＝{xi}是n维向量，X是未知量，B称为右端项。使方程组(2.1)中每一个方程都成立的一组数x1*,x2*,…,xn*称为式，(2.1)的解，把它记为向量的形式，称为解向量。

我们总是希望方程组有解,且有唯一解.由线性代数的克莱姆(cramer)规则可知,如果方程组(2.1)的系数矩阵A的行列式不等于零,那么,这个方程组有唯一解,而且它们可以表示为

xi=Di/D(i=1,…,n)这里,Di是指D中第i列元素用右端b1,…bn代替构成的行列式.

如果方程组(2.1)有唯一解,我们按上面的等式求解,就必须计算n+1个n阶行列式.由行列式的定义,n阶行列式包含有n!项,每一项含有n个因子,计算一个n阶行列式就需要做(n-1)n!次乘法.而我们一共要计算n+1个n阶行列式,共需做(n2-1)n!次乘法.此外,还要做n次除法才能算出xi(i=1,…n).也就是说,用这个办法求解就要做

N=(n2-1)n!+n次乘除法运算,这个计算量是大得惊人的.例如,当n=10(即求解一个含10个未知量的方程组),乘除法的运算次数共为32659210次;

消(元)去法是求解线性方程组(2.2)和满秩矩阵的逆阵A-1的一种直接方法.尽管它比较古老,但它具有演算步骤,推算公式都系统化的特点(对其中主元素消去法,还可以证明是稳定的).因此,它至今仍然是求解方程组的一种有效的方法.

消去法可以引出几种计算方法,下面按三角形方程组和一般线性方程组的顺序来讨论。消去法上三角方程组的一般形式是：对于(2.1)式，有以下运算可简化方程组线性方程组的解法有两类：直接法：即在没有舍入误差的情况下，用有限步的四则运算得出精确解的方法。但实际运算中舍入误差不可避免，此类方法也只能得到近似解。目前常用的是列主元消去法和矩阵三角分解法迭代法：先给一个初始值，按一定法则逐步求解出各个更准确的近似值的方法。目前常用的有Jacobi法、Seidel迭代法、松弛法和梯度法将A和B写在一起，称为增广矩阵将例2.1用增广矩阵的变化写成高斯消去法的求解过程分为两个阶段:首先，把原方程组化为上三角形方程组，称之为“消元”过程；然后，用逆次序逐一求出三角方程组(原方程组的等价方程组)的解，并称之为“回代”过程。2.2高斯消去法(Gaussianelimination)消元：将(2.1)式写成矩阵形式(2.3)(2.4)第1步:若a11(1)

≠0，用第二个方程减去第一个方程乘以a21(1)/

a11(1),用第三个方程减去第一个方程乘以a31(1)/

a11(1)…则有矩阵形式(2.5)(2.6)第2步:若a22(2)

≠0，用第三个方程减去第二个方程乘以a32(2)/

a22(2),用第四个方程减去第二个方程乘以a42(2)/

a22(2)…则有矩阵形式(2.7)(2.8)第k步:若akk

(k)

≠0，用第k+1个方程减去第k个方程乘以ak+1k(k)/

akk(k)…则有矩阵形式(2.9)(2.10)重复n－1次，得到等价的上三角形方程组矩阵形式(2.11)(2.12)以上过程把系数矩阵A(1)变成上三角矩阵A(n)，称之为消元，计算公式可归纳为(2.13)回代(2.12)2.3主元素消去法因此，x1＝0x2＝1因此，x1＝1x2＝1精确解，x1＝10000/9999x2＝9998/9999在做除法运算时，选取绝对值大的作分母。

——主元素消去法的基本思路。列主元素消去法列主元素消去法基本思想1.用高斯消去法求解线性方程组时,应避免小的主元.在实际计算中,进行第k步消去前,应该在第k列元素aik

(i=k,…,n)中找出绝对值最大者,例如

∣a∣=max∣a∣2.再把第p个方程与第k个方程组进行交换,使apk(k-1)成为主元.我们称这个过程为选主元素.由于只在第k列元素中选主元素,通常也称为按列选主元素(或称部分选主元).3.如果在第k步消去前,在第k个方程到第n个方程所有的xk到xn的系数a(i=k,…,n;j=k,…,n)中,找出绝对值最大者,例如 ∣a∣=max∣a∣再交换第k,p两个方程和第k,q两个未知量的次序,使a成为主元素.称这个过程为完全选主元素。4.不论是哪种方式选出主元素,而后再按上面介绍的计算步骤进行消去的计算,一般都称为主元素高斯消去法。在实际计算中,常用按列选主元素的高斯消去法。(k-1)(k-1)pk(k-1)ikk≦i≦n(k-1)pq(k-1)ijk≦i,j≦n(k-1)ij(k-1)pq高斯消去法的乘除总运算分析如下：消元次数k消元乘法次数消元除法次数回代乘除法总次数

1n(n-1)n-12(n-1)(n-2)n-2...k(n-k+1)(n-k)n-k..n-12*11n(n+1)/2

故高斯消去法的计算量为N=n(n2-1)/3+n(n-1)/2+n(n+1)/2=n3/3+n2-n/3

当N充分大时为n3/32.4高斯消去法的计算量分析

方法的特点

由具体计算结果可知,全主元素法的精度优于主元素法,这是由于全主元素是在全体系数中选主元,故它对控制舍入误差十分有效.但全主元素法在计算过程中,需同时作行与列的互换,因而程序比较复杂,计算时间较长.列主元素法的精度虽稍低于全主元素法,但其计算简单工作量大为减少,且计算经验与理论分析均表明,它与全主元素法同样具有良好的数值稳定性,列主元素法是求解中小型稠密性方程组的最好方法之一。

选主元消去法(包括解线性方程组的所有直接的方法)比较适用于中小型方程组.对高阶方程组,即使系数矩阵是稀疏的,但在计算中很难保持稀疏性,因而有存储量大，程序复杂等不足,所幸的是这一缺点可用迭代法解决.

另外,高斯选主元消去法还可技巧性的解决一些特殊线性方程组。由于误差的影响,计算过程中可能出现一些坏现象,要尽可能防止,表现在求解线性方程组的消元法比较上,则应该注意要简化运算,减小运算次数,提高效率;提高数值稳定性等.

在计算过程中，由于计算机字长的有限性，不可避免地产生舍入误差。同时，由于所求问题的初始数据（例如线形方程组的系数矩阵和右端项系数）往往是带有一定误差的。因此计算结果总是不可避免地带有误差，或者说，如果初始数据有扰动，势必将带来具有一定误差的计算结果。就拿Ax=b来说，由于观测或计算等原因，线性方程组两端的系数A和b都带有误差A和b，这样实际建立的方程组是近似方程组(A+A)(x+x)=b+b。对近似方程组求出的解是原问题的真解x加上误差x,即x+x。而x是由A及b引起的，它的大小将直接影响所求解的可靠性。这种解依赖于方程组系数的误差A及b的问题，称为线性方程组解对系数的敏感性。2.5线性方程组解对系数的敏感性方程组

此方程组的准确解为x1=0,x2=-1。现将其右端加以微小的扰动使之变为：

经计算可得准确解为x1=2,x2=-3.

这两个方程组的解相差很大，说明方程组的解对常数项b的扰动很敏感。

病态方程组：如果方程组AX=b由于A或b的