计算方法实例

第1讲计算方法（数值分析）计算方法研究者的鼻祖应该是泰勒等，计算方法的形成与发展应该是与计算机的产生与发展结伴而形的。记得有部电影《横空出世》中描写科技人员为了得到制造原子弹的数据，有成百上千的人在用古老的计算工具一一算盘，昼夜计算。由于计算机只能表示离散的数（且只能表示有限位、有限个数），因此，使用计算机进行计算，通常会产生误差，得到的结果有些是“可接受”的，有些是“不可接受”的。于是，通常要研究以下两个内容：所求结果的近似值，计算过程产生的误差；得到结果所用的计算方法是否可行、快速，节省计算机资源.其内容大致可分为三个部分：数值逼近；数值代数；微分方程的数值求解第一章计算方法的一般概念§1.1引言通过研究下面这个例题给出所要讨论的几个方面的内容。TOC\o"1-5"\h\zEx.1.1.1建立/广』1：5必的递推公式并求当n=0,1,2,…时,/“的值.解I=|1xn+5对T—5对-1dx=j1"1dx—5』1X"Tdx=1-5I,n0x+5 o ox+5n n-1而10=j1—dx=ln6-ln5.I=ln6-ln5(1.1)0(1.1)< 1 _I=—-5I,n=1,2，…nn n-1由补Ex1_1_1(迭代次数)的结果知：(1)用Maple，以系统默认的精度计算10-120时：evalf(Ex1_1_1(0)); ^182321557 \# 系统默认取十位evalf(Ex1_1_1(10)); .02evalf(Ex1_1_1(20)); 100000.(2)用Maple，以较高精度(取20位)，计算10-120有：evalf(Ex1_1_1(0),20); .1823215567939546262evalf(Ex1_1_1(10),20); .015367550448evalf(Ex1_1_1(20),20); .00800TOC\o"1-5"\h\zevalf(Ex1_1_1(26),20); .1evalf(Ex1_1_1(27),20); 0.evalf(Ex1_1_1(28),20); 0.(3)用Maple,以更高精度(取30位),有：evalf(Ex1_1_1(26),30); .00621001212evalf(Ex1_1_1(26),50); .0062100121272123878307412230028即不同精度要求，所得结果相差很大。为什么会这样呢？哪个结果是120的近似值呢？这可由式(1.1)本身的理论研究给出答案。事实上，按理论推导有/1=-5/0+1一一. .. /V .. . . ... ..设计算机计算/°时的结果为/°,则代入上式经计算机计算有/=-5/+1

10记/0-1记/0-10=S，则有同理/-/=-5/-/=-5£

