A第一章誤差與范數數值分析也稱計算方法，它不僅是一種研究并解決數學問題的數值近似解的方法，而且是在計算機上使用的解數學問題的方法.它的計算對象是那些在理論上有解，而無法用手工計算的數學問題.1.1誤差的來源在運用計算方法解決實際問題的過程中，會出現各種各樣的誤差，必須注重誤差分析.否則，一種合理的計算也可能得出錯誤的結果.例1.1.1用差商f'(a)Mf(a+")—f(a)求f(x)=Inx在x=3處導數的近似值.h取h=0.1和h=0.0001，用手工計算，取五位數字計算；取h=0.1，h=0.0001，h=0.000000000000001和h=0.0000000000000001分別用MATLAB軟件計算，取十五位數字計算；比較以上的運算結果，說明是否|h|越小則計算結果越準確.解根據導數定義，可以用差商求f(x)=lnx在x=3處導數的近似值f'(3)mf'(3)mf(3+h)-f(3)ln(3+h)—ln3

h從理論上講，Ih從理論上講，Ih|越小則計算結果越準確.如果用手工計算，取五位數字計算.當h=0.1時，得1.131f'(3)m4-1.09861.131f'(3)m4-1.0986=0.3280，0.1與導數的精確值f'(3)=1=0.33333比較這項計算還是可取的.但是當h=0.0001時，得f'(3)Mf'(3)M1.0986—1.0986=0，0.0001算出的結果反而毫無價值.如果應用MATLAB軟件計算，取十五位數字計算，結果就完全不同了.在MATLAB工作窗口輸入下面程序>>a=3;h=0.1;y=log(a+h)-log(a);yx=y/h運行后得yx=0.32789822822991將此程序中h改為0.0001，運行后得yx=0.33332777790385后者比前者好.再取h=0.000000000000001，運行后得yx=0.44408920985006不如前者好.取h=0.0000000000000001，運行后得yx=0算出的結果反而毫無價值.

例1.1.2分別求方程組AX=b在下列情況時的解，其中A=""1 i.01)（1）b="2.01j解（1）首先將方程組AX=b化為同解方程X=A-1b，然后在MATLAB工作窗口輸入程序>>b=[2,2]';A=[1,1;1,1.01];X=A\b運行后輸出當b=時，AX運行后輸出當b=時，AX=b的解為X=（2）同理可得，當bAX=b的解為X例1.1.3計算e的近似值.解泰勒級數取X=1，得這是一個無限過程，誤差.如果取有限項1_

這是一個無限過程，誤差.如果取有限項1_

