第一章数值计算方法与误差分析.

数值计算方法：研究怎样利用计算工具求出数学问题的数值解，并对算法的收敛性、稳定性和误差进行分析计算的全过程。构建一个完整的数值算法，包含着以下环节：1.提出数值问题（即对对象建立数学模型）2.构思处理数值问题的基本思想（即提出理论）3.列出计算公式4.设计程序框图5.编制源程序并调试6.做出算法的误差分析从工程实际中抽象出来的数学问题往往很复杂，典型的有：)(minxfXx1、数据点的插值2、曲线拟和3、复杂函数的微积分运算4、非线性方程f(x)=0的根的求解5、当n很大时线性方程组AX=B的求解6、常微分方程的求解参考书籍的几种名称：1、数值分析2、数值计算原理3、计算方法4、算法设计5、计算机数值计算方法与程序设计数值计算中的误差1、误差的种类和来源①模型误差②观测误差③截断误差④舍入误差2、误差的有关概念：①绝对误差：xxx)(真值近似值②绝对误差限：1*)(xxx③相对误差：xxxxxxxxxr)()(④相对误差限：2)(xr相对误差限的概念在实际计算中,由于ε(x)与x都不能准确地求得,因此相对误差εr(x)也不可能准确地得到,我们只能估计它的大小范围。即指定一个适当小的正数η，使|εr(x)|=|ε(x)|/|x*|≤η称η为近似值x*的相对误差限。当|εr(x)|较小时，可以用下式来计算η：η=ξ/|x*|有效数字•为了既能表示近似数的大小，又能表示近似数的精确程度，我们下面介绍有效数字的概念（注意：有效数字既能表示近似数的大小，又能表示近似数的精确程度）。半个单位的概念•我们知道，当x有很多位数字时，常常按照“四舍五入”原则取前几位数字作为x的近似值x*。•例1设x=π=3.1415926…取x1*=3作为π的近似值,则|ε1(x)|=0.1415…≤0.5×100；取x2*=3.14，则|ε2(x)|=0.00159…≤0.5×10-2；取x3*=3.1416,则|ε3(x)|=0.0000074…≤0.5×10-4。它们的误差都不超过末位数字的半个单位。有效数字的概念•定义3若近似值x*的绝对误差限是某一位上的半个单位，该位到x*的第一位非零数字一共有n位，则称近似值x*有n位有效数字，或说x*精确到该位。•准确数本身有无穷多位有效数字，即从第一位非零数字以后的所有数字都是有效数字。有效数字举例•如例1中的x*1，x*2,x*3,分别有1,3,5位有效数字。•实际上，用四舍五入法取准确值x的前n位（不包括第一位非零数字前面的零）作为它的近似值x*时，x*有n位有效数字。•例2设x=4.26972，则按四舍五入法，取2位，x1*=4.3有效数字为2位,取3位，x2*=4.27,有效数字为3位，取4位，x3*=4.270，有效数字为4位。特别注意近似值后面的零不能随便省去，如例2中4.27和4.270，前者精确到4.27,有效数字为3位，取4位，x3*=4.270,有效数字为4位。可见,它们的近似程度完全不同,与准确值的最大误差也完全不同。有效数字和绝对误差的关系定义3换一种说法就是：设x的近似值x*=±0.a1a2…an…×10p若其绝对误差|ε(x)|=|x－x*|≤0.5×10p-n则称近似数x*具有n位有效数字。这里p为整数，a1,a2,…,an是0到9中的一个数字且a1≠0。例如，若x*=0.23156×10-2是x的具有五位有效数字的近似值，则绝对误差是|x－x*|≤0.5×10-2-5=0.5×10-7定义3或式|ε(x)|=|x－x*|≤0.5×10p-n建立了绝对误差(限)和有效数字之间的关系。由于n越大，10p-n的值越小，所以有效数字位越多，则绝对误差(限)越小。有效数字与相对误差的关系定理1若近似数x*具有n位有效数字，则其相对误差为|εr(x)|≤1/(2×a1)×10-(n-1)其中a1≠0是x*的第一位有效数字。定理1说明有效数字位越多,相对误差(限)越小。定理2形式如x*=±0.a1a2…an…×10p的近似数x*，若其相对误差满足|εr(x)|≤1/[2×(a1+1)]×10-(n-1)则x*至少有n位有效数字。由此可知，有效数字位数可刻画近似数的精确度，相对误差(限)与有效数字的位数有关。有效数字与相对误差关系举例注意从并不能保证x*一定具有n位有效数字。如x=sin29020´=0.4900设其近似值x*=0.484,其相对误差为我们不能由此推出x*有两位有效数字，这是因为x-x*=0.4900-0.484=0.00600.005即可知近似值x*并不具有两位有效数字。实际上,x*只有一位有效数字。|εr(x)|≤1/(2×a1)×10-(n-1))12(110421104210125.0012397.0484.0484.04900.0数值运算中误差的影响要分析数值运算中误差的传播，首先就要估计数值运算中的误差。数值运算的误差估计情况较复杂,通常利用微分来估计误差。二元函数设数学问题的解y与变量x1,x2有关,y=f(x1,x2)。若x1,x2的近似值为x1*,x2*,相应解为y*,则当数据误差较小时解的绝对误差ε(y)=y－y*=f(x1,x2)－f(x1*,x2*)≈dy=∂f(x1,x2)/∂x1*ε(x1)+∂f(x1,x2)/∂x2*ε(x2)解的相对误差εr(y)≈dy/y=Σ∂f(x1,x2)/∂xi*xi/f(x1,x2)*εr(xi)(i=1,2)利用这两式可得到两数和、差、积、商的误差估计。算法的数值稳定性一、算法的数值稳定性概念所谓算法，是指对一些数据按某种规定的顺序进行的运算序列。在实际计算中，对于同一问题我们选用不同的算法,所得结果的精度往往大不相同。这是因为初始数据的误差或计算中的舍入误差在计算过程中的传播，因算法不同而异，于是就产生了算法的数值稳定性问题。一个算法,如果计算结果受误差的影响小，就称这个算法具有较好的数值稳定性。否则，就称这个算法的数值稳定性不好。算法的数值稳定性概念举例例1一元二次方程X2+2pX+q=0的两个根分别是：x1=–p+(p2–q)0.5，x2=–p–(p2–q)0.5当p=–0.5×105,q=1时，方程的两个根取11位有效数字为：x1=99999.999990,x2=0.000010000000001在高精度的计算机(进制β=10,字长t=8,浮点阶码下限L=–50,浮点阶码上限U=50)上直接用上述公式计算的结果为：x1=100000.00,x2=0可见，结果x1很好，而x2很不理想，这说明直接用上述公式计算第二个根是不稳定的，其原因在于在计算x2时造成相近两数相减，从而使有效数字严重损失。请看下面的求解方法。一元二次方程X2+2pX+q=0的求解方法根据根与系数的关系可知x1x2=q=1所以x2=1/x1因此，如果仍用上述方法算出x1，然后用x2=1/x1计算x2，可得x1=100000.00，x2=0.00001000该结果是非常好的。这就说明这种算法有较好的数值稳定性。一般说来，当|p||q|时，用公式x1=–p–sign(p)•(p2–q)0.5,x2=q/x1来求解方程X2+2pX+q=0是数值稳定的。从而可知，算法数值稳定性的讨论甚为重要。二、设计算法的若干原则为防止误差使计算结果失真(失常）现象发生，要选用数值稳定的计算公式，以保证算法的数值稳定性。下面我们给出设计算法的若干原则，并给出改善算法的例子，这些原则有助于鉴别算法的可靠性并防止误差危害的现象产生。（一）要避免相近两数相减下面再举几个例说明改善算法的方法。例２x充分大时1/x–1/(x+1)=1/[x(x+1)](1+x)1/2–x1/2=1/[(1+x)1/2+x1/2]例３对于小的正数εsin(x+ε)–sinx=2cos(x+ε/2)sin(ε/2)(注：sin(x)–sin(y)=2cos[(x+y)/2]sin[(x-y)/2])例４对于绝对值小的x，可利用泰勒级数ex–1=x+x2/2+x3/6+…取前n项来计算。（二）要防止大数“吃掉“小数，注意保护重要数据在数值运算中，参加运算的数有时数量级相差很大，而计算机位数有限，如不注意运算次序就可能出现大数“吃掉”小数的现象，影响计算结果的可靠性。例5在五位浮点十进制计算机上，计算y=54321+0.4+0.3+0.4如果按从左到右的顺序进行加法运算，后三个数都在对阶过程中被当作零，得出含有较大绝对误差的结果y=54321。要避免这种大数“吃掉”小数的现象，可以调整计算顺序，采用先小数后大数的计算次序，即先将0.4,0.3,0.4加起来，然后再加上54321，结果等于54322。一般情况下，若干数相加，采用绝对值较小者先加的算法，结果的相对误差限较小。（三）注意简化计算步骤、减少运算次数、避免误差积累•同一个计算问题，如果能减少运算次数，不但可以提高计算速度，而且能减少误差的积累。简化计算步骤、减少运算次数、避免误差积累的例子例6计算多项式P4(x)=0.0625x4+0.425x3+1.215x2+1.912x+2.1296的值。如果先计算各项然后相加，需做十次乘法和四次加法。如改用下式计算(((0.0625x+0.425)x+1.215)x+1.912)x+2.1296则只需做四次乘法和四次加法。简化计算步骤、减少运算次数、避免误差积累的例子又如计算1/(1*2)+1/(2*3)+…+1/(1000*1001)的值。若一项一项进行计算，不仅计算次数多，而且误差积累也很大。若简化成1－1/1001进行计算，则整个计算只要一次求倒数和一次减法。（四）要避免绝对值小的数作除数由式ε(x1/x2)≈d(x1/x2)≈[x2ε(x1)－x1ε(x2)]/x22,(x2≠0)可知，当除数x2接近于零时，商的绝对误差就可能很大。因此,在数值计算中要尽量避免绝对值小的数作除数,避免的方法是把算式变形或改变计算顺序。例8当x接近于0时(1-cosx)/sinx的分子、分母都接近0，为避免绝对值小的数作除数，可将原式化为(1-cosx)/sinx=sinx/(1+cosx)例9当x很大时，可化x/[(x+1)0.5－x0.5]=x[(x+1)0.5+x0.5]控制误差传播的例子例10计算积分In=∫01xnex-1dx，n=0,1,2,…,9利用分部积分法，可得In=xnex-1|01–∫01ex-1dxn=1–n∫01xn-1ex-1dx=1–nIn-1从而有递推公式I0=∫01ex-1dx=ex-1|01=1-e-1≈0.6321In=1–nIn-1（n=0,1,2,…,9）计算积分In=∫01xnex-1dx的过程如果直接应用递推公式I0=∫01ex-1dx=ex-1|01=1-e-1≈0.6321In=1–nIn-1（n=1,2,…,9）用四位小数计算依次得到：0.6321,0.3679,0.2642,0.2074,0.17040.1480,0.1120,0.2160,-0.7280,7.5520由此看到I8为负值、I91，显然与一切0In1（由于e-1/(n+1)=min(ex-1)∫01xndxIn(0≤x≤1)max(ex-1)∫01xndx=1/(n+1))矛盾。事实上，从I7开始已经连一位有效数字也没有了（I71/8=0.125,而上面算得I7=0.2160）。是什么原因造成这种结果呢？计算积分In=∫01xnex-1dx的误差分析根据I0=∫01ex-1dx=ex-1|01=1-e-1≈0.6321In=1–nIn-1（n=1,2,…,9）计算In,假定I0与准确值的误差为ε(I0)。容易看出误差传递规律是:ε(In)=–nε(In-1)=…=(-1)nn!ε(I0)（如：I4=1-4I3=1-4(1-3I2)=1-4(1-3(1-2I1))=1-4(1-3(1-2(1-I0))=1-4+12-24+4!I0）这说明由初始值I0的舍入而产生的误差在计算过程中绝对值会迅速扩大。