計算方法緒論與誤差等

計算方法

NumericalAnalysis能源與動力工程學院劉火星82316418liuhuoxing@課程介紹

緒論

誤差分析

線性方程組的解法

常微分方程的初值問題矩陣的特征值和特征向量

插值和擬合

數值微分和數值積分

非線性方程解法CourseOutline徐翠薇、孫繩武，計算方法引論（第二版），高等教育出版社，2002RichardL.Burden&J.DouglasFaires,NumericalAnalysis(SeventhEdition),高等教育出版社，2001數值分析，李慶祥等編，高等教育出版社，2000主要參考書成績評定方法平時成績：20%期末考試:80%先修課程高等數學線性代數程序設計語言緒論Introduction關于計算方法計算數學的過去和未來算法和收斂數學軟件目次

計算機解決實際問題的步驟建立數學模型選擇數值方法編寫程序上機計算實際問題數學問題提供計算方法程序設計上機計算結果分析實際問題工程問題假設實驗模型數學模型

校驗分析/數值解測量物理定律修改

在計算機上是否根據數學公式編程就能得到正確結果?研究例子:求解線性方程組其準確解為x1=x2=x3=1如把方程組的系數舍入成兩位有效數字它的解為x1=-6.222...x2=38.25…x3=-33.65...Numericalanalysisisthestudyofalgorithmsfortheproblemsofcontinuousmathematics----LloydN.Trefethen“計算方法”就是研究在計算機上解決數學問題的理論和數值方法今天的數值計算方法，無論從形式到內容，還是從工具到效果，已遠非半世紀前VonNeumann、Lax等先驅們所處的環(huán)境和條件了，計算機技術和應用軟件的發(fā)展，讓計算數學展開了雙翼。許多迅速發(fā)展的其他學科和社會進步給計算數學的發(fā)展開拓出更為廣闊的新天地計算方法研究的對象研究數值方法的設計、分析和有關理論基礎與軟件實現。計算方法又稱：計算數學、數值方法、數值分析等計算方法的分枝有最優(yōu)化方法、計算幾何、計算概率統(tǒng)計等計算方法的內容

連續(xù)系統(tǒng)的離散化

離散性方程的數值求解

串行計算方法與并行計算方法的關系：并行計算方法70年代初隨并行計算機的出現而產生，是計算數學中最活躍的新領域。計算數學的發(fā)展與科學工程計算是緊密相聯(lián)的，計算數學的發(fā)展歷史也就是與其他學科結合，利用計算機不斷形成新的理論及數值方法并不斷形成新的學科的歷史，例如：“計算物理”“計算流體力學”計算數學發(fā)展的歷史回顧H.Aiken(1900-1973)哈佛大學博士，因做博士論文涉及到空間電荷傳導問題的計算，1937年提出方案，1939年得到IBM資助，1944年建成投入使用。這是繼電式計算機－MarkI三位計算機設計大師的貢獻J.W.Mauchly(1907-1980)賓夕法尼亞物理博士，因從事天氣預報需要想設計計算機，1942年提出計算機方案，1945年底竣工，這就是世界上第一臺電子計算機－ENIAC機三位計算機設計大師的貢獻J.VonNeumann(1903-1957)普林斯頓高級研究所，1945年在普林斯頓研制成MANIAC機，有力地支持美國氫彈研制，稱為計算機之父三位計算機設計大師的貢獻在研制原子彈和氫彈過程中，許多物

理規(guī)律必須通過計算機上的計算摸清

楚。計算物理、理論物理與實驗物理

相輔相成相互促進共同發(fā)展，形成現

代物理學的三大分支由于核武器研制需要，1950年全球只有15臺，到了1962年9月僅美國就有16187臺計算機計算數學發(fā)展的歷史回顧1983年一個由美國著名數學家拉克斯(P.Lax)為首的不同學科的專家委員會向美國政府提出的報告之中，強調“科學計算是關系到國家安全、經濟發(fā)展和科技進步的關鍵性環(huán)節(jié)，是事關國家命脈的大事?！?984年美國政府大幅度地增加對科學計算經費的支持,新建成五個國家級超級計算中心（分別在普林斯頓大學、圣地亞哥、伊里諾大學、康奈爾大學、匹茲堡），配備當時最高性能的計算機，建立NSF-net新網絡