一+3!計算機無法求到精確值.只能在（1.2）（1.2）取有限項時計算，再估計s(1)=1+1+—+—+—+...+n 2! 3! 4!(0<9<1).作為e的值必然會有誤差，根據泰勒余項定理可知其截斷誤差為(0<9<1).Sn^—(n+1)!如果取(1.2)的前九項，輸入程序>>n=8;s=1;S=1;fork=1:ns=s*k;S=S+1/s,ends,S,R=3/(s*(n+1))或>>S1=1+1+1/2+1/(1*2*3)+1/(1*2*3*4)+1/(1*2*3*4*5)+1/(1*2*3*4*5*6)+1/(1*2*3*4*5*6*7)+1/(1*2*3*4*5*6*7*8),R1=3/(1*2*3*4*5*6*7*8*9)運行后結果R=8.267195767195768e-006R=8.267195767195768e-0062.71827876984127因為截斷誤差為所以e的近似值e a3<一a所以e的近似值e a3<一a8.267(8+1)! 9!1S8⑴=1+1+2!+196x10-6(0<9<1),1+3!+6!1-2.71828.8!1.2誤差和有效數字例1.2.1取2.71828作為e的四舍五入近似值時，求其絕對誤差和相對誤差.解在MATLAB工作窗口輸入程序>>juewu=exp(1)-2.71828運行后輸出結果為juewu=1.828459045505326e-006例1.2.2計算』：Sin*dx的近似值，并確定其絕對誤差和相對誤差.解因為被積函數的原函數不是初等函數，故用泰勒級數求之.x(一3<x<3), (1.5)TOC\o"1-5"\h\zsinx x2 x4(一3<x<3), (1.5)=1- + - + +x 3! 5! 7! 9!這是一個無限過程，計算機無法求到精確值.可用(這是一個無限過程，計算機無法求到精確值.可用(1.5)的前四項］一M+x4x±代7!"』：土dx提:(1-土+"』：土dx提:(1-土+0x 0 3!(-)5q 55!x6)dx=壬一一+5! 7! 2 33!根據泰勒余項定理和交錯級數收斂性的判別定理，得到絕對誤差兀n()9R=Iy-yk-^-WU，

9-9!在MATLAB命令窗口輸入計算程序如下：symsxf=1-xA2/(1*2*3)+xA4/(1*2*3*4*5)-xA6/(1*2*3*4*5*6*7)y=int(f,x,0,pi/2),y1=double(y)y11=pi/2-(pi/2)A3/(3*3*2)+(pi/2)A5/(5*5*4*3*2)-(pi/2)A7/(7*7*6*5*4*3*2)infd=double(inf),inf=int(sin(x)/x,x,0,pi/2)WU=(pi/2)A9/(9*9*8*7*6*5*4*3*2),R=infd-y11因為運行后輸出結果為：y=1.37076216815449，y-1.37074466418938，R=1.75039651049147e-005，WU-1.782679830970664e-005<10-4.所以，y的絕對誤差為￡=10-4，故y=J：-S!宜dx*1.3707.y的相對誤差為0x￡ 10-4￡^=——=E0—<0.0073%.1.3誤差估計的基本方法下面的例子揭示了初始誤差的穩定性傳播和不穩定性傳播.例1.3.4設計三種算法求方程2x2+x-15=0在(2,3)的一個正根x*的近似值，并研究每種算法的誤差傳播情況.解為解已知方程，我們可以設計如下三種算法，然后將計算結果與此方程的精確解X*=2.5比較，觀察誤差的傳播.算法1將已知方程化為同解方程x=15-2x2.取初值x0=2，按迭代公式\+1=15-2X2依次計算X],X2,…,Xn,…，結果列入表1-3中.算法2將已知方程化為同解方程x=;%.取初值x0=2，按迭代公式15x= k+1 2x+1依次計算X,X,...,X,...，結果列入表1-3中.1 2 n算法3將已知方程化為同解方程x=x-之X2+X-15.取初值x0=2，按迭代公式為2x2+x-15X+1=Xk- "4x」1依次計算x,x,...,x,...，結果列入表1-3中.12 n我們為這種三算法的計算編寫兩套MATLAB程序如下：(1)MATLAB主程序function[k,juecha,xiangcha,xk]=liti112(x0,x1,limax)%輸入的量--x0是初值,limax是迭代次數和精確值x;%輸出的量--每次迭代次數k和迭代值xk,% --每次迭代的絕對誤差juecha和相對誤差xiangcha，x(1)=x0;fori=1:limaxx(i+1)=fl(x(i));%程序中調用的fl.mjuecha=abs(x(i)-x1);xiangcha=juecha/(abs(x(i))+eps);xk=x(i);k=i-1;[(i-1),juecha,xiangcha,xk]end(2)MATLAB調用函數程序及其計算結果算法2的MATLAB調用函數程序functiony1=fl(x)y1=15/(2*x+1);將MATLAB主程序和調用函數程序分別命名liti112.m和fl.m，分別保存為M文件，然后在MATLAB工作窗口輸入命令>>[k,juecha,xiangcha,xk]=liti112(2,2.5,100)運行后輸出計算結果列入表1-3和表1-4中.將算法2的MATLAB調用函數程序的函數分別用y1=15-2*xA2和y1=x-(2*xA2+x-15)/(4*x+1)代替，得到算法1和算法3的調用函數程序，將其保存，運行后將三種算法的前8個迭代值x『x2,…，x8列在一起(見表1-3)，進行比較.將三種算法的x1,x2,…，x8對應的絕對誤差和相對誤差的值列在一起(見表1—4)，進行比較.表1-3例1.3.4中三種算法的計算結果-算法迭代次數算法1的迭代結果算法2的迭代結果算法3的迭代結果022.000000002.00000000

