計算方法實例

第1講計算方法（數(shù)值分析）計算方法研究者的鼻祖應該是泰勒等，計算方法的形成與發(fā)展應該是與計算機的產生與發(fā)展結伴而形的。記得有部電影《橫空出世》中描寫科技人員為了得到制造原子彈的數(shù)據(jù)，有成百上千的人在用古老的計算工具一一算盤，晝夜計算。由于計算機只能表示離散的數(shù)（且只能表示有限位、有限個數(shù)），因此，使用計算機進行計算，通常會產生誤差，得到的結果有些是“可接受”的，有些是“不可接受”的。于是，通常要研究以下兩個內容：所求結果的近似值，計算過程產生的誤差；得到結果所用的計算方法是否可行、快速，節(jié)省計算機資源.其內容大致可分為三個部分：數(shù)值逼近；數(shù)值代數(shù)；微分方程的數(shù)值求解第一章計算方法的一般概念§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(迭代次數(shù))的結果知：(1)用Maple，以系統(tǒng)默認的精度計算10-120時：evalf(Ex1_1_1(0)); ^182321557 \# 系統(tǒng)默認取十位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 .. . . ... ..設計算機計算/°時的結果為/°,則代入上式經(jīng)計算機計算有/=-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的誤差￡有關,也與迭代次數(shù)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.在只能取到小數(shù)點后面第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在方程中沒有出現(xiàn),即所求的是方程x2+以x=0的根.這雖然不正確.為了避免以上情況并考慮到分子部分可能出現(xiàn)三個相近的數(shù)相減而導致有效數(shù)字的嚴重損失的不利情況,在計算機上求ax2+bx+c=0的根時，是按下列步驟進行的(退化情況a=0或b=0另行考慮)-b-sgn(b)\b2-4ac2ax2a1cx= 2ax1就算b2吃掉了4ac,第一式分子仍是兩個同號數(shù)的和,在求x］時仍用了c;或4ac吃掉了b2,第一式分子第一項還含有b.這樣一般就不再出現(xiàn)本例直接用求根公式計算過程中出現(xiàn)的情況.實驗：法1——用求根函數(shù)；法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紀的數(shù)學家秦九韶所首創(chuàng)的所謂秦九韶算法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位，經(jīng)四舍五入累加到第六位，再四舍五入也能保證最終精度.(因為算到小數(shù)點后的6位再進行四舍五入可能會產生大于10-5的誤差.).?一v6 1「. \evk廣kk=0練習題：一、給出前21個結果(以任何方式算出)補1.1.1：用Maple計算：(1)對于浮點數(shù)，Maple將始終按照輸入時的精度對它進行處理，而且會將計算結果同算式中小數(shù)點后有效位數(shù)最多的一個因子保持精度一致；(2)對于非浮點運算，包括含e,n,sin(x),ln(x)等,Maple會作精確計算(因為它會作符號運算)，必要時用evalf(%,n)給出輸出結果；(3)不加位數(shù)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，用系統(tǒng)默認位數(shù)給出結果為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用后一個結果減去前一個結果時，計算機是先對階計算小數(shù)的差，得0.000,000,001，再乘1014因為這個數(shù)有9位小數(shù)，所以乘后的結果就為10000.2)如果是取20位，分別得.173875386041597883561014.17387538604159796361014后者小數(shù)減去前者小數(shù)后再將小數(shù)點右移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當輸入多項式系數(shù)序列及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