80年代中期我國將“大規(guī)?？茖W與工程計算”列入國家資助重大項目1987年起美國NSF把“科學與工程計算”、“生物工程”“全局性科學”作為三大優(yōu)先資助的領域由于大存儲的高速計算機的使用已導致了科學和技術方面的兩大突出進展：大量用于設計工作的實驗被數學模型的研究逐步取代，如航天飛機設計、反應堆設計、人工心瓣膜設計等能獲取和存儲大量的數據，并能提取隱秘的信息，如計算機層析X射線攝影，核磁共振等1991年以美國總統(tǒng)倡議的形式提出了“高性能計算與通信計劃”。這是為了保持和提高美國在計算和網絡的所有先進領域中的領導地位而制定的。計劃為期五年（1992－1996），投資的重點是發(fā)展先進的軟件技術與并行算法，關鍵技術是可擴展的大規(guī)模并行計算要求到1996年高性能計算能力提高14倍，達到每秒萬億次浮點運算速度（1012

Teraops/S）。計算機網絡通迅能力提高1百倍，達到每秒109位（Gigabits/S）該計劃中列舉的“挑戰(zhàn)”項目有：磁記錄技術、藥物設計、催化、燃燒、海洋模擬、臭氧洞、空氣污染、高速民用運輸機、數字解剖、蛋白質結構設計、金星成像等1993年初美國總統(tǒng)發(fā)布“發(fā)展信息高速公路”（NII）的總統(tǒng)令1994年4月美國總統(tǒng)發(fā)布“建立國家（地球）空間數據基礎實施”（NSDI）的總統(tǒng)令數值方法和數值軟件過去50年的主要進展Before1940Newton’smethod;Gaussianelimination;Gaussquadrature;leastsquaresfitting;AdamsandRunge-Kuttaformulas;Richardsonextrapolation1940-1970floatingpointarithmetic;Fortran;finitedifferences;finiteelements;FFT;simplexalgorithm;MonteCarlo;orthogonallinearalgebra;splinefunction1970-2000quasi-Newtoniterations;adaptivity;stiffODEsolvers;softwarelibraries;Matlab;multigrid;sparseanditerativelinearalgebra;spectralmethods;interiorpointmethods計算數學未來50年的展望將更多的通過聲音，而不是鍵盤向計算機傳遞信息，而計算機將更多地以圖象而不是數字反映結果數值計算將更具有適應性、迭代性、靈活性。計算能力大得驚人數值計算中更具智能性算法：一系列近似計算步驟的組成，目的是找到問題的近似解算法的特征：收斂性、穩(wěn)定性計算量：一個算法所需的乘除運算總次數，單位是flop.計算量是衡量一個算法好壞的重要標準算法和收斂矩陣乘積AB的計算量分析a11a12a13…a1na21a22

a23…a2n...

...

…...am1am2

amm-1…amnb11b12b13…b1sb21b22

b23…b2s...

...

…...bn1bn2

bnn-1…

bns=[cij]ms因為cij=aik

bkj

計算量為n所以上面AmnBns的計算量為N=mn

s誤差分析ErrorAnalysis誤差的來源誤差誤差限有效數字相對誤差和絕對誤差誤差的傳播在近似計算中需要注意的問題目次模型誤差觀測誤差舍入誤差截斷誤差1.1

誤差的來源

計算機字長有限，一般實數不能精確存儲，于是產生舍入誤差。例如：在10位十進制數限制下：1÷3＝0.3333333333本應1÷3＝0.3333333333……1.0000022-1.000004=0本應1.0000022-1.000004

=1.000004000004-1.000004

=0.0000000000041.1.1

舍入誤差(Round-offErrors)舍入誤差很小，本課程將研究它在運算過程中是否能有效控制。1.1.2截斷誤差(TruncationError)用近似的值去代替數學上的準確值帶來的誤差。例如：

