實驗二數值誤差與差商一、目的本實驗之主要目的在於讓同學們認識.docVIP

實驗二 數值誤差與差商 一、目的： 本實驗之主要目的在於讓同學們認識利用計算機進行數值運算時，可能遭遇的數值誤差，以警惕同學們對於計算結果需作審視與分析。另外我們也開始引入差商的定義與應用，以作為日後微分方程之數值解的基礎 二、原理： 儘管隨著電腦由16位元、32位元而至64位元，計算機依然只能處理實數中極有限的有理數。以64位元而言，我們對浮點數只能在位數，的範圍內工作，否則即會出現overflow或underflow的情況。因此當我們在應用計算機從事數值計算，必須適當的地免開這個限制。這包括：選擇適當的物理單位，適當的運算程序，以及小心的測試與檢驗。本實驗即在以一些例子說明這方面的問題，以作為同學們往後對數值運算結果的分析基礎。 (一)捨位誤差(Rounding error) ： 對於單精度的浮點數，計算機是以32個位元中的24個位元來儲存有效數字，及相當於十進位數字的7位數，至於雙精度則是利用64個位元中的53個位元來儲存有效數字，相當於十進位數字的16位數。因此在運算的過程中，若遇到位數超過該變數所宣告的最大有效數字之上限時，其尾數之超出部份將被捨去。因此將造成與精確計算間之誤差，這種誤差及稱為捨位誤差(rounding error)。由此結果我們可以定義機器所能表達之最小有效浮點數ε (即所謂的機器常數，machine epsilon)，凡是小於ε之浮點在加法的運算中都被視為零。關於這一點我們由下面兩個例子中可看出。 Program epsilon Real*4 a4, eps4, b4 Real*8 a8, eps8, b8 ! a4=1.23456789123456789123 a8=1.23456789123456789123 write(*,*) ‘a4=’,a4 write(*,*) ‘a8=’,a8 ! eps4=1. do while (1.+eps4.gt.1.) Eps4=eps4/2. enddo Write(*,*) ‘eps8=’,eps8*2. ! end 由上列的程式我們可以計算出所使用機器在單精度或雙精度下之ε。瞭解了這個限制之後，我們可以來看看下面這個程式中所遇的捨位誤差。 Program exp1 Sum=1. X=0.5 Fac=1. Xp=1. Do I=1,20 Fac=fac*1*1. Xp=xp*x Si=xp/fac Sum=sum+si Write(*,*) i,sum,si Enddo End 我們可以觀察到，當n7之後，及數項的貢獻已小於ε，搶而無法對求如有所貢獻。因此我們將無法再逼近精確值。 (二) Range error 如前面所提，對單精度之浮點變數，計算機以8個位元來表非指數，故其能表示之最大指數為，同理若為雙精度則最大指數為。 因此任何浮點數的計算均不允許超過這個指數限制，否計算機會出現overflow ()的錯誤訊息。 (三) 截斷誤差(Truncation error) 所謂截斷誤差及來自於以有限的計算機時間與空間來近似理論上的精確值時，所必須引入的截斷而造成的誤差。例如泰勤(Taylor)展開之嚴格定義乃一函數之無窮項的級數和，然而實際以計算機來運算時，必須終止於某一有限的項數，如此便截斷了高階對理論值的貢獻。例甘，我考慮以下列的程式來計算的泰勤級數展開時面臨的截斷誤差。 Program exp2 Real*4 sum4,ex4,fac4,xp4,err4 Open(10,file=’exp2.dat’,access=’append’) Write(10,’(///)’) Write(10,’(3x,’’x’’,12x,’’final term’’,9x,’’relative error’’)’) Write(*,’(3x,’’x’’,12x,’’final term’’,9x,’’relative error’’)’) Do ix=5,-5,-1 If(ix.eq.0) cycle X=2.*ix Ex=exp(x) Sum=1. Fac=1. Xp=1. I=0 Do while(.true.) I=I+1 Fac=1.fac/I Xp=xp*x Term=xp*fac If(sum.eq.sum+term) then Error=abs((sum-ex)/ex) Write(*,’(1x,f5.1,3x,2e20.8)’) x,term,error Write(10,’(1x,f5.1,3x,2e20.8)’) x,t