数值算法的稳定性

1/32线性代数一④设x*是准确数,x是x*的近似数，称e=x*-x为近似值x的绝对误差，简称误差。上节课内容回顾——反映是近似值与精确值的绝对差值称***xxxxeer为近似值x的相对误差——反映是近似值与精确值的近似程度通常用百分数来表示，相对误差越小，近似程度越高绝对误差限，相对误差限2/32线性代数一④rrxxxxee|**||*|||则称r为近似值x的相对误差限。|e|=|x*-x|，称为近似值x的绝对误差限，简称误差限或精度3/32线性代数一④如果|e|=|x*-x|1/210m-n称近似数x准确到小数点后第n位,从这小数点后第n位数字直到最左边非零数字之间的所有数字都称为有效数字.有效数字有效数字越多,误差越小,计算结果越精确.近似数x=±0.a1a2…an×10m4/32线性代数一④相对误差与有效数字的关系如下：定理1.1设近似数x=±0.a1a2…an×10m有n位有效数字,则其相对误差限为,102111nra定理1.2设近似数x=±0.a1a2…an×10m的相对误差限为,10)1(2111nra则它至少有n位有效数字。5/32线性代数一④§1.4算法的数值稳定性(数值计算中值得注意的问题)一个算法如果输入数据有误差，而在计算过程中舍入误差不增长，则称此算法是数值稳定的，否则称此算法为不稳定的。换句话说：若误差传播是可控制的,则称此算法是数值稳定的,否则称此算法为不稳定的。6/32线性代数一④见教材第2、10—11页例：计算105dxxxInn(1)In0;(2)In单调递减;0lim)3(nnI)1(151161)4(nnInnIn有以下性质在该例中,用上述公式计算积分的值,I0=ln6-ln5≈0.182322的舍入误差在计算过程迅速传播，每次扩大5倍，致使I12=-0.32902110×10-2严重失真，所以这一公式是不稳定的。,...)2,1(151nnIInn有递推公式7/32线性代数一④nInnInnInnIn10.088392260.0243239110.017324716-10.156920.05803870.018810912-0.003290221750.843330.043138780.021237813-0.093374218-254016140.034306390.017056614-0.395442191270.8650.0284686100.0147169152.0438820-6354.23,...)2,1(151nnIInn舍入误差在计算过程迅速传播，每次扩大5倍.所以此算法是不稳定的。8/32线性代数一④)1,2,...,1,(51511mmnnIInn然后取充分大的m对应的Im的一个估计值为计算初值，再逐步用上式算出Im-1，Im-2，...，I1。用上式计算Im可使计算的误差减少5倍，因而它对应的算法是数值稳定的算法。而将公式变为00873061.022151216120I)1(151161)4(nnInn由:可取自n=20计算到n=1,2)1(51)1(61mmIm,...)2,1(151nnIInn9/32线性代数一④nInnInnInnIn200.00873016150.0105205100.015367650.0284684190.00825397140.011229290.016926540.0343063180.00887552130.012039980.093374230.0431387170.00933601120.012976670.021232620.0580389160.00989750110.014071360.023325010.0883922最后得:I0=0.182322与我们开始计算的I0≈0.182322是一样的)1,2,...,1,(51511mmnnIInn该公式给出的算法就是稳定的下面通过例子给出算法数值稳定的几个原则：10/32线性代数一④一、防止相近的两数相减(会耗失许多有效数字,可以用数学公式化简后再做)xx1xxxx111化成例2:当x较大时，计算0.041只有两位有效数字,有效数字的耗失,说明准确度减小,因此,在计算时需要加工计算公式,以免这种情况发生.控制误差传播的几个原则例1:各有五位有效数字的两个数23.034与22.993相减.23.034-22.993=0.04111/32线性代数一④11759760解:11759760例3:用四位浮点数计算结果只有一位有效数字,有效数字大量损失,造成相对误差扩大。这是由两个比较接近的数相减造成的。结果仍然有四位有效数字。这说明了算法设计的重要性。76017591522100.2100.1316100.1318760175915100.1734100.576817607591760175916注：数值计算中要避免有效数字减少。12/32线性代数一④二、防止大数吃小数当两个绝对值相差很大的数进行加法或减法运算时,绝对值小的数有可能被绝对值大的数吃掉从而引起计算结果很不可靠.例4:求一元二次方程x2-(109+4)x+4×109=0的实数根.采用因式分解法,很容易得到两个根为x1=109,x2=4.如采用字长为8位的计算机来计算,求得根为x1=109，x2=0.(怎样计算可得较好的结果?)两者结果不同,因为计算机计算时做加减法要“对阶”,“对阶”的结果使大数吃掉了小数.产生了误差.为了避免由于上述原因引起的计算结果严重失真,可以根据一些具体情况,有时需要把某些算式改写成另一种等价的形式.13/32线性代数一④四、要控制舍入误差的累积和传播见教材第2、10—11页分母接近0，如何改进？如：|x|1时xxsincos1xxcos1sin三、防止接近零的数做除数分母接近零的数会产生溢出错误,因而产生较大的误差,此时可以用数学公式化简后再做.例5计算105dxxxInn(1)In0;(2)In单调递减;0lim)3(nnI)1(5161)4(1nnInnIn有以下性质也可以用等价无穷小替换14/32线性代数一④在该例中,用上述公式计算积分的值,I0=ln6-ln5≈0.182322的舍入误差在计算过程迅速传播，每次扩大5倍，致使I12=-0.32902110×10-2严重失真，所以这一公式是不稳定的。,...)2,1(151nnIInn有递推公式)1,2,...,1,(51511mmnnIInn然后取充分大的m对应的Im的一个估计值为计算初值，再逐步用上式算出Im-1，Im-2，...，I1。用上式计算Im可使计算的误差减少5倍，因而它对应的算法是数值稳定的算法。而将公式变为15/32线性代数一④在该例中,用上述公式计算积分的值,I0=ln6-ln5≈0.182322的舍入误差在计算过程迅速传播，每次扩大5倍，致使I12=-0.32902110×10-2严重失真，所以这一公式是不稳定的。,...)2,1(151nnIInn有递推公式nInnInnInnIn10.088392260.0243239110.017324716-10.156920.05803870.018810912-0.003290221750.843330.043138780.021237813-0.093374218-254016140.034306390.017056614-0.395442191270.8650.0284686100.0147169152.0438820-6354.2316/32线性代数一④)1,2,...,1,(51511mmnnIInn然后取充分大的m对应的Im的一个估计值为计算初值，再逐步用上式算出Im-1，Im-2，...，I1。用上式计算Im可使计算的误差减少5倍，因而它对应的算法是数值稳定的算法。而将公式变为00873061.022151216120I)1(5161)4(1nnInn由:可取自n=20计算到n=1,)1(2)1(51)1(61mmmIm17/32线性代数一④nInnInnInnIn200.00873016150.0105205100.015367650.0284684190.00825397140.011229290.016926540.0343063180.00887552130.012039980.093374230.0431387170.00933601120.012976670.021232620.0580389160.00989750110.014071360.023325010.0883922最后得:I0=0.182322与我们开始计算的I0≈0.182322是一样的)1,2,...,1,(51511mmnnIInn该公式给出的算法就是稳定的18/32线性代数一④五、简化计算步骤,减小运算次数,避免误差积累例7：设A、B、C、D分别是1020、2050、501、1100的矩阵，试按不同的算法求矩阵乘积E=ABCD.解：由矩阵乘法的结合律，可有如下算法1.E=((AB)C)D.计算量N1=11500flop2.E=A(B(CD)).计算量N2=125000flop3.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=9•9•9•……•9但若写成若逐个相乘要用254次乘法,19/32线性代数一④矩阵乘积AB的计算量分析a11a12a13…a1na21a22a23…a2n......…...am1am2amm-1…amnb11b12b13…b1sb21b22b23…b2s......…...bn1bn2bnn-1…bns=[cij]ms因为cij=aikbkj计算每一个cij的计算量为n所以上面AmnBns的计算量为N=mns计算量：一个算法所需的乘除运算总次数，单位是flop.计算量也是衡量一个算法好坏的重要标准。20/32线性代数一④作业：P139、12实验内容实验1.1、实验1.2（P1312）作业：P12—131、3、4、6、821/32线性代数一④实验1.1（数值微分精度与步长的关系）实验目的：数值计算中误差是不可避免的，要求通过本实验初步认识数值分析中两个重要概念：截断误差和舍入误差，并认真体会误差对计算结果的影响。问题提出：设一元函数f:R→R，则f在x0的导数定义为:hxfhxfxfh)()(lim)('0000实验内容:根据不同的步长h可设计两种算法,计算f在x0处的导数22/32线性代数一④计算一阶导数的算法:hxfhxfxf)()()('000hhxfhxfxf2)()()('000请给出两个计算高阶导数的近似算法,并完成如下工作:(1)(2)1、对同样的h，比较(1)式和(2)式的计算结果；2、针对计算高阶导数的算法，比较h取不同值时(1)式和(2)式的计算结果；实验要求：选择有代表性的函数f(x)（最好多选择几个），利用Matlab提供的绘图工具画出该函数在某区间的导数曲线f(s)(x)，再将数值计算的结果用Matlab画出来，认真思考实验的结果。23/32线性代数一④functiony=dsh1(fu,x,h)y=(feval(fu,x+h)-feval(fu,x))/h;y函数文件(扩展名为.m),可以被其命令调用格式:function输出参数=函数名(输入参数)