泰勒級數?零階近似:?一階近似:?二階近似:完全的泰勒級數:余項(n階近似)::介于

xiand

xi+1

x

=

xi+1-

xi

余項：Taylor級數表示為:截去的部分?零階近似:

截斷誤差

:?

一階近似

Rn:零階近似Rn:斜率:1.2

誤差誤差限有效數字[Def1.1]若用x*表示x準確值的一個近似值。則此近似值x*和準確值x的差稱為誤差，用e*來表示

e*＝x*－x[Def1.2]若

|e*|＝|x*－x|≤ε*ε*稱為近似值x*的誤差限。[例1.2]已知x*=π=3.14159…，求近似值x1=3.14，x2=3.142，x3=3.1416的誤差限。[解]所以誤差限

ε1=0.002，ε2=0.0005，ε3=0.000008有效數字[Def1.3]若用x的近似值x*的誤差限是某一位上的半個單位，該位到x*的第一位非零數字共有n位，則稱x*有n位有效數字若用x*表示x的近似值，并將x*表示成

x*＝±0.a1a2…an×10m若

|x*－x|≤0.5×10m－n則近似值x*有n位有效數字(1.1)[例1.3]

設x*=0.0270是某數x經“四舍五入”所得，則誤差|e(x*)|不超過x*末位的半個單位，即：

|x*－x|≤0.5×10-4

又x*=0.27×10-1,故該不等式又可寫為

|x*－x|≤0.5×10-1-3由有效數字定義可知,x*有3位有效數字,分別是2,7，0。[例1.4]

設x＝32.93,x*＝32.89,則

|x*－x|＝0.04＜0.05＝0.5×10-1即

|x*－x|≤0.5×102-3由有效數字定義可知,x*有3位有效數字,分別是3,2,8。由于x*中的數字9不是有效數字，故x*不是有效數。1.3相對誤差和絕對誤差設x——準確值x*——近似值稱為近似值x*的相對誤差實用中，常用表示近似值x*的相對誤差，稱為相對誤差限相應的，e*稱為絕對誤差，ε稱為絕對誤差限有效數位與誤差的關系有效數位n越多，則絕對誤差|e*|越小形如(1.1)式的近似數x*具有n位有效數字，則其相對誤差限可取為基本算術運算設x*和y*分別是x和y的近似值，把它們的誤差近似地看做是相應地微分，即

dx≈x*－x，dy≈y*－y則

d（x±y）＝dx±

dy d（xy）＝xdy±

ydx d（x/y）＝（－xdy＋ydx

）/y21.4

誤差傳播(1.3)和(1.4)給出了由自變量的誤差引起的函數值的誤差的近似式（誤差傳播）。一元函數設y＝f(x)，若x的近似值是x*，用f(x*)去近似f(x)的誤差可用Taylor公式估計

(1.3)(1.4)多元函數情形由多元函數的Taylor展開公式類似可得

(1.5)

(1.6)

(1.8)

(1.9)[例1.5]測得某桌面的長a的近似值a*=120cm,寬b的近似值b*=60cm。若已知|e(a*)|≤0.2cm,|e(b*)|≤0.1cm。試求近似面積s*=a*b*

的絕對誤差限與相對誤差限。解:面積s=ab,在公式（1.5）中,將y=f(x1,x2)

換為s=ab,則相對誤差限為1.5

在近似計算中需要注意的問題1.盡量簡化計算步驟，減少乘除運算的次數

例如，計算多項式通常運算的乘法次數為若采用遞推算法,

則乘法次數僅為n.又如2.防止大數“吃掉”小數當|a|>>|b|時，盡量避免a+b

。例如，假設計算機只能存放10位尾數的十進制數，則

108+0.04=1083.盡量避免相近數相減例如，當x很大時，應當x接近于0時，應4.避免絕對值很小的數做分母當|b|<<|a|時，應盡量避免a/b

5.

選用數值穩(wěn)定性好的算法，以控制舍入誤差高速增長第2章線性代數方程組解法本章研究的對象是n

