数值算法的稳定性.ppt

1/32 线性代数线性代数一一 设x*是准确数, x是x*的近似数，称e = x* - x 为近似 值x的绝对误差，简称误差。 上节课内容回顾 反映是近似值与精确值的绝对差值 称 为近似值x的相对误差 反映是近似值与精确值的近似程度 通常用百分数来表示，相对误差越小，近似程度越高 绝对误差限，相对误差限 2/32 线性代数线性代数一一 则称r 为近似值x的相对误差限。 |e| = |x* - x| ， 称 为近似值x的绝对误差限，简称误差限或精度 3/32 线性代数线性代数一一 如果|e| = |x* - x| 1/2 10m-n 称近似数x准确到小数点后第n位,从这小数点后第n位 数字直到最左边非零数字之间的所有数字都称为有效 数字. 有效有效数字数字 有效数字越多,误差越小,计算结果越精确. 近似数x=0.a1a2an10m 4/32 线性代数线性代数一一 3.142=0.3142101， 3.1416=0.31416101， 3.141=0.3141101 的近似值如下例： 用四舍五入得到的数都是有效数字. 4位有效数字 5位有效数字 3位有效数字 m=1, n=4 m=1, n=5 m=1, n=3 |e| = |x* - x| 1/2 10m-n 5/32 线性代数线性代数一一 相对误差与有效数字的关系如下： 定理1.1 设近似数x=0.a1a2an10m有n位有效数字,则 其相对误差限为 定理1.2 设近似数x=0.a1a2an10m的相对误差限为 则它至少有n位有效数字。 6/32 线性代数线性代数一一 三、四则运算的误差计算三、四则运算的误差计算 绝对误差：e = x* - x =xdx 相对误差：er =(x*-x)/x* dx/x=dlnx 利用这个关系可以讨论四则运算的误差和函数的误差 例如下列式子说明什么误差结果? d(x+y)=dx+dy dln(xy)=dlnx+dlny dln(xn)=ndlnx 思考: 若则 7/32 线性代数线性代数一一 作业：P1213 1、3、 4、6 四则运算误差限的公式: 8/32 线性代数线性代数一一 近似计算中误差是不可避免的 ，但不能无限 扩大。如何控制误差的传播，是设计算法时必须考 虑的问题。本节将要介绍的内容就是如何控制误差 的传播，给出一般的原则。 解决一个计算问题往往有多种算法，用不同算法 计算的结果其精度往往大不相同。这是因为初始数据 的误差或计算中的舍入误差在计算过程中的传播造成 的，这种传播因算法不同而相异。 9/32 线性代数线性代数一一 1.4 算法的数值稳定性 (数值计算中值得注意的问题) 一个算法如果输入数据有误差，而在计算过 程中舍入误差不增长，则称此算法是数值稳定的 ，否则称此算法为不稳定的。 换句话说：若误差传播是可控制的,则称此算法是 数值稳定的,否则称此算法为不稳定的。 10/32 线性代数线性代数一一 见教材第2、1011页 例：计算 (1)In0; (2)In单调递减 In有以下性质 在该例中,用上述公式计算积分的值,I0=ln6-ln5 0.182322的舍入误差在计算过程迅速传播，每次扩大 5倍，致使I12= -0.3290211010-2 严重失真，所以这一 公式是不稳定的。 有递推公式 11/32 线性代数线性代数一一 nInnInnInnIn 1 0.0883922 6 0.0243239 11 0.0173247 16 -10.1569 2 0.058038 7 0.0188109 12 -0.00329022 17 50.8433 3 0.0431387 8 0.0212378 13 -0.0933742 18 -2540161 4 0.0343063 9 0.0170566 14 -0.395442 19 1270.86 5 0.0284686 10 0.0147169 15 2.04388 20 -6354.23 舍入误差在计算过程迅速传播，每次扩大5倍. 所以此算法是不稳定的。 12/32 线性代数线性代数一一 然后取充分大的m对应的Im的一个估计值为计算初 值，再逐步用上式算出Im-1 ，Im-2 ，.，I1。 用上式计算 Im 可使计算的误差减少5倍，因而它对应 的算法是数值稳定的算法。 而将公式变为 由: 可取 自n=20计算到n=1, 13/32 线性代数线性代数一一 nInnInnInnIn 20 0.00873016 15 0.0105205 10 0.0153676 5 0.0284684 19 0.00825397 14 0.0112292 9 0.0169265 4 0.0343063 18 0.00887552 13 0.0120399 8 0.0933742 3 0.0431387 17 0.00933601 12 0.0129766 7 0.0212326 2 0.0580389 16 0.00989750 11 0.0140713 6 0.0233250 1 0.0883922 最后得: I0=0.182322 与我们开始计算的I00.182322是一样的 该公式给 出的算法 就是稳定 的 下面通过例子给出算法数值稳定的几个原则： 14/32 线性代数线性代数一一 一、防止相近的两数相减 (会耗失许多有效数字,可以用数学公式化简后再做) 例2:当x较大时，计算 0.041只有两位有效数字,有效数字的耗失,说明准确度减 小,因此,在计算时需要加工计算公式,以免这种情况发生. 控制误差传播的几个原则 例1: 各有五位有效数字的两个数23.034与22.993相减. 23.034-22.993=0.041 15/32 线性代数线性代数一一 解: 例3:用四位浮点数计计算 结结果只有一位有效数字,有效数字大量损损失,造成相 对误对误 差扩扩大。这这是由两个比较较接近的数相减造成的。 结结果仍然有四位有效数字。这说这说 明了算法设计设计 的重要性。 注：数值计算中要避免有效数字减少。 16/32 线性代数线性代数一一 二、防止大数吃小数. 当两个绝对值相差很大的数进行加法或减法运算时,绝 对值小的数有可能被绝对值大的数“吃掉“从而引起计算 结果很不可靠. 例4:求一元二次方程x2-(109 +4)x+4109=0 的实数根. 采用因式分解法,很容易得到两个根为x1=109,x2=4.如采用 字长为8位的计算机来计算,求得根为x1=109 ，x2=0.(怎样 计算可得较好的结果?) 两者结果不同,因为计算机计算时做加减法要 “对阶 ”,“对阶”的结果使大数吃掉了小数.产生了误差.为了避免 由于上述原因引起的计算结果严重失真,可以根据一些具 体情况,有时需要把某些算式改写成另一种等价的形式. 17/32 线性代数线性代数一一 四、要控制舍入误差的累积和传播 见教材第2、1011页 分母接近0，如何改进？ 如：|x|0; (2)In单调递减 In有以下性质 也可以用等价无穷替换 18/32 线性代数线性代数一一 在该例中,用上述公式计算积分的值,I0=ln6-ln5 0.182322的舍入误差在计算过程迅速传播，每次扩大 5倍，致使I12= -0.3290211010-2 严重失真，所以这一 公式是不稳定的。 有递推公式 nInnInnInnIn 1 0.0883922 6 0.0243239 11 0.0173247 16 -10.1569 2 0.058038 7 0.0188109 12 -0.00329022 17 50.8433 3 0.0431387 8 0.0212378 13 -0.0933742 18 -2540161 4 0.0343063 9 0.0170566 14 -0.395442 19 1270.86 5 0.0284686 10 0.0147169 15 2.04388 20 -6354.23 19/32 线性代数线性代数一一 然后取充分大的m对应的Im的一个估计值为计算初 值，再逐步用上式算出Im-1 ，Im-2 ，.，I1。 用上式计算 Im 可使计算的误差减少5倍，因而它对应 的算法是数值稳定的算法。 而将公式变为 由: 可取 自n=20计算到n=1, 20/32 线性代数线性代数一一 nInnInnInnIn 20 0.00873016 15 0.0105205 10 0.0153676 5 0.0284684 19 0.00825397 14 0.0112292 9 0.0169265 4 0.0343063 18 0.00887552 13 0.0120399 8 0.0933742 3 0.0431387 17 0.00933601 12 0.0129766 7 0.0212326 2 0.0580389 16 0.00989750 11 0.0140713 6 0.0233250 1 0.0883922 最后得: I0=0.182322 与我们开始计算的I00.182322是一样的 该公式给 出的算法 就是稳定 的 21/32 线性代数线性代数一一 五、简化计算步骤,减小运算次数,避免误差积累 例7：设A、B、C、D分别是1020、 2050、 501、 1100的 矩阵，试按不同的算法求矩阵乘积E=ABCD. 解：由矩阵乘法的结合律，可有如下算法 1. E=(AB)C)D. 计算量N1=11500flop 2. E=A(B(CD). 计算量N2=125000flop 3. E=(A(BC)D. 计算量N3=2200flop 简化计算步骤是提高程序执行速度的关键，它不仅可以节省 时间，还能减少舍入误差。 例6：计算9255的值 9255 = ( 92 )2 )2 )2 ) 2 )2 )2 )2 /9 只需8次乘法和1次除法运算。 9255 = 9 92 94 98 916 932 964 9128 只需做14次乘法运算即可。 9255 = 9999 但若写成 若逐个相乘要用254次乘法, 22/32 线性代数线性代数一一 矩阵乘积AB的计算量分析 a11 a12 a13 a1n a21 a22 a23 a2n . . . am1 am2 amm-1 amn b11 b12 b13 b1s b21 b22 b23 b2s . . . bn1 bn2 bnn-1 bns =cijms 因为 cij=aik bkj 计算每一个cij的计算量为n 所以上面A mn B ns的计算量为N= m n s 计算量：一个算法所需的乘除运算总次数，单位是flop. 计算量也是衡量一个算法好坏的重要标准。 23/32 线性代数线性代数一一 作业：P13 9、12 实验题：实验1.1及 P13 12 作业：P1213 1、3、 4、6、8 24/32 线性代数线性代数一一 实验1.1（数值微分精度与步长的关系） 实验目的：数值计算中误差是不可避免的，要求通过 本实验初步认识数值分析中两个重要概念：截断误差 和舍入误差，并认真体会误差对计算结果的影响。 问题提出：设一元函数f :RR，则f在x0的导数定义为: 实验内容:根据不同的步长h可设计两种算法,计算f在 x0处的导数 25/32 线性代数线性代数一一 计算一阶导数的算法: 请给出两个计算高阶导数的近似算法,并完成如下工作: (1) (2) 1、对同样的h，比较(1)式和(2)式的计算结果； 2、针对计算高阶导数的算法，比较h取不同值时(1) 式和(2)式的计算结果； 实验要求：选择有代表性的函数f (x)（最好多选择几个 ），利用Matlab提供的绘图工具画出该函数在某区间 的导数曲线f (s)(x)，再将数值计算的结果用Matlab 画 出来，认真思考实验的结果。 26/32 线性代数线性代数一一 function y=dsh1(fu,x,h) y=(feval(fu,x+h)-feval(fu,x)/h; y 函数文件(扩展名为.m), 可以被其命令调用 格式: function输出参数=函数名(输入参数) 可以是 多个 函数名 一个 注意: 1、保存的文件名与函数名要一 致，这样才能保证调用成功； 2、调用时的输入输出参数要与 函数中的一样； 3、养成良好的注释习惯，便于 自己或用户调用定义的函数。 注释语句(前面有%) 程序语句 function y=dsh2(fu,x,h) y=(feval(fu,x+h)-feval(fu,x-h)/(2*h); y 调用格式：dsh2(sin,1,0.01) 调用格式：dsh1(sin,1,0.01) 实验报告要上交，不能有相同的报告。 （否则，两个都作废，不予记分） 27/32 线性代数线性代数一一 真解: cos(1)=0.5403 ans = 0.5403 dsh1(sin,1,0.1) y = 0.4974 dsh1(sin,1,0.01) y = 0.5361 dsh1(sin,1,0.001 ) y = 0.5399 dsh1(sin,1,0.0001) y = 0.5403 f(x)=sin(x) -2,2的图形 要求：真解和 近似解做在一 张图上，近似 解的散点图有 红色的*表示, 便于比较. h=0.1 有偏差 28/32 线性代数线性代数一一 29/32 线性代数线性代数一一 dsh2(sin,1,0.1) y = 0.5394 dsh2(sin,1,0.01) y = 0.5403 真解: cos(1)=05403 ans = 0.5403 大家可以找其他 的函数来检验 -2,2的图形 plot(x,y,*r,x,y1,-) y=dsh2(sin,- 2*pi:0.1:2*pi,0.01) ; y1=cos(x) x=-2*pi:0.1:2*pi; h=0.01 几乎没有 偏差 30/32 线性代数线性代数一一 31/32 线性代数线性代数一一 dsh1(sin,1,0.00000001) y = 0.5403 dsh1(sin,1,0.0000000000001) y = 0.5396 步长为10-12 dsh1(sin,1,0.00000000000001) y = 0.5440 dsh1(sin,1,0.000000000000001) y = 0.5551 dsh1(sin,1,0.0000000000000001) y = 0 步长为10-15 32/32 线性代数线性代数一一 截断误差用我们原有的数学思维方式就比较容易 理解的，而舍入误差则是本课程引入的一个新概念。 要真正理解舍入误差，特别是它在计算中的传播及最 终对计算结果的影响，是初步具备科学计算能力的重 要标志。希望大家在完成实验后，认真仔细去体会截 断误差和舍入误差的含义及对计算结果的影响。 图 0 10-710-5 10-5 10-4 10-3 10-2 10-1 10-310-1 误 差 步长h 33/32 线性代数线性代数一一 实验分析：不论采用怎样的算法，计算结果通常都 会有误差。比如算法(1)，由Taylor公式，知： 所以有 利用上式来计算f (x0)，其截断误差为： 所以误差是存在的，并且当步长h越来越小时，(1)的 近似程度也越来越好。 截断误差 34/32 线性代数线性代数一一 类似地可以分析(2)的截断误差为： 截断误差 上述截断误差的分析表明(2)是比(1)更好的算法，因为 对步长h(0,n=1,2,当n=1时，得： 当n2时，由分部积分可得： n=2,3, 另外，还有： 37/32 线性代数线性代数一一 实验内容：由递推关系In=1-nIn-1，可得计算积分 序列In的两种算法： 算法一、直接使用递推公式得： In=1-nIn-1 n=2,3 算法二、利用递推公式变形得： 实验要求：用上述两种算法分别在计算中采用5位、6 位和7位有效数字，请判断哪种算法给出的结果更精确 . 38