第一章誤差與范數173.000000002.555555562-832.142857142.500550063-137632.837837842.500000064-3788403232.246963562.500000005-2.8704x10172.246963562.500000006-1.6478x10352.321774842.500000007-5.4307x10702.657901652.50000000:99:-Inf:2.50000001:2.50000000由表1-3和表1-4可以看出，算法1的初始絕對誤差和相對誤差分別是4.50000000和0.64285714,它們是三種算法中最大者，由算法1得到的迭代序列b.｝對應的絕對誤差是不穩定的，且當nrs時，按指數無界增長，相對誤差趨近于不定值，x.T-3，即ix.｝不收斂；算法2的初始絕對誤差和相對誤差分別是0.50000000和0.16666667,它們是三種算法中居于中間，由算法2得到的迭代序列｛x.｝對應的絕對誤差和相對誤差是穩定的，且當nr3時，按指數緩慢遞減且趨近于0，k｛xk｝緩慢地擺動式收斂于精確解2.5（參見算法2的運行結果）；算法3的初始絕對誤差和相對誤差分別是0.05555560和0.02173913,它們是三種算法中最小者，由算法3得到的迭代序列｛x.｝對應的絕對誤差和相對誤差是穩定的，且當nrs時，按指數快速遞減且趨近于0，｛x.｝快速地遞減且收斂于精確解2.5. *表1-4例1.3.4中三種算法計算結果的誤差算法迭代次數算法1的i沃差算法2的誤差算法3的誤差絕對誤差相對誤差絕對誤差相對誤差絕對誤差相對誤差01234567:990.50000000000000000013765.500003788403262.87039981x10171.64783901x10355.43074680x1070:Inf0.250000000.642857141.030120480.000100021.00000001111:NaN0.500000000.500000000.357142860.337837840.253036440.230287040.178225160.15790165:0.000000010.250000000.166666670.1666666700.119047620.112612620.084345480.076762470.05940839:0.000000000.500000000.055555600.000550100.000000060.00000000000:00.250000000.021739130.000219970.000000020.00000000000:01.4數值計算中應注意的問題從例1.3.4可以看出，一個問題的解決，往往要經過多次運算.每一步運算都可能產生誤差，在反復多次計算的過程中，必然產生誤差的傳播和積累.顯然，當誤差積累偏大時，會使計算結果失真.因此，在每一步計算中，都應該防止產生誤差升級的現象.例1.4.1求數x=715x（5+8-19-1）的近似值.解（1）直接用MATLAB命令>>x=(7A15)*(sqrt(1+8A(-19))-1)運行后輸出結果x=0問題出現在兩個相近的數V'1+8*與1相減時，計算機運行程序>>sqrt(1+8A(-19))-1運行后輸出結果ans=0由于計算機硬件只支持有限位機器數的運算，因此在計算中可能引入和傳播舍入誤差.因為有效數字的嚴重損失，導致輸出司+8-19-1的結果為0，計算機不能再與數715繼續進行真實的計算，所以，最后輸出的結果與X的精確值不符.(2)如果化為: 715x8-19x=715x(K1+8-19—1)= —，<1+8-19+1再用MATLAB命令>>x=(7A15)*((8A(-19))/(sqrt(1+8A(-19))+1))運行后輸出結果x=1.6471e-005這是因為司+8-19-1化為 一后，計算機運行程序V1+8-19+1>>x=(8A(-19))/(sqrt(1+8A(-19))+1)運行后的結果為x=3.4694e-018由于有效數字的損失甚少，所以運算的結果3.4694x10-18再與715繼續計算，最后輸出的結果與X的精確值相差無幾.例1.4.2求數y=ln(30-2302—1)的近似值.解(1)直接用MATLAB程序>>x=30;x1=sqrt(xA2-1)運行后輸出結果x1=29.9833輸入MATLAB程序>>x=30;x1=29.9833;y=log(x-x1)運行后輸出結果y=-4.0923(2)因為ln(30-202-1)中的x=30很大，如果采用倒數變換法z=X-\X2-1= j，即1 : ln(30-*302—1)=ln ； =-ln(30+V900-1).輸入MATLAB程序>>x=30;y=-log(x+sqrt(xA2-1))運行后輸出結果y=-4.0941(3)輸入MATLAB程序>>x=30;y=log(x-sqrt(xA2-1))運行后輸出結果y=-4.0941可見，（2）計算的近似值比（1）的誤差小.參加計算的數，有時數量級相差很大.如果不注意采取相應的措施，在它們的加減法運算中，絕對值很小的那個數經常會被絕對值較大的那個數“吃掉”，不能發揮其作用，造成計算結果失真.例1.4.4請在16位十進制數值精度計算機上利用軟件MATLAB計算下面的兩個數x*=111111111111111+0.1+0.3和y*=1111111111111111+0.1+0.3將計算結果與準確值比較，解釋計算結果.解在MATLAB工作窗口輸入下面程序>>x=11111111111111+0.1+0.3,y=1111111111111111+0.1+0.3運行后輸出結果x=1.11111111111114e+014，y=1.111111111111111e+015從輸出的結果可以看出，x=x*，而y豐y*.為什么y*僅僅比X*多一位1，而y豐y*呢？這是因為計算機進行運算時，首先要把參加運算的數寫成絕對值小于1而“階碼”相同的數，這一過程稱為數的“對階”.在16位十進制數值精度計算機上利用軟件MATLAB計算這兩個數，把運算的數x*寫成浮點規格化形式為X*=0.1111111111111110X1015+0.0000000000000001X1015+0.0000000000000003X1015，在16位十進制數值精度計算機上，三項的數都表示為小數點后面16位數字的數與1015之積，所以，計算機沒有對數進行截斷，而是按原來的三個數進行計算.因此，計算的結果x=x*.而y*=0.1111111111111111X1016+0.00000000000000001X1016+0.00000000000000003X1016三項的數都表示寫成絕對值小于1而“階碼”都為1016的數以后，第一項的純小數的小數點后面有16位數字.但是，后兩項數的純小數的小數點后面有17位數字，超過了16位十進制數值精度計算機的存儲量，計算機對后兩項的數都進行截斷最后一位，即后兩項的數都是16位機上的零，再進行計算，所以計算結果與實際不符.1.5向量和矩陣的范數在研究方程組近似解的誤差估計和迭代法收斂性的過程中，需要對向量和矩陣的“大小”加以度量，向量范數和矩陣范數正是這種度量指標，在數值分析中起著重要作用.1.5.1向量范數與MATLAB命令在3維空間中，將向量x=（x,x,x）T的長度記作||x||，通常用||x||=（蕓x2）：來1 2 3 ii=1定義，向量范數的MATLAB命令及功能列入表1-5：表1-5向量范數的MATLAB命令及功能MATLAB命令等價的命令功能Xp=norm(X,P)sum(abs(X).AP)A(1/P)輸入X為向量，輸出為X的p范數1<p<3

第一章誤差與范數X2=norm(X)或X2=norm(X,2)sum(abs(X).A2)A(1/2)輸入X為向量，輸出為X的2范數X1=norm(X,1)sum(abs(X))輸入X為向量，輸出為X的1范數Xw=norm(X,inf)max(a