階線性代數方程組2.1對象(2.1)線性系統(tǒng)廣泛存在于工程、科學以及社會科學、商業(yè)和經濟問題的定量分析等領域中用矩陣和向量的記法來表示，(2.1)式可寫成對象（續(xù)）(2.2)其中A＝(aij)是方程組(2.1)的系數aij構成的n×n階矩陣，稱為系數矩陣。B＝{bi},X＝{xi}是n維向量，X是未知量，B稱為右端項。使方程組(2.1)中每一個方程都成立的一組數x1*,x2*,…,xn*稱為式，(2.1)的解，把它記為向量的形式，稱為解向量。

我們總是希望方程組有解,且有唯一解.由線性代數的克萊姆(cramer)規(guī)則可知,如果方程組(2.1)的系數矩陣A的行列式不等于零,那么,這個方程組有唯一解,而且它們可以表示為

xi=Di/D(i=1,…,n)這里,Di是指D中第i列元素用右端b1,…bn代替構成的行列式.

如果方程組(2.1)有唯一解,我們按上面的等式求解,就必須計算n+1個n階行列式.由行列式的定義,n階行列式包含有n!項,每一項含有n個因子,計算一個n階行列式就需要做(n-1)n!次乘法.而我們一共要計算n+1個n階行列式,共需做(n2-1)n!次乘法.此外,還要做n次除法才能算出xi(i=1,…n).也就是說,用這個辦法求解就要做

N=(n2-1)n!+n次乘除法運算,這個計算量是大得驚人的.例如,當n=10(即求解一個含10個未知量的方程組),乘除法的運算次數共為32659210次;

消(元)去法是求解線性方程組(2.2)和滿秩矩陣的逆陣A-1的一種直接方法.盡管它比較古老,但它具有演算步驟,推算公式都系統(tǒng)化的特點(對其中主元素消去法,還可以證明是穩(wěn)定的).因此,它至今仍然是求解方程組的一種有效的方法.

消去法可以引出幾種計算方法,下面按三角形方程組和一般線性方程組的順序來討論。消去法上三角方程組的一般形式是：對于(2.1)式，有以下運算可簡化方程組線性方程組的解法有兩類：直接法：即在沒有舍入誤差的情況下，用有限步的四則運算得出精確解的方法。但實際運算中舍入誤差不可避免，此類方法也只能得到近似解。目前常用的是列主元消去法和矩陣三角分解法迭代法：先給一個初始值，按一定法則逐步求解出各個更準確的近似值的方法。目前常用的有Jacobi法、Seidel迭代法、松弛法和梯度法將A和B寫在一起，稱為增廣矩陣將例2.1用增廣矩陣的變化寫成高斯消去法的求解過程分為兩個階段:首先，把原方程組化為上三角形方程組，稱之為“消元”過程；然后，用逆次序逐一求出三角方程組(原方程組的等價方程組)的解，并稱之為“回代”過程。2.2高斯消去法(Gaussianelimination)消元：將(2.1)式寫成矩陣形式(2.3)(2.4)第1步:若a11(1)

≠0，用第二個方程減去第一個方程乘以a21(1)/

a11(1),用第三個方程減去第一個方程乘以a31(1)/

a11(1)…則有矩陣形式(2.5)(2.6)第2步:若a22(2)

≠0，用第三個方程減去第二個方程乘以a32(2)/

a22(2),用第四個方程減去第二個方程乘以a42(2)/

a22(2)…則有矩陣形式(2.7)(2.8)第k步:若akk

(k)

≠0，用第k+1個方程減去第k個方程乘以ak+1k(k)/

akk(k)…則有矩陣形式(2.9)(2.10)重復n－1次，得到等價的上三角形方程組矩陣形式(2.11)(2.12)以上過程把系數矩陣A(1)變成上三角矩陣A(n)，稱之為消元，計算公式可歸納為(2.13)回代(2.12)2.3主元素消去法因此，x1＝0x2＝1因此，x1＝1x2＝1精確解，x1＝10000/9999x2＝9998/9999在做除法運算時，選取絕對值大的作分母。