11 00/广。=一5(、1 -1)=(FWn-2-/n-2)= =(一5"(/0F=(-5)n£.可见,误差/-7的大小既与/-I的误差£有关,也与迭代次数n有关.尽管£可以取得非常nn 0 01小，但误差的传播逐步扩大，结果并不可靠.当£<<—时,才足够精确.5n与书上的思路类似，也可以用下面的关系研究积分的值的近似值：0</二』1—dx〈L』1xndx= T0(nT3).n0x+5 50 5(n+1)有关计算程序参后面的补1.1.1。Ex.1.1.2求方程工2+(a+P)x+109=0的根，其中a=-109,。=-1.解由求根公式得fx1,2-b±\b2-4x1,22a其中a=1,-b=-a-P=109+1=0.1x1010+0.000,000,000,1x1010,c=109=0.1x1010

b2-4ac=0.01x1020-4x1x0.1x1010=0.01x1020-0.000,000,000,04x1020.在只能取到小数点后面第8位的计算机上对阶计算后,将得到-b=0.1x1010+0x1010=0.1x1010，--a吃掉了。-—加减运算时大吃小b2-4ac=0.01x1020-0=0.01x1020--b2吃掉了4ac.这样得到的根将是0.1x1010+0.1x101。 “ 0.1x1010-0.1x1010八=109, x2= 21 =0,这样得到的根将是0.1x1010+0.1x1010小 0.1x1010-0.1x1010八x1= 2_1 =109,x2= 2~1 =0,而这相当于P,c在方程中没有出现,即所求的是方程x2+以x=0的根.这虽然不正确.为了避免以上情况并考虑到分子部分可能出现三个相近的数相减而导致有效数字的严重损失的不利情况,在计算机上求ax2+bx+c=0的根时，是按下列步骤进行的(退化情况a=0或b=0另行考虑)-b-sgn(b)\b2-4ac2ax2a1cx= 2ax1就算b2吃掉了4ac,第一式分子仍是两个同号数的和,在求x］时仍用了c;或4ac吃掉了b2,第一式分子第一项还含有b.这样一般就不再出现本例直接用求根公式计算过程中出现的情况.实验：法1——用求根函数；法2——按上式编程Ex.1.1.3表示p(x)=axn+axn-1+ +a,并求〃=10,a=^—,x=1时的值.0 1 n kk2解通常计算勤⑴要用n+(n-1)+...+2+1=也口次乘法,n次加法而采用我国13世2纪的数学家秦九韶所首创的所谓秦九韶算法p(x)=(…((ax+a)x+a)x+ +a)x+an 0 1 2 n-1 n则只需要n次乘法与n次加法.程序算法：b=a0b^=b^x+a^,k=1,2,—,np(x)=bn n实验程序见后面：实验程序一Ex.1.13__Maple:Ex.1.1.4计算方的近似值，精确到10-5.R(二)= ——，日位于0与1之间n2 (n+1)!2n+1 2解,/ex=2土=X七+R(x),R(x)=——―xn+1,&位于0与工之间.R(二)= ——，日位于0与1之间n2 (n+1)!2n+1 2...=e2=2^^+Rk=0k!.2k n1欲使R(1)<10-5,只要 一一—<——=—1—<-^,故取〃=6.n2 (n+1)!2n+1 (n+1)!-2n+1 (n+1)!-2n 100000(7!*2人6=322560>100000)为保证前7项计算结果满足精度0.00001的要求，每项都要精确到第7位，经四舍五入累加到第六位，再四舍五入也能保证最终精度.(因为算到小数点后的6位再进行四舍五入可能会产生大于10-5的误差.).•一v6 1「. \evk广kk=0练习题：一、给出前21个结果(以任何方式算出)补1.1.1：用Maple计算：(1)对于浮点数，Maple将始终按照输入时的精度对它进行处理，而且会将计算结果同算式中小数点后有效位数最多的一个因子保持精度一致；(2)对于非浮点运算，包括含e,n,sin(x),ln(x)等,Maple会作精确计算(因为它会作符号运算)，必要时用evalf(%,n)给出输出结果；(3)不加位数n给出15位。restart;Ex1_1_1:=proc(n::numeric)ifn<1thenf(x):=xAn/(x+5);int(f(x),x=0..1);else1/n-5*Ex1_1_1(n-1);endif;end;Ex1_1_1:=proc(n::numeric)ifn<1thenf(x):=xAn/(x+5);int(f(x),x=0..1)else1/n-5xEx1_1_1(n-1)endifendprocevalf(Ex1_1_1(0));.182321558evalf(Ex1_1_1(1));.088392212evalf(Ex1_1_1(2));.05803894evalf(Ex1_1_1(3));.0431386evalf(Ex1_1_1(4));.034307evalf(Ex1_1_1(5));.028466evalf(Ex1_1_1(6));.02434evalf(Ex1_1_1(7));.0212evalf(Ex1_1_1(8));.0191evalf(Ex1_1_1(9));.016evalf(Ex1_1_1(10));.02evalf(Ex1_1_1(11));-.02evalf(Ex1_1_1(12));.3evalf(Ex1_1_1(13));0.evalf(Ex1_1_1(14));5.evalf(Ex1_1_1(15));-20.evalf(Ex1_1_1(16));200.evalf(Ex1_1_1(17));-1000.evalf(Ex1_1_1(18));3000.evalf(Ex1_1_1(19));-20000.evalf(Ex1_1_1(20));100000.结果分析：1）计算I20，用系统默认位数给出结果为134922987458706126014977597520>Ex1_1_1(20); 134922987458706126014977597520卜95367431640625ln(2)+95367431640625ln(3)-95367431640625ln(5)134922987458706126014977597520取精确结果134922987458706126014977597520>1349229874587061260149/77597520;>evalf(%); .17387538601014取精确结果后三项单独计算得：>95367431640625*ln(2)+95367431640625*ln(3)-95367431640625*ln(5);95367431640625ln(2)+95367431640625ln(3)-95367431640625ln(5)>evalf(%);>evalf(%);.1738753871014用后一个结果减去前一个结果时，计算机是先对阶计算小数的差，得0.000,000,001，再乘1014因为这个数有9位小数，所以乘后的结果就为10000.2)如果是取20位，分别得.173875386041597883561014.17387538604159796361014后者小数减去前者小数后再将小数点右移14位得0.79636-0.788356=.008004虽然这一结果与结果evalf(Ex1_1_1(20),20);.00800有所差异，但已属正常又：forifrom0to20doevalf(Ex1_1_1(i),20);od;.1823215567939546262.0883922160302268689.058038919848865656.04313873408900506.0343063295549747.0284683522251264.024324905541035.02123261515197.01883692424015.0169264899103.015367550448.014071338668

.01297664000.0120398770.0112291867.010520734.00989634.0093419.0088462.008400.00800实验程序——Ex.1.13__Maple:法1：restart;>with(linalg):Ex1_1_3:=proc(a,x)>localn,b,k;n:=linalg[vectdim](a);>b:=a[1];forkfrom2tondo>b:=b*x+a[k];od;print(p(x)=b);>print(p(x)=simplify(%));end：c:=array(1..5);c:=array(1..5,[])Ex1_1_3(c,x);p(x)=(((cx+c)x+c)x+c)x+c

1 2 3 4 5p(x)=cx4+x3c+x2c+xc+c1 2 3 4 5当输入多项式系数序列及X后就可具体计算出值来.>c:=array(1..5);c:=array(1..5,[])Ex1_1_3(c,x);TOC\o"1-5"\h\zp(X)=cx4+x3c+x2c+xc+c

1 2 3 4 5>c:={seq(1/k人2,k=1・・10)};-11 1111111、c:={1 },,, ,,, ,,64811009162536494d:=convert(c