研究生数值分析课件
河海大学研究生数值分析课件
若 P(x) 是次数不超过n的多项式，即
P( x) a0 a1 x an x n
则称 P(x)为插值多项式。相应的方法称为多项式插值。 若 P(x) 是分段多项式，则称分段多项式插值。 常用的有拉格朗日插值、牛顿插值、埃尔米特插 值、埃特金插值、三次样条插值等。
定义2 称
f ( x1 ) f ( x0 ) f [ x0 , x1 ] x1 x0
为 f (x)关于点
x0 , x1 的一阶均差；称
f [ x0 , x2 ] f [ x0 , x1 ] f [ x0 , x1 , x2 ] x2 x1
为 f (x)的二阶均差；一般的，称
f ( f ( x ,, x )) | ( ) | ( xk ) xk k 1
1 n n
例3 测量得某场地长 l 的值为 110 0.2 ，宽d m 的值为 80 0.1m ，试求面积 s = ld 的绝对误差限与 相对误差限。 （见黑板）
1.3 误差定性分析与避免误差危害
1 ( n1)
若 x 具有n位有效数字，则相对误差限
r
x 10 (a1 a2 10 an 10
) , a1 0
1 ( n 1) 10 2a1 1 ( n 1) 10 ，则 反之，若 x 的相对误差限 r 2a1
至少具有n位有效数字。 （证明见黑板）
其中数值计算方法是数值分析研究的对象。
主要包括：
（1）函数的数值逼近（包括插值法）；
（2）数值微分和数值积分；
（3）非线性方程（组）数值解； （4）数值线性代数（如线性方程组数值解、矩阵 特征值特征向量的计算）； （5）（偏）微分方程数值解。
西安科技大学研究生数值分析课件7章矩阵特征值与特征向量计算
7 矩阵特征值与特征向量地计算设A 为n 阶方阵,所谓A 地特征值问题是求数λ和非零向量x ,使x Ax λ=成立.数λ称作A 地一个特征值,非零向量x 称作与特征值λ对应地特征向量.求给定方阵地特征值与特征向量是先求解特征方程()||0E A ϕλλ=-=然后对应于每一个特征值i λ,再求解退化地齐次线性方程组()0i E A x λ-=从而得到A 地特征值i λ及对应地特征向量x .但是这种方法计算机很大,计算过程复杂,因此有必要研究相对简单地数值解法.本章主要介绍三类计算特征值地方法：计算大型（稀疏）矩阵主特征地幂法与反幂法,计算中小型（实对称）矩阵全部特征值地Jacobi 法,计算中小型矩阵全部特征值地QR 法.7.1 特征值估计在矩阵特征值计算中,有时需要对特征值所在范围给出一个估计.这里介绍一种从矩阵地元素出发,运用较简便地运算估计特征值地方法.定义7-1 设()n m ij A a C ⨯=∈,称由不等式||ii i z a R -≤在复平面上确定地区域为矩阵A 地第i 个盖尔圆（Gerschgorin 圆）,并用i G 表示.其中1||ni ij j j i R a =≠=∑称为盖尔圆i G 地半径(1,2,,)i n =.定理7-1 矩阵()n m ij A a C ⨯=∈地一切特征值均落在它地n 个盖尔圆地并集中,即1(1,2,,)ni jj G i n λ=∈=.证明 设λ是A 地任一特征值,12(,,,)T n x x x x =是λ对应地特征向量.令01||max ||i i i nx x ≤≤=,则00i x ≠.由Ax x λ=,可得001()ni j j i j a x x λ==∑.即∑≠==-ni j j j j i i i i x a x a 000001)(λ于是有 000000011i i jni j j ji ni j j i jji i i R x x ax x aa ≤≤=-∑∑≠=≠=λ这表明任一特征值0i G λ∈,从而也在A 地第n 个盖尔圆地并集中.例7-1 估计矩阵10.10.20.30.530.10.210.310.50.20.30.14A ⎡⎤⎢⎥⎢⎥=⎢⎥-⎢⎥---⎣⎦地特征值范围. 解 A 地4个盖尔圆为：1:|1|0.6G z -≤ 2:|3|0.8G z -≤ 3:|1| 1.8G z +≤ 4:|4|0.6G z +≤画在复平面上其区域如图7-1所示.图7-1 例7-1盖尔圆分布图于是A 地全部特征值就在这4个盖尔圆地并集中.为了更确切地知道某个特征值落在哪个或哪几个盖尔圆地并集中,给出如下第二盖尔圆盘定理.定理7-2 若A 地n 个盖尔圆中,有m 个盖尔圆构成地一个连通域（所谓连通域,是指其中地任意两点都可以用位于该区域内地一条折线连接起来）,且该连通域与其余n m -个盖尔圆严格分离,则在该连通域中恰好有A 地m 个特征值（重特征值按重数重复计算）.特别地,每个孤立地盖尔圆恰有A 地一个特征值（证明从略）.由定理2可知,在例1中2G 与4G 中各有A 地一个特征值,而1G 与3G 构成地连通部分中有两个特征值,但不能确定这两个特征值具体落在哪个盖尔圆中.例7-2 估计矩阵10.80.50A -⎡⎤=⎢⎥⎣⎦地特征值范围. 解 A 地两个盖尔圆为：1:|1|0.8G z -≤,2:|0|0.5G z -≤在复平面上地区域如图7-2所示.图7-2 例7-2盖尔圆分布图此时只能判断A 地两个特征值落在1G 与2G 地并集中,至于是每个盖尔圆中各有一个特征值还是两个特征值都落在其中一个盖尔圆上则无法确定.实际上,由于1,21(12λ=±,1,2||0.5λ=>,所以两个特征值都不会在盖尔圆2G 中,而是落在盖尔圆1G 中.对于某些矩阵,可利用相似变换矩阵具有相同特征值地性质得到更确切地特征值范围.设()ij n m A a ⨯=,取正数12,,,n d d d 构成对角阵12diag(,,,)n D d d d =,对A 作相似变换,令1()iij n n jd B DAD a d -⨯==,由于B 相似于A ,所以B 与A 地特征值完全相同,又由于B 与A 地主对角线元素对应相等,所以B 与A 地盖尔圆圆心相同.这表明,若适当选取正数12,,,n d d d ,可以改变盖尔圆地半径,从而有可能将相交地盖尔圆分离得到仅含一个特征值地孤立盖尔圆.选取12,,,n d d d 地一般方法是：欲使A 地第i 个盖尔圆i G 地半径大而其余盖尔圆变小,就取1i d >,其余1()j d j i =≠.例7-3 求矩阵2050.841011210A j ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地特征值范围. 解 A 地3个盖尔圆为：1:|20| 5.8G z -≤,2:|10|5G z -≤,3:|10|3G z j -≤其中1G 与2G 相交,而3G 孤立.记3G 中所含地一个特征值为3λ,如图7-3所示.为分离2G 与1G ,可以让A 地第3行元素绝对值变大,第3列元素绝对值变小.现取diag(1,1,2)D =,则12050.44100.52410B DAD j -⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦图7-3 例3盖尔圆分布图 图7-4 例7-3分离后盖尔圆分布图其3个盖尔圆分别是：1:|20| 5.4G z '-≤,2:|10| 4.5G z '-≤,3:|10|6G z j '-≤ 显然,B 地盖尔圆是3个孤立地盖尔圆,如图7-4,注意,此情况下,3G '地半径变大了.例7-4 设矩阵()ij n n A a ⨯=按行严格对角占优,则A 可逆.证明 由线性代数知,A 可逆地充分条件是||0A ≠,而1||nj j A λ==∏（其中j λ是A 地特征值）,所以只要证明0j λ≠即可(1,2,,)j n =. 设λ是A 地任一特征值,则必存在某个盖尔圆i G 使∑≠=≤-ij ij i ii a R a λ.若0j λ=,则有∑≠≤ij ij ii a a ,而这与A 按行严格对角占优矛盾,故应有0λ≠,由λ地任意性,得||0A ≠.7.2 幂法与反幂法在线性代数中,设A 是n 阶方阵,若A 存在n 个线性无关地特征向量,则称这n 个特征向量构成A 地一个完全地特征向量组.例如,对矩阵320230005A -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦,110430102B -⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦通过求解特征方程,不难求出A 地三个特征值为1231,5λλλ===,B地三个特征值为1232,1λλλ===.方阵A 可以找到三个线性无关地特征向量,而方阵B 找不到三个线性无关地特征向量.我们称方阵A 可对角化,而B 不可对角化. 7.2.1 幂法幂法地基本思想是构造一个向量序列使之逼近主特征值（矩阵地按模最大地特征值）对应地特征向量,然后求出主特征值.该方法简单易行,但收敛速度较慢.现设()ij n n A a ⨯=有一个完全地特征向量组12,,,n x x x ,其对应地特征值是12,,,n λλλ.已知A 地主特征值是单根1λ,即特征值满足条件12||||||n λλλ>≥≥任取一个非零初始向量0u ,由矩阵A 构造向量序列102210110k k k u Au u Au A u u Au A u++=⎧⎪==⎪⎪⎨⎪==⎪⎪⎩由于A 地完全特征向量组可以作为向量空间n R 地一组基,因此0u 可由12,,,n x x x 线性表示,即有01122n n u a x a x a x =+++ （设10a ≠）于是011122211111121()()k k k k k n n nn kk k i i i k i u A u a x a x a x a x a x a x λλλλλλελ===+++⎡⎤=+=+⎢⎥⎣⎦∑ 其中21()nk i k i i i a x λελ==∑.注意到),,2(11n i i=<λλ,故当k →∞时,0k ε→,因此有111k k u a x λ≈由于1x 是主特征值1λ对应地特征向量,其乘上常数因子11k a λ仍为1λ地特征向量,故当k 充分大时,迭代向量k u 是1λ地特征向量地近似向量.为了利用迭代向量求出主特征值1λ地近似值,设()k i u 表示k u 地第i 个分量,则1111111()()()[]()()()k i i k ik i i k iu a x u a x ελε+++=+ 于是 11()lim()k ik k iu u λ+→∞= 这表明两相邻迭代向量对应分量地比值收敛于主特征值,亦即当k 充分大时,可用两相邻迭代向量地分量比作为主特征值地近似值,即11()()k ik iu u λ+≈若主特征值是A 地r 重实特征值,即12(1)r r n λλλ===≤≤,对应地r 个线性无关特征向量为12,,,n x x x .则有01111()r nkk k i k i i i i i i r u A u a x a x λλλ==+⎡⎤==+⎢⎥⎣⎦∑∑当k 充分大时,11rkk i i i u a x λ=≈∑即k u 仍为主特征值对应地特征向量地近似向量,相邻两迭代向量地分量比仍为主特征值地近似值.综上所述,有定理7-3 设A 是n 阶实矩阵,具有完全地特征向量组,主特征值是r 重根,即112||||||||(1)r r n r n λλλλ++>≥≥≥≤≤则对任意非零初始向量0u ,迭代向量0k k u A u = 满足 111lim(0)rki ikk i u a x a λ→∞==≠∑ ,11()lim ()k ik k iu u λ+→∞= 或 11rk k i i i u a x λ=≈∑,11()()k ik iu u λ+≈ 这样用非零初始向量0u 及矩阵A 构造向量序列{}k u 以计算A 地主特征值1λ及相应地特征向量地方法称为幂法.不过从上面地讨论中可以看到,如果1||1λ>或11<λ,迭代向量k u 当k →∞时,其不为零地分量就会趋于无穷大或趋于零.为克服这个缺点,可以在每步迭代中加上对向量规范化地步骤,使迭代向量地数量级保持在一个稳定地量级上,归纳起来,幂法地计算步骤为：步骤 1 给定非零初始向量0u ,精度12,εε,令00v u =；令(0)10max()v λ=,1=k ；步骤 2 迭代1-=k k Av u ,()1max()k k u λ=,其中)max(k u 表示k u 绝对值最大地分量；步骤3 规范化max()kk k u v u =； 步骤 4 若11k k v v ε--<且()(1)112||k k λλε--<,则k v 即为1λ地近似特征向量,()1k λ即为1λ地近似值；否则,1+=k k ,转步骤2继续迭代.例7-5 用幂法计算1.0 1.00.51.0 1.00.250.50.252.0A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地主特征值和相应地特征向量,结果见表7-1.表7-1而此题地准确值为1 2.5365258λ= 1(0.748221,0.649661,1.000000)T x =7.2.2 幂法地加速幂法地收敛速度由比值21r λλ=来确定,r 越小收敛越快,而当1r ≈时收敛可能很慢.为了克服这一缺点,常采用原点平移法对幂法进行加速.设B A pE =-,其中p 是待定参数.显然,若A 地特征值为12,,,n λλλ,则B 地相应特征值(1,2,,)i k i n =为12,,,n p p p λλλ---,且A .B 地特征向量相同.这是因为对A 有特征方程||0i A E λ-=,而对B 有特征方程|||()|0i i B k E A p k E -=-+=,所以,i i i i p k k p λλ=+=-另一方面,若i x 是A 地对应i λ地特征向量,即i i i Ax x λ=则 ()()i i i i i i Bx A pE x Ax px p x λ=-=-=-原点平移法地思想是引入矩阵B ,恰当地选择参数p ,使11k p λ=-是B 地主特征值,且其速比2211maxB A p r r p λλλλ-=<=-,这样用幂法求B 地主特征值1k 地收敛速度就快于用幂法求A 地主特征值1λ,而一旦1k 求出,由11k p λ+=可得A 地主特征值,达到了加速地目地.但是为了选取恰当地选择参数p ,需要对A 地特征值地分布地大致了解. 例7-6 设4阶方阵A 有特征值15(1,2,3,4)j jj λ=-=其速比210.9A r λλ=≈.作变换 (12)B A pEp =-=则B 地特征值为12k =,21k =,30k =,41k =-,其速比2112B A k r r k ==<. 设A 地实特征值满足121n n λλλλ->≥≥>若2,n λλ地值可大致估计出,若要求1λ,考察待定参数p 地选取. 在原点平移法通过变换pE A B -=后,不论p 如何选取,矩阵地B 主特征值也只能是在n p λ-或 1p λ-.若希望求1λ,就应选择p ,使1p λ-称为B 地主特征值,即1||||n p p λλ->-这时B 地收敛速比B r 是比值21||/||p p λλ--和1||/||n p p λλ--中地较大者,即211||||max ,||||n B p p r p p λλλλ⎧⎫--=⎨⎬--⎩⎭显然B r 依赖于p 地选取,记做()B r p .为了使应用幂法求B 地主特征值地收敛速度尽可能快,我们希望选择最佳参数*p ,使*()min ()B B r p r p =由B r 地表示式（求二者之间地较大值）和)(*p r B 对)(p r B 地最小化要求,只有当2||||n p p λλ-=-时,()B r p 达到最小.由于2n p p λλ-=-会有得到矛盾地结果（2n λλ=）,所以只能是2()n p p λλ-=--即 *22np λλ+=类似地,若用反幂法求最小特征值n λ,若1n λ-,1λ可大致估计,取最佳平移参数*112n p λλ-+=例7-7 取0.75p =,用原点平移法,计算例7-7中矩阵A 地主特征值.解 作变换B A pE =-,则0.2510.510.250.250.50.25 1.25B ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦对B 应用幂法,计算结果见表7-2.即1 1.7865914k ≈,则A 地主特征值1λ为110.75 2.5365914k λ=+=与例7-5比较,上述结果比例7-5迭代15次还好.表7-27.2.3 反幂法设方阵A 按模最小地特征值是n λ,且0n λ≠,则A 可逆.于是,由n n n Ax x λ=,可得11n n nA x x λ-=,这表明1nλ是1A -地主特征值.反幂法就是将幂法应用于1A -,通过求出1A -地主特征值得到A 地按模最小地特征值及其对应地特征向量.定理7-4 设A 是n 阶实矩阵,具有完全地特征向量组,其特征值满足12||||||0n λλλ≥≥≥>则对任意非零初始向量00u v =,按下述方式构造地迭代向量11k k u A v --= ,max()kk k u v u =满足lim max()n k k n x v x →∞=, 1lim max()k k nu λ→∞= /max()k n n v x x ≈,1max()k nu λ≈在实际计算中,可先对A 进行LU 分解,通过求解1k k Ly v -= ,k k Uu y =来求解方程组1k k Au v -=.反幂法地计算步骤为：步骤1 预先取定非零向量00u v =；给定精度12,εε；取(0)0m a x ()nu μ=； 步骤2 对矩阵A 作LU 分解,A LU =；令1=k ；步骤3 求解方程组1k k Ly v -= ,k k Uu y = 得到迭代向量k u ； 步骤4 规范化max()kk k u v u =步骤5 若11k k v v ε--<且()(1)2||k k n n μμε--<,则k v 即为A 地对应于n λ地近似特征向量,()1k nμ即为n λ地近似值；否则,令1+=k k ,转步骤3继续迭代.7.3 矩阵地两种正交变换本节先介绍镜面(初等)反射变换和平面旋转变换,它们是QR 算法和Jacobi 算法地基础.7.3.1 豪斯荷尔德（House holder ）变换定义7-2 设有方阵B ,若当1i j >+时(,1,2,,)i j n =,0ij b =,则称B 是上Hessenberg 矩阵,即1112121222,1n n n n nn b b b b b b B b b -⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦定义7-3 设向量ω满足21ω=,矩阵2T H E ωω=- （ω是列向量）称为初等反射矩阵,又称House holder 矩阵,记为()H ω,即211212212221212222122()2212n n n n n H ωωωωωωωωωωωωωωωω⎡⎤---⎢⎥---⎢⎥=⎢⎥⎢⎥---⎢⎥⎣⎦其中(1,2,,)i i n ω=是ω地分量.可以证明初等反射阵是对称阵()T H H =.正交阵()T H H E =. 例7-8 设向量0α≠,试证矩阵222TH E ααα=- 是一个初等反射阵. 证明 令2αωα=,则 222221||||||||1αωααα=== 由定义7-3,2222TTH E E ααωωα=-=-是初等反射阵.定理7-5 设,x y 是两个不相等地n 维列向量,且22||||||||x y =,则存在一个初等反射阵H,使得Hx y =证明 令2||||x yx y ω-=-,由例7-8可知22()()22||||T T Tx y x y H E E x y ωω--=-=-- 是一个初等反射阵.由于22||||()()T T T T Tx y x y x y x x y x x y y y -=--=--+ 注意到22||||||||x y =,即T T x x y y =,又()T T T T x y x y y x == ,故22||||2()T Tx y x x y x -=-从而22()()2||||T T x y x x y x Hx x x y --=--y y x x =--=)(. 例7-9 设1(1,2,2),(1,0,0)T T x e ==,用Householder 变换将x 化为与1e 同方向地向量.解 因为2||||3x =,可设13y e =,则22||||||||x y = 取21,1,1)||||T x y w x y -==--,构造Householder 矩阵[]11122212111,1,12123311221T H E ww -⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥=-=--=-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥-⎣⎦⎣⎦⎣⎦则13Hx e =推论 设向量12(,,,)0T n x x x x =≠,12()||||r sign x x =,且1x r ε≠-,则存在初等反射阵1222||||T T uu H E E uu u ρ-=-=- 使1Hx r ε=- .其中,1(1,0,,0)T ε=,1u x r ε=+,22||||/2u ρ=.设12(,,,)T n u u u u =,则12(,,,)T n u x r x x =+22222122222112111||||[()]221(2)2()n n u x r x x r rx x x x r r x ρ==++++=+++++=+引入初等反射阵地目地,是设法用一系列初等反射阵将原始矩阵约化成上Hessenberg 阵.由于约化过程是逐列进行地,我们先给出计算Hx 地算法步骤,该算法算出H 及r ,使Hx r ε=-,u 地分量冲掉x 地分量.（1）1max ||i i nx η≤≤=；（2）(1,2,,)ii i x x u i n η←==,此步规范化是为避免计算r 时产生溢出；（3） 12211()()nii r sign x x ==∑；（4）11u u r ←+；（5） 1ru ρ=； （6） r r η←；于是初等反射阵1T H E uu ρ-=-,1Hx r ε=-.如果要将H 作用于矩阵A ,设i a 是A 地第i 列向量,则12(,,,)n A a a a =,12(,,,)n HA Ha Ha Ha = 其中,11()()(1,2,,)T T i i i i Ha E uu a a u a ui n ρρ--=-=-=.下面讨论用初等反射阵约化原始矩阵A 称为上Hessenberg 阵地步骤.11121(1)(1)2122211121(1)(1)212212n n n n nn a a a a a a a A A A a a a a a ⎡⎤⎢⎥⎡⎤⎢⎥===⎢⎥⎢⎥⎣⎦⎢⎥⎣⎦步骤1 不妨设(1)210a ≠（否则这一步不需约化）,选择初等反射阵1R ,使(1)12111R a r ε=-,其中： 1(1)(1)2212112(1)1211112111121211111()(())(1)1()2ni i T r sign a a u a r n u r r a R E u u εερρ=-⎧=⎪⎪⎪=+-⎨⎪==+⎪⎪=-⎩∑是维单位坐标列向量 令11100U R ⎡⎤=⎢⎥⎣⎦则(2)(2)(2)(1)111213111212111(2)(2)(1)(1)222312112210A a A a A R A U AU a A R a R A R ⎡⎤⎡⎤===⎢⎥⎢⎥⎣⎦⎣⎦其中,(2)11A 是21⨯阵,(2)22a 是2n -维列向量,(2)23A 是2n -阶方阵.步骤k 设对A 已进行了1k -步约化,即111(2)()()()()11121,111,11(2)()()()1222,12,2()()()1,1,()()()1,1,11,()()(),1(2,3,,1)k k k k k k k k k k k n k k k k kn k k k k kk k k k n k k k k k k k k nk k k nkn k nnA U A U k n a a a a a a r a a a a r a a a a a a a a a ----+--++++++==-⎡⎢-⎢=-⎣()()()111213()()22230k k k k k A a A a A ⎤⎥⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎦⎡⎤=⎢⎥⎣⎦其中,()11k A 是(1)k k ⨯-阵,()22k a 是n k -维列向量,()23k A 是n k -阶方阵.设()220k a ≠,选初等反射阵()k R n k -阶,使()221k k k R a r ε=-,其中1ε是n k -维单位坐标向量,可得1()()221,1()221()1,1()(())()nk k k k k ik i k k kk k k k kk nT k k k k r sign a a u a r r r a R E u u ερρ+=++-⎧=⎪⎪⎪=+⎨⎪=+⎪⎪=-⎩∑ 令 00k k E U R ⎡⎤=⎢⎥⎣⎦则 ()()()1112131()()2223()()()111213()12300k k k k k k k k k k k k k k k k k k k k k A a A R A U A U R a R A R A a A R r R A R ε+⎡⎤==⎢⎥⎣⎦⎡⎤=⎢⎥-⎣⎦ 可见1k A +地左上角1k +阶子阵为上Hessenberg 阵,从而约化又进了一步.重复此过程,直到122112211(2)122(3)233(1)1n n n n n nn A U U U AU U U a r a r a r a -----=⨯⨯⨯⎡⎤⎢⎥-⨯⨯⎢⎥⎢⎥=-⨯⎢⎥⎢⎥⎢⎥-⎣⎦使原始矩阵A 在一系列初等反射阵地作用下,约化为上Hessenberg 阵.综上所述,有定理7-6.定理7-6 如果A 是n 阶实矩阵,则存在初等反射阵122,,,n U U U -,使221122n n U U U AU U U C --=（上Hessenberg 阵）例7-10 试证矩阵A 与其约化成为地上Hessenberg 阵C 有相同地特征值.证明 记221n P U U U -=,由于初等反射阵是正交对称阵,故122T n P U U U -=,且P 是正交阵,故T PAP C =.于是||||||||||||T T C E PAP E P A E P A E λλλλ-=-=-=-其中T PP E =,||||1T P P =.这表明A 与C 具有相同地特征多项式,即两者有相同地特征值.进一步,设y 是C 地对应于特征值λ地特征向量,即Cy y λ=,则有T PAP y y λ= ()()T T A P y P y λ=这表明T P y 为A 地对应于λ地特征向量,于是求原始矩阵A 地特征值与特征向量可转化为求上Hessenberg 阵C 地特征值和特征向量.定理7-7 若A 是实对称矩阵,则存在初等反射阵122,,,n U U U -使2211221112211()n n n n n U U U AU U U c b b c b C b b c ----⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦对称三对角阵 证明 由定理7-6,存在初等反射阵可使A 约化为上Hessenberg 阵C ,当A 是对称矩阵时,C 亦为对称阵,即T C C =,且T C 亦为上Hessenberg 阵,故C 是对称三对角阵.例7-11 用豪斯荷尔德方法将下述矩阵化为上Hessenberg 阵.1437232427A A ---⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦解 （1）对1k =,确定变换阵111000U R ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,(1)2124a ⎡⎤=⎢⎥⎣⎦ 其中1R 为初等反射阵,使(1)121110R a r ⎡⎤=-⎢⎥⎣⎦(1)12124.472136r a ==≈(1)12111 6.472136244u a r ε⎡⎡⎤+=+=≈⎢⎢⎥⎣⎦⎣⎦11121()2)28.94427r r a ρ=+≈[]1111110 6.4721361 6.472136401428.944270.4472070.8944230.8944230.447216TR E u u ρ-=-⎡⎤⎡⎤=-⎢⎥⎢⎥⎣⎦⎣⎦--⎡⎤=⎢⎥-⎣⎦（2）计算(1)122R A .记(1)221232(,)27A a a ⎡⎤==⎢⎥⎣⎦,于是 (1)1221112 3.1304967.155419(,) 1.788855 1.341640R A R a R a --⎡⎤==⎢⎥-⎣⎦其中,111111111()()(1,2)T T i i i i R a E u u a a u a u i ρρ--=-=-=（3）计算(1)121A R 及(1)1221()R A R ,即求 1(1)121211(1)1223373.1304967.1554191.788855 1.341640T T T b A R b R R R A b ⎡⎤--⎡⎤⎡⎤⎢⎥⎢⎥==--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-⎣⎦⎣⎦7.6026340.4472127.800030.3999990.399999 2.200000-⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦其中,11111()(1,2,3)T T T Ti i i b R b b u u i ρ-=-=（4）计算2111A U AU =.(1)12121(1)1221447.6026340.4472124.4721367.8000030.39999900.399999 2.2000000A R A r R A R ⎡⎤--⎡⎤⎢⎥⎢⎥⎢⎥==--⎢⎥⎢⎥-⎢⎥-⎢⎥⎣⎦⎢⎥⎣⎦为上Hessenberg 阵.7.3.2 平面旋转变换 定义7-4 称矩阵111(,)111i j csi P i j scj ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥-⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦第列第列第行第行 为平面旋转矩阵,又称Givens 矩阵,其中cos c θ=,sin s θ=.平面旋转阵(,)P i j 是一个正交阵,与单位阵只有在(,),(,),(,i i i j j j j i四个位置上地元素不一样,用其左乘矩阵A 只改变A 地第i 行和第j 行元素.设12(,,,)T n x x x x =则平面旋转变换Px y =地结果为⎪⎩⎪⎨⎧≠=+-=+=ji k x y cx sx y sx cx y k kj i j j i i ,若令/i c x =,j s x =, 则平面旋转变换向量y 地第i个分两为22j i x x +,第j 个分量为0,其余分量即为x 对应地分量.和初等反射变换一样,用平面旋转变换也可以将一个方阵化为上Hessenberg 矩阵,也可以将将一个方阵化为上三角矩阵.7.4 QR 算法7.4.1 矩阵地QR 分解定理7-8 设A 是可逆矩阵,则存在正交矩阵121,,,n P P P -使121()n P P P A R -=上三角矩阵且R 地主对角元素0(1,2,,1)ii r i n >=-.证明 若10(2,3,,)j a j n ==,则A 地第一列不需约化.若有某个 10(2)j a j n ≠≤≤,则可选择1(1,)j P j P =使A 地第一列中第j 个元素变为零.一般地,可设平面旋转矩阵12131,,,n P P P ,使(2)(2)11121(2)(2)222113122(2)(2)200nn nn nn r a a a a P P P A A a a ⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦记111312nP P P P =,则12P A A =.同理,若(2)20(3,4,,)j a j n ≠=,可选取23242,,,n P P P 使(2)(2)(2)1112131(3)(3)22232(3)(3)2212323333(3)(3)3nn n n n n nn r a a a r a a P P P A A a a a a -⎡⎤⎢⎥⎢⎥⎢⎥==⎢⎥⎢⎥⎢⎥⎣⎦记2223nP P P =,则213P P A A =.重复上述过程,可得一系列正交阵121,,,n P P P -使11121222121n n n nn r r r r r P P P A R r -⎡⎤⎢⎥⎢⎥==⎢⎥⎢⎥⎣⎦ 定理7-9 （矩阵地QR 分解）如果n 阶实矩阵A 可逆,则A 可分解为一正交阵Q 和上三角阵R 地乘积,即A QR =,且当R 地对角元素都为正数时分解唯一.证明 由定理8知存在正交阵11,,n P P -使121n P P P A R -=为上三角阵,记121T n Q P P P -=,于是T Q A R =由于(1,2,,1)i P i n =-是正交阵,则T Q 亦为正交阵,故A QR =. 若A 有两种QR 分解,记为1122A Q R Q R ==其中12,R R 为上三角阵且主对角元素都为正数,12,Q Q 为正交阵,于是12121T Q Q R R -=注意121R R -是上三角阵地乘积,结果仍为上三角阵,而12,TQ Q 是正交阵,所以121R R -也应是正交阵.若记121D R R -=,由其上三角性T D 应是下三角阵,1D -应是上三角阵；由其正交性由1T D D -=,故D 只能是对角阵,且有2T D D D E ==.又因12,R R 地主对角元素都为正数,即有222212diag[,,,]diag[1,1,,1]n D d d d E ===故1(1,2,,)i d i n ==,则D E =,于是12R R =,12Q Q =.例7-12 求矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=212240130A 地QR 分解. 解 方法1：利用初等反射阵进行QR 分解令(0)1(0,0,2)T a =,取(0)112||||2d a ==,则)2,0,2(81211)0(111)0(11-=--=e d ae d a u1110012010100TH E u u ⎡⎤⎢⎥=-=⎢⎥⎢⎥⎣⎦,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=1302402121A H 再令(0)2(4,3)T a =,取(0)222||||5d a ==,则(1)2212(2)22121,3)||||T a d e u a d e -==--2224312345TH E u u ⎡⎤=-=⎢⎥-⎣⎦令2210014305534055H H ⎡⎤⎢⎥⎢⎥⎡⎤⎢⎥==⎢⎥⎢⎥⎣⎦⎢⎥⎢⎥-⎣⎦于是21212051002H H A R ⎡⎤⎢⎥=-=⎢⎥⎢⎥-⎣⎦故123405521243005155002100T TA H H R ⎡⎤-⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥==-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦方法2：利用平面旋转阵进行QR 分解. 取1202,0100221221=+==+=s c ,则130********T ⎡⎤⎢⎥=⎢⎥⎢⎥-⎣⎦,132********T A ⎡⎤⎢⎥=-⎢⎥⎢⎥--⎣⎦再取53)3(43,54)3(44222222-=-+-==-+=s c ,则231004305534055T ⎡⎤⎢⎥⎢⎥⎢⎥=-⎢⎥⎢⎥⎢⎥⎣⎦,2313212051002T T A R ⎡⎤⎢⎥=-=⎢⎥⎢⎥-⎣⎦ 故13233405521243005155002100T T A T T R ⎡⎤-⎢⎥⎡⎤⎢⎥⎢⎥⎢⎥==-⎢⎥⎢⎥⎢⎥-⎢⎥⎣⎦⎢⎥⎢⎥⎣⎦例7-13 求矩阵⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡-=110133044A 地QR 分解,使得R 地对角线元素为正数.解 A A =1地第一列T x ]0,3,4[1=,521=x .用1x 构造镜面反射阵1H ,使得T y x H ]0,0,5[111==,令T y x u ]0,3,1[111-=-=,有⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-=-=10005453053542221111u u u E H T ,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-==11054005355112A H A 2A 地第2列对角线以下为T x ]1,0[2=,122=x .用2x 构造镜面反射阵2~H ,使得T y x H ]0,1[~222==,令T y x u ]1,1[222-=-=,易得 ⎥⎦⎤⎢⎣⎡=-=01102~222222u u u E H T,⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=⎥⎦⎤⎢⎣⎡=010100001~122H H 于是有R A H A =⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡-==54001105355333,⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡-==010540535305421H H Q容易验证,QR A =.请读者用平面旋转变换对本例地矩阵A 进行QR 分解.7.5.3 QR 算法QR 算法就是利用QR 分解构造一个矩阵序列{}k A ,当k 充分大时,k A 是近似地上三角矩阵,而该上三角阵地对角元素便是原始矩阵A 地全部特征值.设1()n n ij n n A A a R ⨯⨯==∈,对A 做QR 分解,即A QR =其中R 为上三角阵,Q 为正交阵.利用这个分解可得新矩阵（对QR 交换乘积）2T A RQ Q AQ == 由于2A 是1A 经过正交相似变换得到地,因此2A 与1A 有相同地特征值.再对2A 做QR 分解,按上述方式又可得新矩阵3A ,且3A 与2A 也具有相同地特征值.具体地说,其步骤为：设1A A =,做QR 分解111A Q R =求矩阵211111T A R Q Q A Q ==求得k A 后对k A 作QR 分解k k k A Q R =求矩阵1Tk k k k k k A R Q Q A Q +==只要A 可逆,由定理9可知,按上述方法可唯一确定矩阵序列{}k A ,且序列中任意k A 与原始矩阵有相同特征值.因此只要恰当选择正交相似变换阵12,,,k Q Q Q ,使1111111T T TT TT T k k k k k k k k k k k k k k A Q A Q Q Q A Q Q Q Q Q A Q Q Q +----====当k →∞时,逼近一个上三角阵,便可求出A 地全部特征值（为所逼近上三角阵地主对角元素）.可见,QR 算法地关键在于选择正交变换阵(1,2,)k Q k =.从定理7-8地证明看到,正交变换阵k Q 是一系列平面转换矩阵地乘积,这些平面旋转矩阵是用来将k A 地主对角线以下元素约化为零地.如果将QR 算法直接应用于原始矩阵,计算量很大,所以在实际计算中,总是先将原始矩阵用豪斯赫尔德方法约化为上Hessenberg 阵,而后再对上Hessenberg 阵应用QR 算法.可以证明,由上Hessenberg 阵用QR 算法生成地矩阵序列中地每个矩阵仍为上Hessenberg 阵.7.5 雅可比方法雅可比方法是用来计算实对称矩阵地全部特征值及特征向量地一种有效方法.它地基本思想是,通过一组正交相似变换对称矩阵A 化为对角矩阵,得其全部特征值.定理7-10 设A 为n 阶对称矩阵,T C PAP =,其中P 为正交矩阵,则22||||||||F F C A = 证明 一方面2222111||||()()nnnFiji i j i A a tr A A λ======∑∑∑另一方面2221||||()()()nTFi i C tr C C tr C C λ====∑由假设()()i i A C λλ=,故22||||||||F F C A =.设n n A R ⨯∈为对称矩阵,(,)P i j 为一平面旋转矩阵,则T C PAP =（其中()ij n n C c ⨯=）地元素计算公式为：（1）22cos sin 2sin cos ii ii jj ij c a a a θθθθ=++22sin cos 2sin cos jj ii jj ij c a a a θθθθ=+-（2）1()sin 2cos 22ij ji jj ii ij c c a a a θθ==-+ （3）第i 行元素和第j 列元素cos sin (,)ik ki ik jk c c a a k i j θθ==+≠ （4）第j 行元素和第i 列元素 cos sin (,)jk kj jk ik c c a a k i j θθ==+≠（5）(,,)lk lkc a l k i j =≠这说明,当A 经过一初等正交相似变换化为C 时,只需按上述公式计算C 地第i 列.第j 列元素,由对称性可得第i 行和第j 行元素,C 地其余元素与A 地对应元素相同.设A 地非对角元素0ij a ≠,我们可选择平面旋转阵(,)P i j ,使T C PAP =地非对角元素0ij ji c c ==.由定理11可选择(,)P i j ,使sin 2cos 202jj iiij ji ij a a c c a θθ-==+=即选择θ,使22(||)4ij ii jja tg a a πθθ=≤-其中定理7-11 设n n A R ⨯∈为对称阵,0ij a ≠为A 地一个非对角元素,则可选择一平面旋转阵(,)P i j ,使T C PAP =地非对角元素0ij ji c c ==且T C PAP =与A 地元素满足下述关系（1）2222(,)ik jk ik jkc c a a k i j +=+≠（2）222222ii jj ii jj ij c c a a a +=++ （3）22(,,)iklk c a l k i j =≠证明 由上面地计算ij c 公式直接计算可知（1）成立.由（1）及定理7-10可证（2）.如果用()S A 表示A 地非对角线元素地平方和,()D A 表示A 地对角线元素平方和,则2()()2ijD C D A a =+ ,2()()2ij S C S A a =- 这说明C 地对角线元素平方和比A 地对角线元素平方和增加了22ij a ,C 地非对角线元素平方和比A 地非对角线元素平方和减少了22ij a .下面介绍雅可比方法.首先在A 地非对角元素中选择绝对值最大地元素（称为主元素）,如11||max ||i j lk l ka a ≠=可设110i j a ≠,否则A 已经对角化了.由定理12,选择一平面旋转矩阵111(,)P i j ,使111TAP AP =地非对角元素11110i j j i c c ==. 再选(1)1()lkn n A a ⨯=地非对角元素中地主元素,如 22(1)(1)||max ||0i j lk l ka a ≠=≠由定理12,又可选择一平面旋转矩阵222(,)P i j ,使2212T A P A P =地非对角元素2222(2)(2)i j j i a a ==（注意上次消除了地主元素这次又可能变为不是零）. 继续这个过程,连续对A 实行一系列平面旋转变换,消去非对角线绝对值最大地元素,直到将A 地非对角元素全化为充分小为止,从而求得A 地全部（近似）特征值.定理7-12 （雅可比方法地收敛性）设()ij n n A a ⨯=为实对称矩阵,对A 施行上述一系列平面旋转变换1(1,2,)Tm m m mA P A P m -==则 lim ()m m A D→∞=对角矩阵证明 记()()m m lk n n A a ⨯=,()2()m m lk l kS a ≠=∑由定理7-11地（2）可得()212()m m m ij S S a +=-其中 ()()||max ||m m ijlk l ka a ≠= 又由于()2()2()(1)()m m m lk ij l kS a n n a ≠=≤-∑即()2()(1)m m ij S a n n ≤- 由以上得12(1)(1)m m S S n n +≤-- 反复应用上式,即得1102(1)(2)(1)m m S S n n n ++≤->-故 lim 0m m S →∞= 可以证明()lim m ll m a →∞存在(1,2,,)l n =. 下面介绍特征向量地计算.由雅可比收敛定理知,当m 充分大时2112T TTmm P P P AP P P D ≈记12T T T T m m R P P P =,则T m R 地列向量就是A 地近似特征向量.计算Tm R 可采用累积地办法,用一数组R 保存Tm R ,开始时R E ←,以后对A 每进行一次平面旋转变换,就进行计算Tm R RP ←用初等正交阵T m P 右乘R 只需计算R 地两列元素,若记(,)m m P P i j =,则Tm RP 地计算公式为()()cos ()sin (1,,)()()sin ()cos li li lj li li lj l n θθθθ←+⎧⎪=⎨←+⎪⎩R R R R R R关于sin θ和cos θ地计算如下.由定理7-11知,当0ij a ≠时,可选θ满足2tg2ij ii jja a a θ=-方ii jj a a ≠时,由22tg 1tg21tg dθθθ=≡- 得到tg θ地二次方程2tg 2tg 10d θθ+-=解得tg θ=选取tg 0d d θ>=<由此得 |tg |1θ≤可由集合{},,ii jj ij a a a 来计算sin ,cos θθ,设0,||max ||ij ij lk l ka a a ≠≠=,则210tg ,()10cos sin cos ii jj ija a d a d t s d d c t ct sθθθθ-⎧=⎪⎪⎪≥⎧⎪=≡=⎨-<⎨⎩⎪⎪=≡⎪⎪=⋅=≡⎩如果jj ii ij a a a -<<,则12ij ii jja t d a a ≈=-,将c,s 代入定理7-9地（1）中可得ii ii ij jj jj ij ij ji c a ta c a ta c c ⎧=+⎪=-⎨⎪==⎩ 每迭代一次地主要工作是选m A 地非对角线元素中地主元素与计算T 111m m m +++=A P AP .首先计算sin ,cos ,θθ,只需计算1m +A 地第i 列,第j 列元素,再算对称元素,不用做3个矩阵地乘法.计算机计算时,需要两组工作单元,以便存储A （或m A ）和R .可用()2()m m lk l ka ε≠=<∑S 控制迭代终止,其中ε是要求地精度.例7-14 用雅可比方法计算对称阵210121012⎡⎤-⎢⎥--⎢⎥⎢⎥-⎣⎦A = 地特征值.解 第1步0=A A ,选非对角线元素中地主元素121(1,2)a i j =-==0,1,1/0.7071068,1/0.7071068d t c s ======T 111100.7071068030.70710680.70710680.70710682⎡⎤-⎢⎥==-⎢⎥⎢⎥--⎣⎦A P AP第2步 在1A 中选非对角元素地主元素(1)130.7071068(1,3)a i j =-==0.7071068,0.5176381,0.8880738,0.4597008d t c s ====T 22120.63397460.325057600.325057630.627963000.62798302.366025-⎡⎤⎢⎥=--⎢⎥⎢⎥-⎣⎦A P A P 第3步 在2A 中选非对角元素地主元素(2)230.627930(2,3)a i j =-==0.5047869,0.6153960,0.8516540,0.5241045d t c s =-=-==-T 33230.63397460.27683660.17036420.27683663.38644600.170364201.979579⎡⎤--⎢⎥=-⎢⎥⎢⎥-⎣⎦A P A P 第4步 在3A 中选非对角元素地主元素(3)120.2768366(1,2)a i j =-==4.971292,0.09958013,0.9950785,0.09909004d t c s ====T 44340.606407200.169525803.4140130.016881400.16952580.016881401.979579⎡⎤-⎢⎥=⎢⎥⎢⎥-⎣⎦A P A P 第5步 在4A 中选非对角元素地主元素(4)130.1695258(1,3)a i j =-==4.050038,0.1216293,0.9926842,0.1207395d t c s ==== 2T 255450.58578790.20382521000.203825210 3.4140130.0167579000.016757902.000198--⎡⎤⨯⎢⎥=⨯⎢⎥⎢⎥⎣⎦A P A P 于是A 地特征值为1233.414013, 2.000198,0.5857879λλλ===A 地精确特征值为12(1 3.414214λ=≈,22λ=,32(10.585786λ=-≈ 且可逐步求出412345T T T T T T R P P P P P =地列向量,即得A 地近似特征向量.雅可比方法是一个求对称矩阵A 地全部特征值及特征向量地迭代方法,精确度较高,但计算量较大,对稀松带状矩阵经过平面旋转变换后其稀松带状将被破坏,所以很少使用.习题71．设911203111(2102113810A j j B ⎡⎤⎡⎤⎢⎥⎢⎥===⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦试估计它们地特征值所在地范围.2．编写幂法程序,并求矩阵732341213A -⎡⎤⎢⎥=⎢⎥⎢⎥--⎣⎦地主特征值及对应地特征向量（准确到小数点后3位）.3．若p 是A 地特征值j λ地一个近似值,且||||()j i p p i j λλ-<-≠则1j pλ-是1()A pE --地主特征值.试用反幂法求矩阵134231111A ⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦地最接近于6地特征值及对应地特征向量.4．设有向量(2,1,2)Tx=,试构造初等反射阵H,使(3,0,0)THx=.5．设(2,3,0,5)Tx=,(1,0,0,5)Te=,用Householder变换化x为与e同方向向量.6．设031042212A⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦,求其QR分解.7．设221022212A⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦,求其QR分解.8．利用初等反射阵将134312421A⎡⎤⎢⎥=⎢⎥⎢⎥⎣⎦正交相似约化为对称三对角阵.9．试用平面旋转变换阵对矩阵A作QR分解,其中111021245A⎡⎤⎢⎥=-⎢⎥⎢⎥-⎣⎦.10．按下列要求编写程序框图.（1）将一般矩阵用豪斯赫尔德方法约化称上Hessenberg阵.（2）对矩阵作QR分解.（3）对上Hessenberg阵应用QR算法求全部特征值及相应地特征向量.11．用QR算法求矩阵120211013A⎡⎤⎢⎥=-⎢⎥⎢⎥⎣⎦地全部特征值.12．设A是对称矩阵,λ和(1)x x=是A地一个特征值及相应地特征向量.又设p是一个正交阵,使1(1,0,0,,0)Tpx e==证明T=是第一行和第一列除了λ外,其余元素均为零.B PAP。
数值分析课件第9章
同法解得
y(0.4) y2 2.020118, y(0.8) y4 2.8565830,
y(0.6) y3 2.451578 y(1.0) y5 3.243224
工科研究生公共课程数学系列
机动 上页 下页 首页 结束
4、单步法的局部截断误差与阶
单步法的一般表示形式
yn1 yn h(xn , yn , yn1, h) 其中与f (x, y)有关，称为增量函数,当含有yn1时，方法是隐式的，
f
( xn
ih,
y( xn
ih))
或表示为 yn1 yn h(xn , yn , h)
其中
r
(xn, yn, h) ci Ki i 1
K1 f (xn , yn )
i 1
Ki f (xn ih, yn h ij K j ) j 1
这里ci , i , ij均为常数。称为r级显式龙格 - 库塔(Runge Kutta)法，
简称R - K方法。
工科研究生公共课程数学系列
机动 上页 下页 首页 结束
2、二阶显式R-K方法
r 2 的R - K方法, 计算公式如下 :
yn1 yn h(c1K1 c2K2 ) K1 f (xn , yn )
K2 f (xn 2h, yn 21hK1)
这里c1, c2 , 2 , 21均为待定常数。
上的近似值 y1, y2 , , yn , yn1, 。相邻两节点间的间距 hn xn1 xn称为步长。假定hi h为定数，这时节点 为xn x0 nh, n 0,1, 2,
初值问题数值解法的基本特点：它们都采用 “步进式”，即求解 过程顺着节点排列的次序一步一步地向前推进。描述这类算法， 只要给出已知信息yn , yn-1, yn-2, 计算yn1的递推公式。
西安科技大学 研究生 数值分析课件
➢优化计算(jì suàn)步骤，提高计算(jì suàn)效率。 ➢严格控制递推公式中误差的传播
简化
(jiǎnhuà)
计 算步骤
n
Pn (x) ak xk x(x((an x an1) an2 ) a1 a0 k 0
1.4.2 算法的稳定性与病态数学问题 1 算法的稳定性
➢ 算法的数值稳定性是指算法对误差的传播（或积累）是否受 到控制(kòngzhì)的问题。如果算法的计算结果对初始数据的误 差及计算中的舍入误差不敏感，则称该算法是稳定的；否 则该算法是不稳定的。
➢ 由于误差不可避免，算法的稳定性在数值计算中是不可回避 的重要问题。
z* e5 0.0067379 , z 0.00673
共二十页
1.2 误差(wùchā)与有效数字
1.2.2 误差的度量
3. 有效数字
• 这种定义方法实际上给出了有效数字与绝对误差的关系。下
面的定理(dìnglǐ)揭示了有效数字与相对误差的关系。
• 定理1-1 设近似值x可写成（1-1）的规格化形式，若x至少有n
1.2.3 误差的传播(chuánbō)
1. 函数的误差
f (x*) f (x) f (x)(x* x)
例1-5
( f (x)) f (x) (x)
r ( f
(x))
( f (x))
f (x)
f (x) (x)
f (x)
P13习题(xítí)3 (1 5)
f (x*, y*) f (x, y) f (x, y) (x* x) f (x, y) ( y* y)
研究生数值分析 讲义
第一章绪论上次课要点：§1 数值分析的几个基本问题一、用数学方法解决科学与工程问题的步骤二、研究对象三、研究内容四、研究数值计算方法的意义五、算法设计的基本思想六、算法应具备的特性§2 数值计算的误差一、误差的分类 1．截断误差 2．舍入误差二、误差的概念 1．绝对误差x x x E -=**)(2．相对误差xx E x E r )()(**=（其中0≠x ）本次课继续。
三、数值运算的误差当自变量有误差时，一般地，其函数值也有误差。
误差——可能是截断误差——也可能是舍入误差1．一元函数的误差设*x 是准确值x 的近似值，则函数)(x f 的近似值为)(*x f 。
由于))(()()(**x x f x f x f -'=-ξ，ξ介于x 与*x 之间，所以)()()()(**x x f x f x f -'=-ξ从而)()())((**x f x f εξε'≈2. 多元函数的误差对于多元函数),,,(21n x x x f ，设自变量的近似值分别为**2*1,,,nx x x ，则),,,(),,,(),,,((21**2*1**2*1n nnx x x f x x x f x x x f E -=)(|)(|*),,,(*1),,,(1**2*1**2*1nx x x nx x x x e x f x e x f nn∂∂++∂∂≈于是误差限),,,((**2*1nx x x f ε∑=∂∂≈nk kx x x kx x f n1*),,()(|**2*1ε特别)()()(*2*1*2*1x x x x εεε+≈±)()()(*1*2*2*1*2*1x x x x x x εεε+≈2*2*1*2*2*1*2*1)()()/(xx x x x x x εεε+≈四、病态问题与条件数一个工程或科学计算问题：——往往需要巨量的机器运算——每次运算都可能产生误差——这些误差有正有负，绝对值有大有小误差积累的结果很难定量分析。
数值分析ppt课件
数值积分与微分
数值积分
通过数值方法近似计算定积 分，如梯形法则、辛普森法 则等。
数值微分
通过数值方法近似计算函数 的导数，如差分法、中心差 分法等。
常微分方程的数值解法
通过数值方法求解常微分方 程，如欧拉方法、龙格-库塔 方法等。
03
数值分析的稳定性与误差分析
误差的来源与分类
模型误差
由于数学模型本身的近 似性和简化，与真实系
非线性代数方法
非线性方程组的求解
通过迭代法、直接法等求解非线性方程组，如牛顿法、拟牛顿法 等。
非线性最小二乘问题
通过迭代法、直接法等求解非线性最小二乘问题，如GaussNewton方法、Levenberg-Marquardt方法等。
多项式插值与逼近
通过多项式插值与逼近方法对函数进行近似，如拉格朗日插值、 样条插值等。
机器学习与数值分析的交叉研究
机器学习算法
利用数值分析方法优化和改进机器学 习模型的训练和预测过程，提高模型 的准确性和效率。
数据驱动的模型
通过数值分析方法处理大规模数据集 ，提取有用的特征和模式，为机器学 习模型提供更好的输入和输出。
大数据与数值分析的结合
大数据处理
利用数值分析方法处理和分析大规模数 据集，挖掘其中的规律、趋势和关联信 息。
数值分析PPT课件
contents
目录
• 引言 • 数值分析的基本方法 • 数值分析的稳定性与误差分析 • 数值分析的优化方法 • 数值分析的未来发展与挑战
01
引言
数值分析的定义
数值分析
数值分析是一门研究数值计算方法及 其应用的学科，旨在解决各种数学问 题，如微积分、线性代数、微分方程 等。
研究生数值分析课件ch
数值分析是数学的一个重要分支，主要研究如何利用数值方法求解数学问题和近似计算 实际问题的数值解。它为科学研究、工程技术和实际应用等领域提供了重要的数学工具。 数值分析的重要性在于它能够将许多抽象的数学概念和理论转化为具体的数值计算方法，
使得我们能够更加方便地解决各种复杂的实际问题。
数值分析的应用领域
在金融领域，数值分析也被 广泛应用于风险评估、投资 组合优化、期权定价等方面 。通过数值分析的方法，我 们可以更加准确地评估投资 风险和收益，从而做出更加 明智的决策。
数值分析的发展历程
总结词
数值分析的发展历程可以追溯到上世纪初，随着计算 机技术的不断发展，数值分析的理论和方法也在不断 更新和完善。
05
数值积分与微分
牛顿-莱布尼兹公式与复化求积法
牛顿-莱布尼兹公式
该公式是微积分中的一个基本定理，用于计算定积分。 通过将积分区间分成若干小区间，并在每个小区间上应 用微积分基本定理，再利用定积分的线性性质进行求和 ，最后取极限得到定积分的值。
复化求积法
当被积函数是复杂函数或者积分区间是复杂形状时，直 接应用牛顿-莱布尼兹公式可能会遇到困难。此时，可以 采用复化求积法，即将积分区间分成若干个小区间，然 后在每个小区间上应用牛顿-莱布尼兹公式，最后将所有 的结果相加得到定积分的近似值。
改进欧拉法
为了提高欧拉方法的精度，可以对欧拉方法进行改进。一种常见的改进方法是使用二阶 欧拉方法，它考虑了更多的函数值，从而提高了逼近的精度。
龙格-库塔方法
龙格-库塔方法是一种高阶数值方法，用于求解常微分方程。它基于泰勒级数的思想，通过迭代的方式逐步逼近方程的精确解 。与欧拉方法相比，龙格-库塔方法具有更高的精度和更好的稳定性。
《数值分析教程》课件
一种适用于大规模计算的数值方法
详细描述
谱方法适用于大规模计算，通过将问题分解为较小的子问 题并利用多线程或分布式计算等技术进行并行计算，可以 有效地处理大规模的计算任务。
感谢您的观看
THANKS
具有简单、稳定和可靠的优点。
05
数值积分与微分
牛顿-莱布尼兹公式
要点一
总结词
牛顿-莱布尼兹公式是数值积分中的基本公式，用于计算定 积分。
要点二
详细描述
牛顿-莱布尼兹公式基于定积分的定义，通过选取一系列小 区间上的近似值，将定积分转化为一系列小矩形面积之和 ，从而实现了数值积分。
复化求积公式
总结词
算机实现各种算法，为各个领域的科学研究和技术开发提供了强有力的支持。
数值分析的应用领域
总结词
数值分析的应用领域非常广泛，包括科学计算、工程 、经济、金融、生物医学等。
详细描述
数值分析的应用领域非常广泛，几乎涵盖了所有的科学 和工程领域。在科学计算方面，数值分析用于模拟和预 测各种自然现象，如气候变化、生态系统和地球科学等 。在工程领域，数值分析用于解决各种复杂的工程问题 ，如航空航天、机械、土木和电子工程等。在经济和金 融领域，数值分析用于进行统计分析、预测和优化等。 在生物医学领域，数值分析用于图像处理、疾病诊断和 治疗等。总之，数值分析已经成为各个领域中不可或缺 的重要工具。
03
线性方程组的数值解法
高斯消去法
总结词
高斯消去法是一种直接求解线性方程组的方法，通过一系列 行变换将系数矩阵变为上三角矩阵，然后求解上三角方程组 得到解。
详细描述
高斯消去法的基本思想是将系数矩阵通过行变换化为上三角 矩阵，然后通过回带求解得到方程组的解。该方法具有较高 的稳定性和精度，适用于中小规模线性方程组的求解。
研究生数值分析(1)48页
x * 有 2 位有效数字，绝对误差限为
(x*)1102 0.005
2
相对误差限为 r(x*)(xx**)409 .0100250.0102
y * 有 3 位有效数字，绝对误差限为0.00005，
相对误差限为0.00102。
（4）舍入误差：由于计算机字长有限，只 能对有限位进行运算，因而往往进行四舍五入， 这样产生的误差称为舍入误差。
误差是不可避免的，要做到与实际问题的绝 对准确，是办不到的。因此，在计算方法里讨论 的问题就是怎样尽量设法减少误差，提高精度。
在四中误差中，模型误差和观测误差是客 观存在的，截断误差和舍入误差是由计算方法和 计算工具引起的，我们在研究数学问题的数值解 法时，主要是分析讨论计算方法的截断误差和舍 入误差。
函数近似值 y * 的相对误差
e r *(y)e* y (* y) ( x f1)*x y 1 * *e r *(x 1 ) ( x f2)*x y 2 * *e r * (x2) （2）
利用（1）、（2）两式，可以得到两数 和、差、积、商的绝对误差与相对误差传播 的估计式.
e* ( x1 x2 ) e* ( x1) e* ( x2 )
一个近似值的准确程度的。
(1)绝对误差与绝对误差限：
若 x * 为准确值x的一个近似值，则称 x x *
为近似值 x * 的绝对误差，用 e * ( x ) 表示，
即 e*(x)xx*
实际问题中，由于无法知道准确值 x 因 而无法计算绝对误差的大小，只能根据具体 情况估计绝对误差的上限使
e*(x)xx* *
如果
e*(x) 110mn 2
数值分析第一章基础知识优秀课件
16 周二 3课时 第八章 常微分方程初值问题数值解法[1] 17 周二 3课时 第八章 常微分方程初值问题数值解法[2] 18 周二 3课时 习题课 19 周二 3课时 总复习
注：数值算法演示主要用Matlab和C语言实现，有时采用
Mathematica
实8/7现6 。课郑后州实大验学题201可4-用20任15何学年一硕种士计研算究生工课具程完成数值。分析 Numerical Analysis
4/76
郑州大学2014-2015学年硕士研究生课程 数值分析 Numerical Analysis
预备知识
➢ 微积分和常微分方程； ➢ 线性代数； ➢ 数值计算程序设计
（C/Matlab和Mathematica）
5/76
郑州大学2014-2015学年硕士研究生课程 数值分析 Numerical Ana.1 教学内容时间安排
周次 2 3 4 5 6 7 8 9 10 11
课次 周二 周二 周二 周二 周二 周二 周二 周二 周二 周二
课时 3课时 3课时 3课时 3课时 3课时 3课时 3课时 3课时 3课时 3课时
教学内容 第一章 基础知识 第二章 代数插值[1] 第二章 代数插值[2] 第三章 数据拟合的最小二乘法[1] 第三章 数据拟合的最小二乘法[2] 第四章 数值微分与数值积分[1] 第四章 数值微分与数值积分[2] 习题课 第五章 解线性代数方程组的直接法[1] 第五章 解线性代数方程组的直接法[2]
参考教材
教材
李庆扬,王能超,易大义.数值分析(第五版).北京:清华大学出版社,2008 李清善，宋士仓. 数值方法. 郑州：郑州大学出版社，2007.
参考资料
1.关治，陈景良. 数值计算方法. 北京：清华大学出版社，1990. 2.周铁，徐树方等. 计算方法. 北京：清华大学出版社，2006. 3.徐翠微，孙绳武. 计算方法引论. 北京：高等教育出版社，2005. 4.John H.Mathews, Kurtis D.Fink. 数值方法(MATLAB版). 北京：电子
第二章 距离空间 研究生 数值分析 教学课件(共43张PPT)
闭集：设 E 是距离空间，A E ，若 A 的补集 AEC E A 为开集，则称 A 为 E 中的闭集。
第十八页，共43页。
极限点（聚点）、导集：设 E 是距离空间，A E, x0 E ， 若在O(x0, ) 内都含有属于 A 而异于 x0 的点，则称 x0 为 A 的一个极限点（或聚点）。 A 的极限点的全体称 为 A 的导集。记作 A 。
1) 定义（收敛点列） 设 X 是一个距离空间，{x n}是
X 中点列， x X 。若 n 时, (xn, x) 0 （即 0, N, 当n N时, (xn, x) ）
则称点列 x n 在 X 中按距离 收敛于 x，记作
lim
n
xn
x
或 xn
x(n
)
此时，称 x n 为收敛点列，x 为 x n 的极限或极限点。
(x, y) max x(t) y(t) t[ a ,b ]
则C[a,b]在 下是距离空间。
若 1(x, y)
b
x(t) y(t) dt
a
， 则C[a,b]在 1下也
是距离空间
第十页，共43页。
例4 设 Lp[a,b] (P 1) 表示[a,b]上 p 方可积的所有函数的
全体，即
2) 稠密性
定义（稠密性）设 X 是距离空间， A, B X 。若
x A，总存在 B 中的点列 x n 收敛于 x, 则称 B 在 A 中
稠密。
（即x
A,
{xn}
B ，使
lim
n
xn
x
）
例 1 有理数集 Q 与无理数集 QRC 都在 R1 中稠密。
第二十三页，共43页。
《数值分析》》课件
遗传算法
模拟生物进化过程的搜索算法，通过优胜略汰 的方式找到最优解。
模拟退火法
模拟金属退火过程的搜索算法，通过随机性和 温度控制来逼近最优解。
粒子群优化
模拟粒子群行为的算法，通过粒子之间的合作 和个体经验找到最优解。
截断误差
使用有限项进行级数展开时未考虑所有无穷项导致的误差。
舍入误差
由于数学运算符的近似计算和截取，导致了计算结果与真实结果之间的差距。
插值和拟合方法
插值和拟合方法是数值分析中常用的技术，用于根据已知数据点推导出未知数据点的值或找到拟合曲线或曲面。
插值方法
利用已知数据点之间的关系推导出处于数据点之间 位置的值。
2 物理学
求解量子力学方程、天体力学模拟和粒子物 理实验结果分析。
3 金融
风险评估、期权定价和投资组合优化。
4 医学
数值模拟手术、疾病预测和药物研发。
数值分析的历史和趋势
数值分析起源于古代文明对数学问题的解决方案。如今，随着计算机技术进步，数值分析在各个领域的 应用呈指数级增长。
1
古代
古埃及的巴比伦人使用分段直线插值法求解方程。
《数值分析》PPT课件
本课程介绍《数值分析》的学习目标，定义和应用领域。深入探讨数值分析 的历史、发展和误差分析。了解插值和拟合方法，数值微积分和数值积分。
数值分析的应用价值
数值分析在工程、物理学、金融等领域扮演着重要角色。通过数值模拟和优化算法，我们能够解决复杂问题并 做出准确的预测。
1 工程
计算结构力学、流体力学和电磁场分析，优 化设计和仿真。
2
20世纪
计算机的发明使数值分析成为可能，并发展了更高精度和快速的算法。
数值分析PPT课件
03
数值分析的方法和技巧广泛应用于科学计算、工程、经 济、金融等领域。
主题的重要性
随着计算机技术的不断发展， 数值计算已经成为解决实际问 题的重要手段。
数值分析为各种数学问题提供 了有效的数值计算方法和技巧， 使得许多问题可以通过计算机 得以解决。
掌握数值分析的知识和方法对 于数学建模、科学计算、数据 分析等领域具有重要意义。
意义。
未来数值分析的发展方向
随着计算机技术的不断发展，数值分析 将更加依赖于计算机实现，因此数值算 法的优化和并行化将是未来的重要研究
方向。
随着大数据时代的到来，数值分析将更 加注重对大规模数据的处理和分析，因 此数据科学和数值分析的交叉研究将成
为一个新的研究热点。
随着人工智能和机器学习的发展，数值 分析将更加注重对非线性、非平稳问题 的处理，因此新的数值算法和模型将不
数值积分和微分
矩形法
将积分区间划分为若干个小的矩形区域，求 和得到近似积分值。
辛普森法
梯形法
利用梯形公式近似计算定积分，适用于简单 的被积函数。
利用三个矩形区域和一个梯形区域的面积近 似计算定积分。
02
01
高斯积分法
利用高斯点将积分区间划分为若干个子区间， 通过求和得到近似积分值。
04
03
矩阵的特征值和特征向量
数值分析ppt课件
目录
• 引言 • 数值分析的基本概念 • 数值分析的主要算法 • 数值分析的误差分析 • 数值分析的实例和应用 • 结论
01
引言
主题简介
01
数值分析是数学的一个重要分支，主要研究如何利用数 值计算方法解决各种数学问题。
02
它涉及到线性代数、微积分、微分方程、最优化理论等 多个数学领域。
研一 数值分析课件
数值分析Numerical Analysis数值分析是学习和了解科学计算的桥梁！数学的一种分类基础数学（理想化的）计算数学（实用化的）随机数学（圆滑的）数值分析学习方法1.注意掌握各种方法的基本原理2.注意各种方法的构造手法3.重视各种方法的误差分析4.做一定量的习题5.注意与实际问题相联系6.了解各种方法的算法与程序实现●教材与参考书1.《数值分析简明教程》，北京交通大学自编 , 教材科,20112.《数值分析》第四版，李庆杨等 , 清华大学出版社,20013.数学实验基础，王兵团，清华大学出版社，2008●考试方法1.闭卷考试（可带一张A4纸资料）第1章绪论本章主要介绍科学计算的特点、数值分析基本知识和概念，它们对学习数值分析、了解科学计算原理，以及进行科学计算都是很有帮助的。
1.1 学习数值分析的重要性思考：用一种计算机语言正确编程，计算机就一定能给出正确的结果，问题是这样简单吗？例 1.1 将数列105nn x I dx x =+⎰写成递推公式形式，并计算数列12,,I I 的值。
解：因为111011111005551555n n n n n n n x x x I dxx x x dx dx I x n-----+-=+=-=-+⎰⎰⎰ 得到计算I n 的递推公式()1151,2, 1.1n n I I n n-=-=由10016ln 55I dx x ==+⎰由递推公式（1.1）可依次算出I 1，I 2，……。
实际中，计算时一般需要具体的数据，若取0I 为准确到小数点后8位的近似值作为初始值，在字长为8的计算机上编程计算，可出现2120.3290211010I =-⨯的结果，这显然是错误的！（为什么？）用计算机解决实际问题的四个步骤1.建立数学模型；2.选择数值方法；！3.编写程序；4.上机计算。
1.2计算机中的数系与运算特点1.计算机的数系 ● 数学中的实数123100.cx a a a =±⨯⋅⋅⋅其中}{9,4,3,2,1,0⋅⋅⋅∈i a ，c 为整数。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
取对数函数 y ln(x)
C(x) 1/ ln(x)
C(0.1) 0.43,C(0.9) 0.95,C(0.99) 99.5,C(1.01) 100.5 C(10) 0.43,C(100) 0.22
x 越接近1条件数越大， 求对数函数值的相对误差越大。
1.3.3 避免误差危害的若干原则
2
n!
截断误差:
ex xn1 (n 1)!
例2 用0.3333近似代替1/3.
误差:1/3-0.3333=0.000033….
【注】 数值分析中主要关心方法误差和舍入误差。
1.2.2 误差的概念
设x*为精确值,a为x*的一个近似数。 一、绝对误差
数值分析——绪论
【定义1】 近似数a的绝对误差(误差): E(a) x* a 如果 E(a) 则δ为a的绝对误差限(误差限)
数值分析——绪论
1.2 数值计算的误差 1.2.1 误差的来源
（1）模型误差
（2）观测误差（测量误差）
（3）方法误差(截断误差)（4）舍入误差
例1 由Taylor公式求 ex 的近似值。
ex 1 x x2 xn ex xn1
2
n! (n 1)!
取n项近似: ex 1 x x2 xn
条件数C衡量问题的病态程度：
C越大病态可能越严重。
例8 y=f(x)的条件数定义为：
xf (x)
C(x) cond( f (x)) f (x)
因为
f (x) f (a) f (x)(x a) xf (x) x a
r f (a)
f (x)
f (x)
f (x)
x
cond( f (x))ra
乘法：n次 加法：n次 例11 用克莱姆（Cramer）法则求解线性方程组Ax=b
当n=20时 总计算量≈21×19×20！约307816年 天啊，30多万年的时间！ 高斯消元法：2千多次乘除运算 【注】空间复杂度：算法的复杂度、时间复杂度。
三、要避免用绝对值很小的数除数
数值分析——绪论
四、两数相加要防止大数“吃”掉小数
误差 3.1416 1 104
2
误差 3.14 1 102
2
【注】四舍五入后得到的近似数，从第一位非零数开 始直到最末位，有几位就称该近似数有几位有效数字
数值分析——绪论
【定义3】 设x*的近似值为 a 10m 0.a1a2 an a1是1到9中的一个整数，
算法2误差分析
i1
1 10
i
,i
n,
n
1,
1
0
1 10
n
n
数值稳定
蝴蝶效应
1.3.2 病态数学问题与条件数
数值分析——绪论
数值稳定性是对算法而言的。 病态是数学问题即数学模型本身的性质, 与算法无关。
病态数学问题：当输入数据有微小摄动时， 会引起解的大扰动。 相反的问题为良态数学问题。
第1章
立体化教学资源系列——数值分析
绪论
理学院应用数学系
数值分析——绪论
1.1 数值分析的研究对象与特点
数值分析(Numerical Analysis): 研究用计算机求解各种数学问题的数值计算方
法，对求得的解的精度进行评估，以及如何在计算 机上实现求解等。
数值分析特点：
（１）严格的数学理论，实用性很强。 （２）提供实际可行的且计算复杂性好的有效算法。 （３）可靠的理论分析与数值试验。
算法1是不稳定的 算法2是稳定的对于某个算法，若输入数据的误差在计 算过程中迅速增长而得不到控制，则称该算法是数 值不稳定的，否则是数值稳定的。
算法1误差分析
En*
1 n
10En*1, En
1 n
10En1
n 10 n1 (10)n0 , n 1,2, 数值不稳定
二、注意简化计算步骤，减少运算次数
数值分析——绪论
例10 计算多项式 Pn (x) an xn an1xn1 a1x a0 的值。 (1)若直接计算 ak xk 再逐项相加
乘法：n (n 1) 2 1 n(n 1)
2
加法：n (2)秦九韶算法 Pn (x) ( (an x an1)x a1)x a0
E(n)=1/n-10*E(n-1); %递推计算 end E
%用quadl计算积分值E8
E8=quadl('x.^8./ (x+10)',0,1)
n
En(算法1)
E
(算法2)
n
En*(精确值)
0
0.09531
0.0953
1
0.0469
0.0469
2
0.0310
0.0310
3
0.0233
0.0232
实际运算 Er (a) (x* a) / a
例4 a=3.14是π的近似值。
r / a
E(a) 3.14 0.002
Er
(a)
0.002
0.002 3.14
6.36942104
三、有效数字 例如 3.14159265
取3位,a=3.14,δ≤0.002 取5位,a=3.1416.δ≤0.000008
4
0.0167
0.0185
5
0.0333
0.0154
6
-0.1667
0.0131
7 1.8095
0.0115
8
-17.9702 0.0102
0.09531018 0.04689820 0.03101798 0.02315353 0.01846471 0.01535290 0.01313766 0.01148056 0.01019439
数值分析——绪论
一、避免两相近的数相减
例9 计算 1 cos2
1 cos2 1 0.9994 0.0006
只有一位有效数字
利用三角恒等式 1 cosx 2sin2 x
2
1 cos2 2sin2 1 0.000613
具有三位有效数字
如果无法改变算式，则采用增加有效位数进行运算。 在计算机上则采用双倍字长的高精度运算。
例3 绝对误差的局限性例子。 δ光速=４公里／秒，δ跑速=0.0１公里／秒 1)数量:δ光速>δ跑速 2)二者比较，光速的测量更准确。 【注】绝对误差（或误差限）不能充分说明近似数的 精确程度。
二、相对误差
数值分析——绪论
Er (a) r(相对误差限)
【定义2】 近似数a的相对误差: Er (a) (x* a) / x*
1 xn dx, n 0,1, 0 x 10
,8
得两个递推算法:
算法1
En
1 n
10En1, n
1,2,
,8
算法2
11
En1
10
( n
En ),n
8,7,
,1
数值分析——绪论
利用MATLAB程序递推计
算，结果见表
两算法计算值与精确值比较
MATLAB计算程序： E0=0.09531;%初始值 E(1)=1-10*E0; for n=2:8
a2,…,an为0到9中的任意整数。
m为整数，若使 E(a) x* a 1 10mn 成立,
2
则称a近似 x *有n位有效数字。
例5 设 x* 0.002567 , a 0.00256 102 0.256 则 x* a 0.00005 1 104
2
因为m=-2,所以n=2， 即a有2位有效数字。 例6 设x* =8.00001，则a=8.0000具有5位有效数字。
【注】 近似数的有效数字不但给出了近似值的大小, 而且还指出了它的绝对误差限。
数值分析——绪论
1.3 误差定性分析与避免误差危害
1.3.1 算法的数值稳定性
用一个算法进行计算，初始数据误差（舍入误差造 成）在计算中传播使计算结果误差增长很快，则称 该算法是数值不稳定的，否则是数值稳定的。
例7
计算积分 En
12346+0.6+0.6-12345 =0.12346×105+0.000006×105
+0.000006×105-0.12345×105 =0.12346×105-0.12345×105
=1