——主元素消去法的基本思路。列主元素消去法列主元素消去法基本思想1.用高斯消去法求解線性方程組時,應避免小的主元.在實際計算中,進行第k步消去前,應該在第k列元素aik

(i=k,…,n)中找出絕對值最大者,例如

∣a∣=max∣a∣2.再把第p個方程與第k個方程組進行交換,使apk(k-1)成為主元.我們稱這個過程為選主元素.由于只在第k列元素中選主元素,通常也稱為按列選主元素(或稱部分選主元).3.如果在第k步消去前,在第k個方程到第n個方程所有的xk到xn的系數a(i=k,…,n;j=k,…,n)中,找出絕對值最大者,例如 ∣a∣=max∣a∣再交換第k,p兩個方程和第k,q兩個未知量的次序,使a成為主元素.稱這個過程為完全選主元素。4.不論是哪種方式選出主元素,而后再按上面介紹的計算步驟進行消去的計算,一般都稱為主元素高斯消去法。在實際計算中,常用按列選主元素的高斯消去法。(k-1)(k-1)pk(k-1)ikk≦i≦n(k-1)pq(k-1)ijk≦i,j≦n(k-1)ij(k-1)pq高斯消去法的乘除總運算分析如下：消元次數k消元乘法次數消元除法次數回代乘除法總次數

1n(n-1)n-12(n-1)(n-2)n-2...k(n-k+1)(n-k)n-k..n-12*11n(n+1)/2

故高斯消去法的計算量為N=n(n2-1)/3+n(n-1)/2+n(n+1)/2=n3/3+n2-n/3

當N充分大時為n3/32.4高斯消去法的計算量分析

方法的特點

由具體計算結果可知,全主元素法的精度優(yōu)于主元素法,這是由于全主元素是在全體系數中選主元,故它對控制舍入誤差十分有效.但全主元素法在計算過程中,需同時作行與列的互換,因而程序比較復雜,計算時間較長.列主元素法的精度雖稍低于全主元素法,但其計算簡單工作量大為減少,且計算經驗與理論分析均表明,它與全主元素法同樣具有良好的數值穩(wěn)定性,列主元素法是求解中小型稠密性方程組的最好方法之一。

選主元消去法(包括解線性方程組的所有直接的方法)比較適用于中小型方程組.對高階方程組,即使系數矩陣是稀疏的,但在計算中很難保持稀疏性,因而有存儲量大，程序復雜等不足,所幸的是這一缺點可用迭代法解決.

另外,高斯選主元消去法還可技巧性的解決一些特殊線性方程組。由于誤差的影響,計算過程中可能出現一些壞現象,要盡可能防止,表現在求解線性方程組的消元法比較上,則應該注意要簡化運算,減小運算次數,提高效率;提高數值穩(wěn)定性等.

在計算過程中，由于計算機字長的有限性，不可避免地產生舍入誤差。同時，由于所求問題的初始數據（例如線形方程組的系數矩陣和右端項系數）往往是帶有一定誤差的。因此計算結果總是不可避免地帶有誤差，或者說，如果初始數據有擾動，勢必將帶來具有一定誤差的計算結果。就拿Ax=b來說，由于觀測或計算等原因，線性方程組兩端的系數A和b都帶有誤差A和b，這樣實際建立的方程組是近似方程組(A+A)(x+x)=b+b。對近似方程組求出的解是原問題的真解x加上誤差x,即x+x。而x是由A及b引起的，它的大小將直接影響所求解的可靠性。這種解依賴于方程組系數的誤差A及b的問題，稱為線性方程組解對系數的敏感性。2.5線性方程組解對系數的敏感性方程組

此方程組的準確解為x1=0,x2=-1?，F將其右端加以微小的擾動使之變?yōu)椋?/p>

經計算可得準確解為x1=2,x2=-3.

這兩個方程組的解相差很大，說明方程組的解對常數項b的擾動很敏感。

病態(tài)方程組：如果方程組AX=b由于A或b的