数值分析也称计算方法

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