TOC\o"1-5"\h\z第1章誤差分析與數(shù)值計算 3引言 3絕對誤差與相對誤差、有效數(shù)字…9近似數(shù)的簡單算術(shù)運算 12數(shù)值計算中誤差分析的一些原則，4第2章非線性方程（組）的近似解法........17引言 17根的隔離 18對分法 19迭代法 21弦截法 24弦截法 25§1.7用牛頓法解方程組 26本章小結(jié) 28第3章線性方程組的解法 29引言 29高斯消去法 31矩陣的LU分解 35對稱矩陣的LDlJ分解 36線性方程組解的可靠性 37簡單迭代法 40本章小結(jié) 46第4章矩陣特征值與特征向量的計算…47引言 47基法和反暴法 48雅可比方法 49QR方法* 51本章小結(jié) 52第5章插值與擬合 53引言 53插值多項式的存在和唯一性 54拉格朗日插值多項式 55均差插值公式 57差分等距結(jié)點插值公式 59愛爾米特插值公式 61分段低次插值 62三次樣條函數(shù) 63曲線擬合的最小二乘法 6871第6章數(shù)值積分和數(shù)值微分 72引言 72牛頓一科特斯型積分公式 73復(fù)合求積公式 75龍貝格求積公式 79高斯求積公式 80二重積分的數(shù)值積分法 84數(shù)值微分 8587第7章常微分方程的數(shù)值解法 88引言 88歐拉法和改進的歐拉法 89龍格-庫塔方法 90線性多步法 96算法的穩(wěn)定性與收斂性 98微分方程組和高階微分方程解法“99本章小結(jié) 101第1章誤差分析與數(shù)值計算§1.1引言1、課程任務(wù)和目的：在第七屆國際軟件工程學(xué)術(shù)會議上，“計算方法”被列入應(yīng)用方法學(xué)的研究領(lǐng)域，強調(diào)了計算方法的研究應(yīng)用與軟件方法學(xué)的研究密切結(jié)合。這就說明了計算方法與軟件之間的聯(lián)系以及在應(yīng)用軟件研制中的地位與作用，計算方法是研究各種數(shù)學(xué)問題求解的數(shù)值計算方法。在計算機成為數(shù)值計算的主要工具的今天，則要求研究適合于計算機使用的數(shù)值計算方法。計算方法就是研究用計算機解決數(shù)學(xué)問題的數(shù)值方法及其理論，它的內(nèi)容包括函數(shù)的數(shù)值逼近、數(shù)值微分與數(shù)值積分、非線性方程值解、線性方程組數(shù)值解、常微和偏微數(shù)值解等，即都是以數(shù)學(xué)問題為研究對象的。因此，計算方法是數(shù)學(xué)的一個分支，只是它不象純數(shù)學(xué)那樣只研究數(shù)學(xué)本身的理論，是把理論與計算緊密結(jié)合，著重研究數(shù)學(xué)問題的數(shù)值方法及其理論，計算方法是計算機應(yīng)用和軟件研制開發(fā)的重要組成部分，通過本課程的學(xué)習(xí)和上機實習(xí)，使學(xué)生掌握利用計算機進行科學(xué)計算的基本理論和基本方法，并且學(xué)會將基本理論和基本方法應(yīng)用于軟件開發(fā)以及軟件研制。2、本課程基本要求(1)掌握方法的基本原理和思想。(2)掌握方法處理的技巧及與計算機的結(jié)合。(3)掌握誤差分析，收斂性及穩(wěn)定性的基本理論。(4)學(xué)會進行可靠的理論分析，對近似計算要確保精度要求，要進行誤差分析。(5)通過例子，學(xué)習(xí)使用各種計算方法解決實際計算問題。(6)通過上機實踐，能編寫算法和實現(xiàn)算法。(7)掌握數(shù)值計算中一些最基本、最常用的計算方法和算法。3、本課程與各課程的關(guān)系：由于本課內(nèi)容包括了微積分、代數(shù)、常微分方程的數(shù)值方法，學(xué)生必須掌握這幾門課的基本內(nèi)容才能學(xué)好這一課程，同時，學(xué)習(xí)此課程還必須具備計算機系統(tǒng)的初步知識，掌握一門常用的高級語言，如：BASIC,PASCAL,C語言等，并須具備一定的編程能力。4、本課程的特點：(1)面向計算機，要根據(jù)計算機特點提供實際可行的有效算法。即算法只能包括加、減、乘、除運算和邏輯運算，是計算機能直接處理的。(2)有可靠的理論分析，能任意逼近并達到精度要求，對近似算法要保證收斂性和數(shù)值穩(wěn)定性，還要對誤差進行分析，而且都是建立在相應(yīng)數(shù)學(xué)理論基礎(chǔ)上的。(3)有好的計算復(fù)雜性。時間復(fù)雜性好是指節(jié)省時間：空間復(fù)雜性好是指節(jié)省存儲量。這也是建立算法時要研究的問題，因為它關(guān)系到算法能否在計算機上完成。(4)要有數(shù)值實驗。即任何一種算法除了從理論上要滿足上述三點外，還要通過數(shù)值實驗證明是行之有效的。計算方法最基本的立足點是容許誤差，在誤差容許的范圍內(nèi)對某一數(shù)學(xué)問題進行近似計算，得到能滿足耍求的近似結(jié)果?，F(xiàn)實世界中誤差是普遍存在的，由于世界上沒有絕對精確的量具（絕對精確的量具是沒有刻度的），因此人類通過量具采集的數(shù)據(jù)都是近似值，另一方面，我們的生產(chǎn)、實驗工具都不是絕對精確的，這就使得人類在生產(chǎn)和科學(xué)實驗中必需容許誤差。計算機的應(yīng)用可以分為二個方面，即數(shù)值計算和非數(shù)值計算。利用計算機進行數(shù)值計算的過程如下圖所示:實際問題―?數(shù)學(xué)模型―?計算方法―?程序設(shè)計——?上機求解在上圖中，計算方法的任務(wù)是：由建立的數(shù)學(xué)模型給出可編程并由計算機能完成的計算方法，然后編程和上機求解。由于計算方法是編程后可由計算機求解的近似計算方法，如何確保近似解的精度顯得尤為重要，必須深入討論有關(guān)誤差的基本概念和基本理論，為近似計算的精度分析打下基礎(chǔ)。1、誤差的來源（種類）誤差的來源主要有以下四種（1）模型誤差：建立數(shù)學(xué)模型時的誤差。例如：在求重量的數(shù)學(xué)模型G=m*g中，重量G不是僅與質(zhì)量和重力加速度有關(guān)，它還與溫度、測量地點的海拔、地層結(jié)構(gòu)等眾多因素有關(guān)，為了使模型較為簡單和實用，采用抓住主要矛盾的方法，去掉了大量對重量影響不大的次要因素，建立了上述重量的近似模型，由此產(chǎn)生了模型誤差。（2）觀測誤差：采集數(shù)據(jù)時的誤差。采集數(shù)據(jù)時，通常是依靠儀器和量具，由于沒有絕對精確的儀器和量具，因此采集的數(shù)據(jù)有誤差，此誤差稱為觀測誤差。（3）舍入誤差：由于計算機字長有限而產(chǎn)生的誤差。硬件再發(fā)展，計算機的字長總是有限的，在計算過程中，當數(shù)據(jù)的長度超過了計算機的字長時，計算機就會進行四舍五入，由此產(chǎn)生的誤差稱為舍入誤差。（4）截斷誤差：無限形式的有限化而產(chǎn)生的誤差。在計算中有時會運用無限形式的計算公式，例如臺勞公式：f(x)=f(x0f(x)=f(x0)+f-^X°\x-x0)+---+f(n)(x0)n!(x-x0)n+???顯然此公式無法進行計算，因此必需根據(jù)實際需要，從某?項起將后面的各項截斷，即f(x)=f(xo)+彗以(x—xo)+…+吧，。％—xo)n

1! n!由此產(chǎn)生的誤差稱為截斷誤差?！?.2絕對誤差與相對誤差、有效數(shù)字為描述方便，首先約定x*是精確值x的近似值。引入誤差的概念，其目的是為了衡量近似值x*的好壞。(1)絕對誤差：x*-x由于精確值X通常無法確定，因此絕對誤差無法計算，由此引入絕對誤差限的概念。絕對誤差限：絕對誤差的一個上界。即：若|x*-x|<e,則稱e為X*的絕對誤差限。絕對誤差限的性質(zhì)是：A.不唯一這是因為|x*-x|的上界是不唯一的。B.可確定只要我們對x*的實際背景有一定的了解，就不難確定|x*-xI的上界。例如，X*表示身高，則Ix*-xI的上界可為3米。當X*是你求出的，那么為了說明你的工作認真，你一定會將|x*-x|的上界估計得盡量小，因此在這種意義上絕對誤差限可用來衡量x*的好壞。由于絕對誤差限沒有考慮問題的規(guī)模，因此有時它也不能衡量x*的好壞。例如：x是地球與太陽的距離，y是分子中二個原子間的距離，若|x*-x|41公里，|y*-y|4l厘米，則并不能說y*比x*精確。由此引入相對誤差和相時誤差限的概念。(2)相對誤差：(x*-x)/x*相對誤差限：相對誤差絕對值的一個上界。3、有效數(shù)字這里我們必須搞清楚什么是有效數(shù)字以及如何確定X*有幾位有效數(shù)字。(1)有效數(shù)字的定義若|x*-x|<x*的某一位的半個單位，則稱X*精確到這,位，并從這一位開始，一直到前面第一個不為零的數(shù)都是X*的有效數(shù)字。此定義實際上定義了什么叫精確到某一位和什么叫有效數(shù)字。例如:若x*精確到小數(shù)點后第3位，即指|x*-x|<0.5x10-3o(2)有效數(shù)字的判定方法方法一：四舍五入此方法首先確定x*是由x的哪一位四舍五入產(chǎn)生的，然后從這一位的前一位開始一直到前面第一個不為零的數(shù)都是x*的有效數(shù)字。例1若x=0.872596,x*=0.87,求x*的有效位數(shù)。解：???x*是由x的小數(shù)點后第三位四舍五入產(chǎn)生的，所以x*有二位有效數(shù)字。注意，方法一判定有效數(shù)字很簡單，但有時會失效。例如，若x=0.272987x*=0.273102,此時無法用方法-確定x*的有效位數(shù)，原因是x*不是由x四舍五入產(chǎn)生的，在這種情況下，必須用有效數(shù)字的定義來確定X*的有效位數(shù)。即方法二：用定義此方法首先計算Ix--x|,再判斷它小于等于X*的哪一位的半個單位，然后從近一位開始，一直到第一個不為零的數(shù)都是有效數(shù)字。例2若x=0.62073,x*=0.6207,確定x*的有效位數(shù)。解：因為|x*-x|40.000340.5x10-4,x*精確到小數(shù)點后第4位，所以x*有四位有效數(shù)字。例3若x=0.080199,x*=0.802,確定x*的有效位數(shù)。解：因為|x*-x1=0.0000140.5x10-3,所以推出x*有三位有效數(shù)字。例4若x=6.28936,x*=7.3132,確定x*的有效位數(shù)。解：|x*-x|=0.02357<0.5xl0-1,所以x*有二位有效數(shù)字。近似數(shù)的簡單算術(shù)運算近似數(shù)的加法設(shè)有n的近似數(shù)Xk*>0(k=l,2,...n)淇準確值為xk>0,(k=l,2,...n)x*=fx；的絕對誤差限E(x*)k=\E")=f々-f工這(々-x；)=f )k=\k=lk=\ k=l|￡(x)|<^|￡(x；)|A=1(1)和的絕對誤差等于各項絕對誤差之和。(2)和的絕對誤差限不超過各項絕對誤差限之和。類似的可以得到：和的相對誤差限介于各項相對誤差限的最小者與最大者之間。近似數(shù)的乘法結(jié)論：乘積的相時誤差限不超過各項相對誤差限之和。似數(shù)的除法結(jié)論：商的相對誤差限不超過被除數(shù)與除數(shù)相對誤差限之和。134近似數(shù)的鼎和根(見教材p9)近似數(shù)的對數(shù)(見教材p9)近似數(shù)的減法結(jié)論：差的絕對誤差等于各項絕對誤差之和。注意兩個幾乎相等的近似數(shù)相減會使結(jié)果的有效數(shù)字損失，影響整個計算工作的準確性，應(yīng)盡量避免。7101-7^00=^=2^^=l-cos(x)=2sin2(-)ln(x)-ln(x,)=ln(^-)V3.01+V3.00 2 x2(當|x|很小時，X］與X2很接近時)數(shù)值計算中誤差分析的一些原則為保證計算結(jié)果的高精度，在進行數(shù)值計算時應(yīng)遵循下述幾個原則。(1)在進行除法時，要避免除數(shù)的絕對值＜＜被除數(shù)的絕時值。①為什么要“避免”？若不“避免”，則除出的結(jié)果很大，由于計算機字長有限，它裝不下，因此會進行四舍五入，一個很大的數(shù)進行四舍五入時舍去的部分也會很大，這會使舍入誤差變大。②怎樣“避免”？因為用戶只關(guān)心最后的計算結(jié)果，當中間計算過程中出現(xiàn)了除數(shù)的絕對值＜＜被除數(shù)的絕對值時，就應(yīng)該換一種計算方法，以避免這種情況的發(fā)生，以后我們將會針對具體的計算問題來討論“避免”的方法。(2)在進行減法時，要避免二個相近的數(shù)相減。①為什么要“避免”？若不“避免”，就可能失去大量的有效數(shù)字，例如：若a=30001和b=30000都有五位有效數(shù)字，因為a-b=l,所以結(jié)果至多有1位有效數(shù)字。②怎么“避免”？“避免”的思路與第1個原則中“避免”的思路相同，須針對具體計算問題來討論。（3）要防止“大數(shù)吃小數(shù)”①什么是“大數(shù)吃小數(shù)”？我們用一個例子為說明。n計算8756294874+Zaj，其中『IO?。,（KasKT6。i=l此題是一個很大的數(shù)與很多很小的數(shù)相加，若采用將大數(shù)依次與a1,a2,…,an相加，由于計算機字長有限，因此在與內(nèi)相加時會進行四舍五入將為舍去，這樣，最后的結(jié)果仍是大數(shù)，這就是大數(shù)將ai,a2,…,a”吃掉了。②為什么要“避免”？盡管每個小數(shù)都很小，但它們很多，可能它們的和比大數(shù)還大，而最后計算工結(jié)果為大數(shù)，顯然誤差可能很大。③怎樣“避免”？有的同學(xué)提出先將小數(shù)相加，然后再與大數(shù)相加，這個思路是對的，但有一個漏洞，因為小數(shù)相加到一定程度也會變成大數(shù)，它也開始吃小數(shù)了?？梢圆扇》植肯嗉拥姆椒ń鉀Q。第2章非線性方程（組）的近似解法引言方程Rx）=O的解稱為方程的根。也叫做函數(shù)Rx）的零點。方程求根大致包括三個問題（1）方程有沒有根？如果有根，有幾個根？（2）哪里有根？求有根的區(qū)間，區(qū)間內(nèi)的任意一點作為根的近似值。（3）根的精確化，已知一個根的近似值后設(shè)法逐步把根精確化，直到足夠精確為止。本課程主要研究問題（2）和（3）.根的隔離求方程f(x)=()的解的近似值時，首先要確定若干個區(qū)間，使每個區(qū)間內(nèi)只有的一個根，這個步驟稱為根的隔離。對一般的方程，根的隔離有兩種方法(1)試值法。求出f(x)在若干點上的函數(shù)值，觀察函數(shù)值符號變化的情況，從而確定隔根區(qū)間。(2)作圖法。畫出尸f(x)的草圖，觀察曲線產(chǎn)f(x)與x軸交點的大致位置，從而確定隔根區(qū)間。例1.2.1討論方程f(x)=2x*-4x-+4x+2=0的根的位置。fl=inline('2*xA3-4*xA2+4*x+2'),fplot(fl,[-1,1]), f2=inline('log(x)-1/x'),fplot(f2,[l,2])例1.2.2將方程xlog(x)=1的根進行隔離。設(shè)有方程1x)=0在(ab)內(nèi)有且僅有?個根a,這時有f(a)f(b)<0可用對分法求a的近似值，方法如下(1)準備：計算區(qū)間(ab)兩個端點的函數(shù)值f(a),Rb)(2)對分：取c=(a+b)/2為(ab)的中點，計算的)(3)判斷:如果f(c)=0,則c為f(x)=0的根，否則檢驗:若fi［c)f(a)<0,則方程的根位于［ac］內(nèi),用c代替b,若f(c)f(b)<0,則方程的根位于［cb］內(nèi),用c代替a。(4)檢驗：若|b-aRe(e為精度要求)此時計算結(jié)束，取a=c,否則轉(zhuǎn)(2)。例2.3.1用對分法求方程Rx)=x3+2x-5=0在［12］內(nèi)的根，［e=10力。有根區(qū)間f=inline(,xA3+2*x-51)1.00002.0000ep=le-3;a=l;b=2;1.00001.5000whileabs(b-a)>ep1.25001.5000c=(a+b)/2;1.25001.3750iff(c)*f(a)>01.31251.3750a=c;1.31251.3438else1.32811.3438b=c;1.32811.3359end1.32811.3320[a,b,abs(b-a)],pause1.32811.3301end1.32811.3291fplot(f,[l2])方程的解x=1.3286§2.4迭代法設(shè)有方程f(x)=O在［ab］上有且僅有一個根a,可用迭代法求a的近似值，方法如下(1)將方程電)=0寫成迭代形式x=<p(x)(2)在［ab］上任取一個初始值xo.(3)計算xi=<p(xo)(4)若|xLX()|<e(e為精度要求)，此時計算結(jié)束a=x”否則令x()=X|轉(zhuǎn)(3)。定理241(收斂定理)設(shè)函數(shù)<p(x)在［ab］上連續(xù)，在(ab)內(nèi)可導(dǎo)，且滿足(1)當x￡［ab］時,(p(x)e［ab］?(2)當xG［ab］時，|(p'(x)|<m<1,m為一個常數(shù)。則以下結(jié)論成立：(1)在［ab］上<p(x)有且僅有一個根a。(2)對任意x()e［ab］,由迭代公式Xn+i=(p(Xn),n=0,l,2,…產(chǎn)生的數(shù)列對仁［ab］,且有xn-?a(n->°°)otn' m(3)成立誤差估計式區(qū) 1王一與11。一工〃區(qū) |xn-xn_}\\—m \—m討論：(1)此定理說明只要|Xn-Xn"足夠小，就可以保證X。充分接近方程的根電所以在實際計算時，用條件|Xn-Xn-^e控制迭代過程的結(jié)束。(2)定理2.4.1指當對任意x°w[ab]作為初始條件，迭代過程都收斂，這種形式的收斂定理稱為大范圍收斂定理，是一個充分條件。當此條件不成立時，往往可以取初始值X。接近于方程的根a,使迭代過程收斂。(見P21定理2.4.2)。(3)從定理2.4.1中可以看出，收斂速度與m值有很大關(guān)系，m越小，收斂速度越快。當m接近于1時，收斂速度很慢。根據(jù)卬(x)￡m<l,所以采用迭代法求解方程的根時，應(yīng)適當?shù)臉?gòu)造(p(x)使|(p'(x)|在區(qū)間[ab]內(nèi)盡量的小。(見下面的例子)例2.4.1用迭代法解方程*=10x-2,x0=l分別采用迭代格式x=10x-2和x=k>g(x+2),觀察兩個計算過程的區(qū)別。e=le-3迭代過程：1.00000.4771 0.3939 0.3791 0.3764 0.3759迭代6次x=0.3759Mnlme('logl0(x+2)')X=1x=f(x)例2.4.2用迭代法求方程f(x)=x3+2x-5=0的根，x0=l[xn+|=^5-2xn迭代過程:1.00001.44221.28371.34491.32201.33061.32741.32861.3281迭代9次乂=1.3281Mnline(,(5-2*x)A(l/3y)x=lx=Rx)§2.5牛頓迭代法牛頓法是解方程f(x)=O的重要方法，它也是一種迭代法。設(shè)有方程f(x)=O在［ab］上有且僅有一個根a,可用牛頓法Xn+產(chǎn)X「f(Xn)/f'(Xn)求a的近似值，具體計算步驟如下(1)求出f(x)的導(dǎo)數(shù)f'(X),在［ab］上任取一個初始值Xo。(2)計算xi=x<)-fi［xo)/f'(xo)(3)若|X|-x()|ve(e為精度要求)，此時計算結(jié)束a=X］,否則令x()=X|轉(zhuǎn)(2)。例2.5.1用牛頓法解方程f(x)=x3-2x2>4x-7=O在［34］內(nèi)的根［x0=4］,迭代過程： 4.0000 3.6786 3.6329 3.6320迭代4次x=3.6320f=inline('xA3-2*xA2-4*x-7')fd=inline('3*xA2-4*x-4')x=4x=x-f(x)/fd(x)§2.6弦截法弦截法也是一種是解方程f(x)=O的迭代法，它的特點是不需要計算f(x)的導(dǎo)數(shù)f'(x),且收斂速度也相當快,是工程計算中常用的算法之一。設(shè)有方程4x)=0在［ab］上有且僅有一個根a,可用弦截法求a的近似值，方法如下(1)求函數(shù)fi[x)在區(qū)間[ab]的兩個端點的函數(shù)值Nxo),f(X|)，其中a=x(),b=X](2)計算x2=Xi-f(x()[X!-x0]/[Rxi)-f(x?)](3)若|X2-X]|<e(e為精度要求)，此時計算結(jié)束a=X2，否則令x()=X|，Xi=x?轉(zhuǎn)(2)。例2.6.1用牛頓法解方程f(x)=x3-2x2-4x-7=0在[34]內(nèi)的根。f=inline(*xA3-2*xA2-4*x-7,)x0=3,xl=4x2=xl-f(xl)*(x1-x0)/(fi(xl)-f(x0))x0=xl,xl=x2§1.7用牛頓法解方程組設(shè)有非線性方程組u(x,y尸0,v(x,y)=0,在(xn,yn)按臺勞級數(shù)展開，取展開式的前兩項得到卜(—+也產(chǎn)**°)+吟山0?)=0卜(—+也產(chǎn)**°)+吟山0?)=0

ok dyv(xn,yn)+aV(^,yn)(x-xn)+aV(^,yn)(y-yn)=0

ox dy加(Xn?n)

dx

加加―)

dxHu(Xn，yn)HyAv(Xn，yn)

Hy其中(Xn,yn)是根的第n次近似值，如果Jn#0方程組的第n+1次近似值(Xn+I,%+|)可用以下公式計算1n+1=Xn+~JjlHu(Xn，yn)

Hy1n+1=Xn+~JjlHu(Xn，yn)

Hy

報小0)

打u(xn,yn)v(xn,y?)1yn+i=yn+—Jnu(xn,yn)v(xn,yn)^u(xn,yn)dx次心’)

dxx0=l;y0=1.4;x0=l;y0=1.4;JO=det([3*xOA2,3*yOA2;4*xOA3,2*yO]);fbrk=l:10uO=xOA3+yOA3-4;vO=xO八4+yO人2?3;Jx=det([3*yOA2,uO;2*yO,vO]);Jy=det([uO,3*xOA2;vO,4*xOA3]);xl=xO+Jx/JO;yl=yO+Jy/JO;xO=xl;yO=yl;[xOyO],pauseend例1.7.1用牛頓法解方程組u(x,y)=x3+y3-4=0v(x,y)=x4+y2-3=0(3u/3x=3x2,du/dy=3y2,dv/dx=4x3,9v/dy=2y)迭代初值x0=l,yo=l.4例172用牛頓法解方程組u(x,y)=2x3-y2-l=0v(x,y)=xy3-y-4=0(du/3x=6x2,3u/3y=-2x,dv/3x=6y3,3v/dy=2xy2-1)迭代初值Xo=l.2,yo=l.7o例1.7.3用牛頓法解方程組u(x,y)=x-cos(y)=0,v(x,y)=y-sin(x)=0(du/dx=1,du/3y=sin(y),dv/dx=-cos(x),3v/3y=l)迭代初值xo=O,yo=0o本章小結(jié)為了比較各種迭代方法的收斂速度，我們引入收斂階的概念。設(shè)迭代過程X,收斂于方程X=(p(X)的根a,令en=Xn-a,e1,稱為迭代誤差，如果存在實數(shù)P21和非零常數(shù)K,使得lim&J=K,則稱該迭代過程為P階收斂的。P=1稱為線性收斂，P>1稱為超線性收斂，P=2稱為平方收斂，顯然P越大，迭代過程收斂的越快?？梢宰C明當a是方程f(x)=O的單根時，牛頓法是平方收斂的。當a是方程4x)=0的重根時，牛頓法僅為線性收斂。對分法：收斂速度與公比為1/2的等比級數(shù)相同。收斂速度最慢，但簡單有效，不存在發(fā)散問題。弦截法：弦截法的收斂階P=1.618。收斂速度次之，不需要計算f(x)的導(dǎo)函數(shù)，有發(fā)散問題。牛頓法：牛頓法是平方收斂，收斂速度最快，但要計算f(x)的導(dǎo)函數(shù)，有發(fā)散問題。第3章線性方程組的解法§3.1引言在科學(xué)實驗和工程設(shè)計中，經(jīng)常用到解線性方程組的問題。本章討論用計算機求解線性方程組的兩類主要方法：I'［接法和迭代法。解線性方程組的一般表達式aMi+ag+…+amXn=5a21X|+a22x2+---+a2nxn=b2根據(jù)矩陣的運算的性質(zhì)有3b,an|Xi+an2X2+…+annXn=b簡記為Ax=b其中A=4b,方程組AxH有唯一解的充分必要條件是det(A/O。理論上，可以用克萊姆法則x產(chǎn)AJA,(k=l,2..n)其中A=det(A),是A中第k列換成向量所得到的矩陣的行列式。計算需耍的乘法次數(shù)為(n-l)(n+l)!,當n較大時,計算量很大。這種方法效率低，在實際應(yīng)用中很少使用。解線性方程組的方法可以分為兩類：一類是直接法，其特點是只包含有限次的四則運算，在每次運算都無舍入誤差的情況下，所得到的是方程組的準確解。由于實際計算中總是有舍入誤差，所以實際得到的也是近似解。令一類是迭代法，它首先選擇一組初始值，再運用同樣的計算步驟，重復(fù)計算，得到近似解。由于這類方法中出現(xiàn)了極限過程，必須研究迭代過程的收斂性。本章主要介紹：直接法中的高斯消去法和主元高斯消去法。迭代法中的簡單迭代法和塞德爾單迭代法。3.2.1順序高斯消去法以n=4為例說明高斯消去法的計算過程，設(shè)有線性方程組anxia2ixia31XIaanxia2ixia31XIa川X]+a12x2+a13x3+a14x4=bj+ +a+a24X4=b2~~=>+a32x2+a33x3+a34x4=b3+a42x2+a43x3+a44x4=b4a11x,+a12x2+a13x3+a14x4=b,a，*⑴32⑴42+a^Xj+&22X4=b?+a；?X3+a；?X4=b,=+a；?X3+a；：,=b?a“X|+al2x2+a13x3+a14x4=b.aux】+a12x2+a13x3+a14x4=b,a^x2a^x2+a^x3+a^x4=b^a22）X2+a22）X3+a22）X4=a*+a（2）d33X3a*+a（2）d33X3十d34分（2）乂+分（2）d43X3十d44X4=b;a數(shù)3+a數(shù)4=中a^x4=b￡）經(jīng)過3次消元步驟，得到以上形式。從最后一個方程中解出山依次回代得到方程組的全部解。計算過程見教材(P37)算法3.2』3.2.23.2.2主元消去法例3.2.2解線性方程組《0.0003X,+3.0000x2=2.0001例3.2.2解線性方程組《1.0000%,+1.0000x2=1.0000解：按4位小數(shù)計算，得至iJxi=0.6667,孫=0.0000, 準確解(勺=1/3,*2=2/3),誤差很大。如果將方程組的順序?qū)Q，得到《LOOOOX]+LOO。。％=L0000

0.0003^+3.0000x2=2.0001如果將方程組的順序?qū)Q，得到《按4位小數(shù)計算，得到占=0.3334,*2=0.6667,準確解(勺=1/3,*2=2/3),誤差很小。在順序主元消去法中，如果遇到方程組Ax斗中4的主對角線元素為零時，計算過程不能進行。如果A的主對角線元素的絕對值很小，也會產(chǎn)生較大的舍人誤差，使得最后得到的解可準確解有較大的誤差。工程中使用較多的是列主元高斯消去法，計算過程見教材(P42)算法3.4(6xi+3x2+2x3=6例2.2.3用順序消去法解方程組〈10X1+5X2+6X3=0I8xi+5x2+3x3=0方程組的增廣矩陣[A|b]63 26105 6085 30消元6.00003.00002.00006.0000002.6667-10.000001.00000.3333-8.0000方程組系數(shù)矩陣主對角線元素為零，消元過程無法進行!例2.2.4用主元高斯消去法解例2.2.3中的方程組。方程組的增廣矩陣[A|b]6326105608530選主元1056063268530消元10.00005.00006.0000000-1.60006.000001.0000-1.80000選主元10.00005.00006.0000001.0000-1.8000000-1.60006.0000消元10.00005.00006.0000001.0000-1.8000000-1.60006.0000回代得到方程組的解5.6250 -6.7500 -3.7500定理3.3.1設(shè)矩陣4的各階順序主子式不等于O"0(k=l,2,…n),則A有唯一的LU分解，A=LU,其中L為單位下三角矩陣，。為上三角矩陣。對于線性方程組Ax=b如果已經(jīng)得到了的LU分解，方程組的解可以通過以下方法得到LUx=b，令I(lǐng)7x=y，則有D="所以，求解Ax=6如的問題等價于求解兩的系數(shù)矩陣為三角方陣的方程組。X=4k-1憶=2,3,…〃y=lXn=yn/UnnX*=(”-Z k=n-\,n-2---\六2+1求解方程組AxR的計算變得非常簡單。例3.3.1解線性方程組例3.3.1解線性方程組A尸力PA=LU.LUx=Pb,Ly=Pb,Ux=yA=magic(3),b=[l;2;3],[L,U,P]=lu(A),y=L\(P*b),x=U\y對稱矩陣的LDI?分解設(shè)矩陣4的各階順序主子式不等于。卜和(k=l,2,…n),則A有唯一的LV分解，其中L為單.位下三角矩陣，U為主對角線元素不為0的上三角矩陣。M11W12 .0?.0]1u=0〃22 ?0“22 ,-0°1?“"2"/M22=DV00.00.■:°000,,Unn.00?unnJ［。0 ??1從而A=L。匕若為對稱矩陣A=〃，則有4=〃=(LOV)T=戶(。①，其中為L為單位下三角矩陣，(。〃)為主對角元素不為0的上三角矩陣，由LU分解的唯一性得到，L=L,于是定理3.4.1設(shè)對稱矩陣4的各階順序主子式不等于A#。(k=l,2,...n),則唯一確定一個單位下三角矩陣L和主對角元素不為0的對角矩陣。，使例3.4.1A=[5105-5;1024-2-6;5-244 -6-1123],[L,P]=chol(A),L'*L如果線性方程組Ax=b系數(shù)矩陣是對稱正定的，A=lJL,lJlx=b，令Lx=y,LTy=b例3.4.3b=[510-1-19],,y=L'\b,x=L\y線性方程組解的可靠性向量范數(shù)和誤差向量定義3.5.1設(shè)是定義在W上的實值函數(shù)，如果對于W中的任意向量及任意非負實數(shù)A,滿足⑴（非負性）對任意向量x都有|閨|>0且想|=0的充分必要條件是尸0。⑵（正齊性）對任意實數(shù)A和任意向量x都旬"|曰|國|。（3）（三角不等式）對任意向量版旬歸|r||+|刖稱||.||為向量范數(shù)。定理3.5.1對瞪上的任意向量x=（xg…x“）T,記IM=同+上|+…聞Wl2="：+￥+.?一|x|L=max|xj||x||p||x||2,||Xb都是向量范數(shù)。定義352設(shè)II.||是定義在/?而上的實值函數(shù)，如果對于Rg中的任意矩陣4s及任意非負實數(shù)%滿足（1）（非負性）對"任意矩陣A都有且M||=0的充分必要條件是4=0。（2）（正齊性）對任意實數(shù)k和任意矩陣4都有||必||=川川|。（3）（三角不等式）對任意矩陣4,5都有M+用區(qū)MH+II訓(xùn)，||4陽區(qū)||4||旭||稱||.||為矩陣范數(shù)。定理3.5.2對Rmn上的任意矩陣A記1閥=,黑 ||<=maxJ|a..|Ji=l j=lIIAIh，||A||2,IIAIL都是矩陣范數(shù)。,貝”|A卜臀刈Ax|是一種矩陣范數(shù)，并且它與給定的向量范數(shù)相容。定理3.5.3對N,貝”|A卜臀刈Ax|是一種矩陣范數(shù)，并且它與給定的向量范數(shù)相容。A=[-5-3l;401;41-2],x=[2-51]'norm(x,l),norm(x,2),norm(x,inf)nonn(A,l),nonn(A,2),norm(A,inf)norm(A*x,l),norm(A*x,2),norm(A*x,inf)計算各種范數(shù)。定義3.5.2設(shè)是定義在R2上的范數(shù)，對于火的中的任意矩陣A稱cond(A月|A||||A4||為矩陣A的條件數(shù)。例3.5.3A=[11.0001;11],b=[2;2],bl=[2.0001;2],A\b,A\bl當一個線性方程組AxW的系數(shù)矩陣4的條件數(shù)很大時，稱方程組是病態(tài)的，判斷線性方程組Ax2是否為病態(tài)的方法是：1、當|det(A)|很小或A的一些行或列近似相關(guān)時，方程組可能是病態(tài)的；2、使用列主元消去法，出現(xiàn)主元的絕對值很小時，方程組可能是病態(tài)的；3、當A的元素在數(shù)量級上有很大差別，且無一定規(guī)律時，方程組可能是病態(tài)的。對于病態(tài)的方程組可以采取下面的方法：1、采用高精度的算術(shù)運算，例如采用雙精度數(shù)運算。2、當A的元素在數(shù)量級上差別很大時，采用行平衡或列平衡的方法可以降低A的條件數(shù)。A=[l1;1le4],b=[2;le4],cond(A,inf),s=max(A),d=diag(1./s),cond(d*A,inf)clear,n=4; %l[Zn=6,8,10等進行比較b=ones(n,1)x=hilb(n)\bxO=invhilb(n)*bmax(abs(x-xO))cond(hilb(n))rats(hilb(4))%字符串rats(invhilb(4))%字符串設(shè)有方程組AxH，變?yōu)榈问?，x=Mx+f,或x(k+i)=Mx")+f,任取初始值x(°)程迭代得到￡0) ?(i)1r⑵Xfx>Ax3￡0) ?(i)1r⑵Xfx>A以n=4為例XT”廠Y(k+1)X2Y(k+1)A3Y(k+DLA4」=m21m3i_m4im12m22m32m421%m?3m43m；m24m34m人工廠Y(k)X2xf)[x]+■f.f2f3J4.TT?Tx(k+l)yx(k)f寫成分量形式(\(k+l)一mx(k)X1 -milXl+m摟x.+m13X3k)+m14x*lk>+f.X并)iX,+m22X2k)+m23x<k)+m24x(4k)+f2x<k+l)=m31x{k}+m32X2k)+m33x<k)+H134X，+0xa)mx；k)+m42X2k)+II143X，+m44x^)+f4定理3.6.1若存在某種范數(shù)||.||,使|| 則簡單迭代法收斂。定理362迭代公式x"T)=Mx(k)+f,對任意初始值x(°)和f都收斂的充分必要條件是矩陣M的各個特征值的模都小于1。在簡單迭代法的基礎(chǔ)上作改進x(k+l)=M,x(k+1)+M2x(k)+f,以n=4為例X]1）A2Y（^+1）y（H1）_A4 _二-0m21m3l_m4i00m32m42000m21o-000Yd）一Yd）A2W+l）W+1）LX4.+00Omi2mi3m22m230 m330 0mjm24m34A1x"）A2y（K）A3X⑹LX4.+-f;f?f3.f4.TTTTTTx（k+1）M,x(k+】)m2x（k）f§3.7雅可比迭代法與高斯-塞德爾迭代法雅可比迭代法設(shè)方程組Ax斗，的系數(shù)矩陣主對角線元素不為零，A可以分解為A=O+L+U,其中O,L,U分別取A的主對角線，下三角和上三角部分的元素。^\(D+L+U)x=b,0*=a+。比+瓦*=0“億+。比+。%,由此可以得到迭代形式x<k+1)=-Dl(L+t7)r(k)+D'b (k=0,l,2…) (3-67)雅可比迭代法的迭代矩陣M}=-D\L+U).寫成分量形式 %*+4)/?！?尤;'+優(yōu))/a22(3-68)靖+"=(-。|2球,%2*+，>=(~a2\X %*+4)/。“ 尤;'+優(yōu))/a22(3-68)產(chǎn)/，=(-。.2以-。"3靖 外”-出+:)/ann高斯-塞德爾迭代法在計算(3-68)式時，用已經(jīng)計算出的與進行后續(xù)的計算，可以得到xf+"= )—。13衛(wèi)" a\nX^}+4)/。11(3-69)E*+l)=(_%x尸＞一。23石" 。2況*'+2)/a22(3-69)x”=(r2鏟)-心鏟)-…-*咸?+么)/%稱式(3-69)為高斯-塞德爾迭代法。式(3-69)的矩陣表達形式為x(k+i)=D\Lx*+i)_n'U)xw+D'b (k=0,l,2…)x(k+1)="(D+L)-'Uxw+(D+L)'b (k=0,l,2…) (3-70)高斯-塞德爾迭代法的迭代矩陣mgs==-(d+ly'u.雅可比迭代法和高斯-塞德爾迭代法的收斂性定理3.7.1雅可比迭代法和高斯-塞德爾迭代法收斂的充分必要條件是它們的迭代矩陣的各個特征值的模都小于1.定理3.7.2若存在某種范數(shù)||.||,使則雅可比迭代法收斂。若存在某種范數(shù)使||Mgs||<1,則高斯-塞德爾迭代法收斂。例3.7.1分別用雅可比迭代法和塞德爾迭代法解方程組1Ox1—￡—2工3=72<-%1+10x2-2x3=83誤差e<10-6-%)一9+513—42雅可比迭代法迭代公式X(k+D=D"(L+次(?%D-\h高斯?塞德爾迭代法迭代公式x(k+,)=-(D+L)-lt/x(k)+(D+L)?A=[1O-1-2;-l10-2;-l-15]b=[728342"D=diag(diag(A));U=triu(A,l);A=[1O-1-2;-l10-2;-l-15]b=[728342"D=diag(diag(A));U=triu(A,l);L=tril(A,-l);M=-inv(D)*(L+U);b=inv(D)*b;x=b;fork=l:l0,x=M*x+b,endA=[10-1-2;-l10-2;-l-15]b=[728342了；D=diag(diag(A));U=triu(A,l);L=tril(A,-l);M=-inv(D+L)*U;b=inv(D+L)*b;x=b;fork=l:5,x=M*x+b,end例3.7.2考察用雅可比迭代法和塞德爾迭代法解方程組的收斂性。雅可比迭代法迭代公式x<k+,)=-D-'(L+(Or例3.7.2考察用雅可比迭代法和塞德爾迭代法解方程組的收斂性。雅可比迭代法迭代公式x<k+,)=-D-'(L+(Or(k)+D'h高斯-塞德爾迭代法迭代公式x<k+1)=-(D+L)-1t7x(k)+(D+L)'bA=[2-11;11l;l1-2]D=diag(diag(A));U=triu(A,l);L=tril(A,-l);MJ=-inv(D)*(L+U);MS=-inv(D+L)*U;abs(eig(MJ))abs(eig(MS))例3.7.3用雅可比迭代法解方程組雅可比迭代法迭代公式x(M)=-D'(L+U)xw+D'bA=[l-22;-l1-l;-2-21]b=[l0-2"D=diag(diag(A));U=triu(A,l);L=tril(A,-l);M=?inv(D)*(L+U);b=inv(D)*b;x=b;x=M*x+b本章小結(jié)本章討論了解線性方程組的直接解法和迭代解法。直接解法比較適用與系數(shù)矩陣稠密（既零元素較少）的中、小型線性方程組，但對系數(shù)矩陣是帶狀或近似帶狀的大型線性方程組也適用。直接解法中的列上元高斯消去法具有精度較高和省時的優(yōu)點，是計算機中常用的算法。迭代解法中主要介紹了雅可比迭代法、高斯-塞德爾迭代法和松弛法。迭代法具有計算公式簡單、程序設(shè)計容易、占用計算機內(nèi)存較少的優(yōu)點。適用于解大型稀疏矩陣（既零元素較多）線性方程組。高斯-塞德爾迭代法是在雅可比迭代法的基礎(chǔ)上改進得到，在很多情況下可以加快收斂速度，但它的收斂域與雅可比迭代法不同，因此不能互相取代。松弛法可以加速迭代過程的收斂速度，但要適當選擇松弛因子（（KM2）。在選擇迭代法時，要特別注意檢驗方法的收斂性問題。第4章矩陣特征值與特征向量的計算§4.1引言求矩陣的特征值和特征向量，是代數(shù)計算中的重要問題。在自然科學(xué)和工程中的許多問題，例如電磁振蕩、橋梁的振動，機械振動等都可以歸結(jié)為求矩陣的特征值和特征向量問題。矩陣A的特征值和特征向量是指，如果數(shù)入和非零列向量x滿足關(guān)系式Ax=M,則數(shù)入稱為A的特征值非零列向量x稱為A的與特征值人對應(yīng)的特征向量。計算n階矩陣A的特征值，就是求特征方程|A-Q|=O的根%(i=l,2,...,n)=齊次線性方程組(A-XiI)x=O的非零解即是無對應(yīng)的特征向量。本章討論一些在計算機上計算矩陣的特征值和特征向量的較為穩(wěn)定的數(shù)值算法。1、’幕法：計算n階矩陣A的模最大的特征值（主特征值）及對應(yīng)的特征向量。任取n維列向量x（°）,用迭代公式xMLax00計算得到x（°）,x（1）,P，…設(shè)x<0）=aiV|4-a2V2+…+anvn,因為AVj=%Vj所以x（l）=Ax""=a仇］v）+a2九2V2+...+anX.nvnx（2>=Ax"，=aili2Vi+a2九22V2+…+anXn2vn一般地有x"T）=Ax'"=ai?i1kV1+a2X2kv2+...+anXnkvn=Xik［aivi+a2（X2/X|）kv2+...+an（Xn/A，i）kvn］當k充分大時x*+i）=a仇J+/產(chǎn)儲W）向量x*+i）與x,k）向近似地只差一個倍數(shù)，這個倍數(shù)就是模最大的特征值儲。1,反'募法：Ax=?ix,A-'Ax=A-'Xx,A-'x=X-'x,即A的特征值的倒數(shù)M是A的逆矩陣A」的特征值。用事法求A-1的模最大的特征值，它的倒數(shù)就是A的模最小的特征值。對于2階方陣ut（0）au（g）=A=aHcosy（a2_4aH 「cos。sin01丁 「10[a：, U（e）=Lm。cosoj口（。川⑼=[030+a22sin20+a21sin20—（a22-au）sin20+a21cos20aHsin20+a22cos20-a21sin202-aiJsinZe+a2cos20令3（%2-a11）sin20+a2Icos20=0tan20=2a2]an-a22uT（e）Au（e）=入o'o入2a1la21a31-對于n階方陣（以n=3為例）A=a21a22a32a31COS0a32a33sin0c一X0X令tan20=all-J作變換矩陣a22u（e）二-sin0cos00 001則有uT（e）Au（e）=0XXXXXcos00sin9xx0令tan2。=——作變換矩陣U(0)=0 1 0則有Ut(9)AU(0)=XXXa”一a33-sin00cosO0xx1 0 0XXX令tan2。= 作變換矩陣U(0)=0cos0sin0則有uT(e)AU(e)=xx0a22一a330-sin0cos0x0x2a..一般地說，令tan29=——可以將A中的元素的ag和細變?yōu)?。在實際計算中采用以下公式aii-ajjsin0=.f2(1+4一心)cos0=>/l-sin20-1例4.3.1用雅可比求對稱矩陣A=§4.4QR方法*-I-I-1的特征值和特征向量。（見教材P86）QR方法是求一般矩陣A的全部特征值和特征向量的?種迭代方法。其基本思路是利用矩陣A的QR分解,其中，R是一個上三角矩陣，Q是正交矩陣。通過迭代格式A=QkRkAk1將A化為相似的上三角矩陣（或213例4.4.1用QR方法求矩陣A=123的特征值。012分塊上三角矩陣），從而求出A的全部特征值和特征向量。-0.6215-0.62150.47700.6972-0.6760-0.6760-0.29354.3028特征向量-0.8165-0.40820.4082特征值1.0000a=[213;123;012],[q,r]=qr(a),a=r*q本章小結(jié)本章介紹了求矩陣的特征值和對應(yīng)的特征向量的幾種方法。鼎法可以求出矩陣的主特征值和時應(yīng)的特征向量，優(yōu)點是算法簡單，但當\XA2\=1時，收斂速度很慢。反惠法可以求出矩陣的模最小的特征值和對應(yīng)的特征向量。雅可比方法是利用一系列正交相似變換（即平面旋轉(zhuǎn)變換）把實對稱矩陣A化為對角陣（近似），從而求出實對?稱矩陣全部特征值。QR方法是用鏡向反射陣將矩陣A作QR分解，是一種求矩陣的全部特征值的有效方法。第5章插值與擬合§5.1引言己知表格函數(shù)y=f(x)XiXoX]X2Xn-1Xnf(Xi)yoyiY2yn-iYn構(gòu)造一個公式p(x)近似地表示f(x),解決這個問題的方法有兩類：一類是插值法，另一類是擬合法，又稱為逼近法。已知函數(shù)y=f(x)在互異點xo,X,,x2,xn.hxn上的函數(shù)值外,yi,y*yn，構(gòu)造一個函數(shù)p(x)使得p(Xi)=y這樣的問題稱為插值問題。y=f(x)稱為被插值函數(shù)，［Xoxn］稱為插值區(qū)間，p(x)稱為插值函數(shù)，Xo,X1,X2,…,X?l,Xn稱為插值點，在插值區(qū)間內(nèi)部用p(x)代替f(x)稱為內(nèi)插，在插值區(qū)間外部用p(x)代替f(x)稱為外推，R(x)=f(x)-p(x)稱為插值函數(shù)p(x)的誤差。定理：給出n+1個插值點及函數(shù)值XiXoX1X2Xn.1Xnf(Xi)yoyiY2yn-iyn求一個n次多項式pn(x)=ao+aix+a2X2+...+anxn1、給出2個插值點(xo,y0),(xi,y?？梢缘玫揭淮味囗検絧1(x)= -y0+ 里y1x0-X| x,-x02、給出3個插值點(xo,y°),(xi,yi)，僅2,丫2)可以得到二次多項式(x)=(x-xJCx-x^?(x-x0)(x-x2) 1(x-x0)(x-x,)2 (X0-X,)(X0-X2)° (X1-xo)(x,-x2)'(X2-X0)(X2-xj：不難驗證P2(x)滿足插值條件P2(Xo)=yop2(X1)=yip2(x2)=y23、給出n+1個插值點(x(),yo),(x〕，y]),…,y?)可以得到一個n次多項式Pn(x)=L0(x)y0+L,(x)y,+---+Ln(x)yn其中Lk(x)=-~~^4;_r?n(x)=(x-x0)(x-x1)---(x-xn)(x-xK)con(xK)例5.3.1按下列表格求y(-0.5)和y(0.5)的值。xl1 2 3y|71423解：插值多項式L,(x)=(x_x」(x_X2)方+(x_x())(x_X2)%+(x_x°)(x_X|)2 (X。-X|)(Xo-X2)O(x^XoXXj-xJ'(Xj-XqXXj-xJ=(x-2)(x-3)7+(x-l)(x-3)u+(x-l)(x-2)23"(l-2)(l-3) (2-l)(2-3) (3-l)(3-2)=x2+4x+2L2(-0.5)=0.2500 L2(0.5)=4.2500例4.3.2按下列表格求y(2.5)、y(4.5)、y(5.5)的值。xl1234 5y|0215 3解:插值多項式L4(x)=-0.7917x4+9.25x3-37.21x2+60.75x-32L4(2.5)=0.9180,L4(4.5)=6.1323,L4(5.5)=-8.9637§5.4均差插值公式已知函數(shù)Rx)在互異點Xo,X1, Xn上的值為Rxo),Rx1),f(xn)稱f(Xi，Xj)=，^—-(jHi)為函數(shù)Hx)在點Xi,Xj處的一階均差。Xj-Xjf(xj,Xk)-f(X;,Xi)稱f(Xj,Xj,Xk)=—— —出")為函數(shù)電)在點*[,\),*|(處的二階均差。Xk-Xii般地稱f(Xo，Xi，…,Xn)=fg'X2'…'Xn)f(X0，Xl，…,Xn-l)為函數(shù)付在點x°,x】，…，X。上的n階xn-X0均差。牛頓插值公式Nn(x產(chǎn)f(x0)+(x-x0)耳xo,Xi)+(x-Xo)(x-Xi)f(Xo,X1,X2)+...+(x-X0)(X-Xi)…(X-Xn-l)f(XQ ,Xn)牛頓插值公式具有遞推關(guān)系Nk+i(x尸Nk(x)+(x-Xo)(x-X])...(x-xk)f(xo,Xi,…，Xk+i)新增加一個節(jié)點，只需要增加計算一項(X-Xo)(X-Xi3..(X-Xk)f(Xo,X|”..,Xk+l)1、差分已知函數(shù)Rx)在等距節(jié)點Xo,X1,…，Xn上的值XXox+hx+2hx+nhf(x)yoyiyaYn稱△yk=yk+i-yk為函數(shù)Rx)在點xk處的一階差分。稱yk=Z^yk+i-Aiyk為函數(shù)f(x)在點Xk處的二階差分。一般地稱yk=^Am'y^i-A1"'yk為函數(shù)f(x)在點Xk處的m階差分。牛頓向前插值公式. tt(t-1)a2 t(t—1)…(t—n+1)AnTOC\o"1-5"\h\zNn(xo+th)=yo+-Ay0+---A~y0+???+ : △Yo1! 2! n!牛頓向后插值公式、r/ ［、t口t(t+l)v?2 t(t+1)…(t+n—1)7nNn(x0+th)=y0+-Vy0+---Vy0+???+ : VYo1! 2! n!定理：給出n+1個插值點上的函數(shù)值和導(dǎo)數(shù)值XiXoXlX2Xn.1Xnf(Xi)yoyiY2yn-iynf'(Xi)y。'y/y2'yn-izy；求一個2n+l次多項式H(x)滿足Hn(Xi尸y,Hzn(Xi)=y\(i=0,l,2,...,n)oh(x)=zyiAi(X)+Zy：Bi(x)i=0 i=0其中Ai(x)=[1—2(x—Xj)玄———rBi(x)=——fj=Oxi-xj(x-Xj)?(Xi)] (x-Xj)?(Xi)「j=i5.7.1龍格現(xiàn)象例5.7.1給定函數(shù)/(x)=―J-1<X<1,取等距節(jié)點，構(gòu)造f(x)的n次拉格朗日插值多項式，并l+25x畫出圖形。(n=4,10)fMnline('l./(l+25*x.A2)‘)x=-l:0.5:l;a=polyfit(x,f^x),4);xl=-l:0.2:l,al=polyfit(x 1),10);xi=-l:0.01:l;plot(xi,f(xi),xi,polyval(a,xi),xi,polyval(al,xi))gridon,legend('f(x)7n=4Vn=10,)在xoy平面上給定n+1個點(x0,y0),(xhyi),...,(xn,yn)＞構(gòu)造一個函數(shù)S(x)滿足以下條件(1)S(xj)=yi(i=0,l,2,...,n)(2)在區(qū)間(x0,xn)內(nèi)S(x)具有連續(xù)的二階導(dǎo)數(shù)(3)在每個子區(qū)間［xi,xj上S(x)(表達式是Si(x))是一個三次多項式。滿足以上條件S(x)稱為三次樣條插值多項式。三次樣條插值多項式的推導(dǎo)設(shè)在子區(qū)間［Xi_bXi］上S(x)=Si(x)(i=1,2，…,n)由條件(1)得SGi一Ly1,S((xj)=yi設(shè)S(x)在節(jié)點x-i處的二階導(dǎo)數(shù)為Mi，在節(jié)點Xi處的二階導(dǎo)數(shù)為Mi，即S\(xz尸Mi，S”i(Xi尸Mi，顯然S'(x)是x的線性函數(shù)，根據(jù)拉格朗日插值公式有S；(x)=?二i2泊+x-Xi』Mj,記xi-xi=hi有S：(x)==^Mi_i+=^」MiXi-i-Xj x,-x,.1 -hj %將上式積分兩次得到 Sj(x)=M,.1a：?)+Mj+C1X+C2-6hj 6hj利用Si(xi")=yz,Si(xi尸yi定出積分常數(shù)C］和C?Mi-7-+C1xMi-7-+C1xi_1+C2=Yj-j6 解此方程組得到Mj-y-d-CjXj+C2=Yj6C2Yi-Yi-i+?(Mt-MJ

oh,代入上式整理后得到3 3+ +gr+s；Jm,h,,s；Jm,h,, ,,+U(Mi-MJ

oS；(xJ=M奇+Yi-Yi-ihih+U(M1-MJ6類似地有S：+1(x)S：+1(x)=Mi(X-Xj+])-2hi+I+Mj(x-xj31yi+l-yj

2hj+iLiS；M(xJ=Mj與+白?工+零(M「MQ-2hi+I6因為S'i+i(Xi尸S'“Xi)所以有Jkm,,+/乂+h±LMi+l=Yi+I_yi_yi-yi-'整理后得到關(guān)于位知數(shù)6 3 ' 6 1+1hi+1hiMo,Mi,...,Mn 的線性方程組為Mw+2Mi+biMi+i=&(i=l,2,...,n-l)其中a.=-%—b= %_.』=」一必口_立組)'瓦+hi+i' %+hi+| 'hi+%+ihi+I h邊界條件：n-1個方程組，n+1個位知數(shù)，所有需要補充兩個條件，常見的是以下兩種(1)給定端點的二階導(dǎo)數(shù)Mo=a,Mn=b,特別地，當a=b=O時稱為自然樣條。(2)給定端點的一階導(dǎo)數(shù)S'(xo)=a,S/(xn)=b,這時有2M0+M,=—(Y1-y°-a) Mn, +2Mn=—v1i i ' n—i niv i z%hi hn hn(2)M(2)Mo=1,M3=O例5.7.1給定插值條件，和兩種邊界條件(1)項=1,013=0xIO12 3y|000 0分別求出兩種條件下的三次樣條插值多項式的表達式。m3=0sl(x)=0.73333x3-1.7333xM0=l,M3=0sl(x)=-0.21111xM0=l,M3=0sl(x)=-0.21111x3+0.5x2-0.28889x0<x<ls2(x)=0.055556x3-0.3x2+0.511llx-0.26667l<x<2s3(x)=-0.011111x3+0.1x2-0.28889x+0.266672<x<3s2(x)=-0.2x3+1.0667x2-1.8x+0.93333l<x<2s3(x)=0.066667x-0.53333x2+1.4x-1.22<x<3§5.9曲線擬合的最小二乘法給出表格函數(shù)XiXoX1X2Xn.1Xnf(Xi)yoyiY2yn-iyn求一個m(m<n)次多項式為Pm(x)=a()+aix+a2x2+…+amX111逼近這組數(shù)據(jù)。pm(x)的系數(shù)的確定:aoao+a1xo+a2Xo+---+amXo-y0=R()a0+a,xl+a2x?+-+ainx?-y,=R]即gajX；-yj=Rj(i=0,1,2,...,n)j=oa0+a山”+a?x：+…+amX：—y”=Rnn nmO(ao，a”???am)=Z&2=￡(ZajX：-丫。2由微分學(xué)知，使°(即,ai,am)達到極小值的即，a1，…,i=0 i=0j=0am滿足必要條件 .(a。：!?…a/=2次(￡ajX,一yjx：=0(k=0,1,2,m)^aK i=0j=0寫成分量形式%(〃+1)+為E若+出EX；+…+4Z町=2yi=0 i=0 /=0 i=0。。2＞，+%?/1＞，3+"”，“2＞產(chǎn)=2＞，凹 -于柳石超如1 i=0 i=0 i=0 i=0 i=0 .止現(xiàn)力在組41＜+ ＜'+a2tx/2+…+a，5丁=t＜'X、 i=0 i=0 i=0 i=0 i=0解此方程組得到徹，ai,…,am即得到所要的多項式。例5.8.1有數(shù)據(jù)表如下，分別用二次和三次多項式逼近這組數(shù)據(jù)。x1-3 -2 -1 0 1 23y|1 0 0 0 0 12P1（x）=0.1786x4-0.5714P2（x）=0.17857x2+0.17857x-0.14286p3（x）=3.2895*10l7x3+0.17857x2+0.17857x-0.14286一次多項式逼近誤差平方和2.8214,二次多項式逼近誤差平方和0.1429,三次多項式逼近誤差平方和0.1429o本章小結(jié)本章介紹了兩種構(gòu)造函數(shù)Rx）的逼近函數(shù)的方法一插值法和曲線擬合。關(guān)于插值法，討論的是多項式插值，它要求所構(gòu)造的多項式函數(shù)嚴格地通過給定的所有數(shù)據(jù)點。由Lagrange插值基函數(shù)的討論，導(dǎo)出了Lagrange插值公式及其余項公式。作為對Lagrange插值的改進，建立了具有遞推性的牛頓插值公式。Hermite插值是一種在插值節(jié)點上插值函數(shù)與被插值函數(shù)不僅有相同的函數(shù)值，而且還有相同的一階導(dǎo)數(shù)值的多項式插值。由于高次插值會產(chǎn)生龍格現(xiàn)象，因此討論了分段插值法，重點介紹了三次樣條插值公式。關(guān)于曲線擬合，討論的是曲線擬合的最小二乘法，它不要求所構(gòu)造的逼近函數(shù)嚴格地通過給定的所有數(shù)據(jù)點，只是在多項式中，以使得殘差的平方和最小為標準，選擇多項式，以作為被逼近函數(shù)的近似替代。本章最后介紹了數(shù)值微分，其中5點公式是常用的。第6章數(shù)值積分和數(shù)值微分§6.1引言b計算定積分的牛頓-萊布尼茲公式Jf(x)dx=F(b)-F(a),其中F(x)是函數(shù)[x)的原函數(shù)。在實際應(yīng)用中a(1)常遇到某些函數(shù)f(x)的原函數(shù)不能用初等函數(shù)表示，例如sin*),cos(x2),sin(x)/x,l/log(x)等。(2)有些函數(shù)fi[x)是用表格形式表示，無法得到它們的原函數(shù)。(3)有些函數(shù)f(x)的原函數(shù)十分復(fù)雜，不利于工程上的使用。因此，要研究計算定積分的近似方法：數(shù)值積分。b, n定義6.1.1若一個求積公式=Z例〃七)，對于任何次數(shù)不超過m的多項式，等式都準確成立,a=0而對于某一個m+l次多項式不能精確成立，則稱這個求積公式的代數(shù)精確度為m。牛頓一科特斯型積分公式由上一章知，任意函數(shù)f(x)可以用一個拉格朗日插值多項式6(x)近似表示，因此Rx)的定積分也可以用TOC\o"1-5"\h\zb bd(x)的定積分近似表示Jf(x)dx=J①“(x)dx以此為基礎(chǔ)得到的數(shù)值積分公式稱為牛頓?科特斯型積分公式。b b bJf(x)dx=j①n(x)dx=jL0(x)y0+L1(x)y,+???+Ln(x)yndxa a ab b b=y0fLo(x)dx+yifLi(x)dx+--.+ynjLn(x)dx=A()y0+A,yl+…十人一口a a a其中 Lk(x)=--(on(x)=(x-x0)(x-xl)---(x-xn)Ak=[LK(x)dx(x-xK)(on(x) ，MXa=x0<X|<...<xn=b為一組等距節(jié)點，☆x=a+bt,h=(b?a)/n得到2 (0(x) f_n(n-K)knfAK=fLK(x)dx=f-——\ = ^ft(t-l)^(t-K+l)(t-K-l)-(t-n)dt； i(x-XK)0)n(x) K!(n-K)!,a a U(_])(n-K)(_])(n-K)

nK!(n-K)!AK=(b-a)C^Jt(t-1),,?(1—K.+1)(t—K—1),,,(t—n)dx,于是AK=(b-a)C^0b所以得jf(x)dx=C<n)y0+C*n)yi+-+C<n,yn對不同的n系數(shù)如下n牛頓一科特斯型積分公式系數(shù)余項（誤差）11/2,1/2 （梯形公式）-(1/12)h3f⑵⑴) a<T|<b21/6,4/6,1/6 （拋物線公式或辛卜生公式）-(1/90)h5f<4)(T])31/8,3/8,3/8,1/8-(3/80)11”⑷⑺47/90,32/90,⑵90,32/90,7/90 (柯特斯公式)-(8/945)h7f⑹(r|)519/288,75/288,50/288,50/288,75/288,19/288-(275/12096)h7f⑹⑴)641/840,216/840,27/840,272/840,27/840,216/840,41/840-(9/1400)h3f⑵⑴)復(fù)合求積公式1、復(fù)合梯形公式取a=x（）vx|V??.VXn=b為一組等距節(jié)點將積分區(qū)間［ab］分為n個子區(qū)間［xkxk.）］" hXk-Xk-i=(b-a)/n=h,在每個子區(qū)間［xkXk」］上應(yīng)用梯形公式得：Jf(x)dx=—(yk_1+yk)Xi 2bh于是得到復(fù)合梯形公式 Jf(x)dx=-［y0+2(Y1+y2++yn_,)+yn］2、復(fù)合辛卜生公式在積分區(qū)間［ab］取點2=\0〈乂1〈?..〈乂211.1＜乂211=1＞將［ab］分為2n個子區(qū)間，令x2k-2-x2k=(b-a)/n=2h,在每個子區(qū)間［x2k.2x2k］上應(yīng)用辛卜生公式得：X2k hJf(x)dx=§(y2k_2+4Y2k-l+Y2k)(k=l,2,n)X2k-2于是得到復(fù)合辛卜生公式b 卜jf(x)dx=-［y0+4y,+2y2+4y3+2y4+---+2y2n_1+4y2n_,+y2n］例6.3.1f^inline(,exp(l./x),；x,);%2.0200586244339...a=l;b=2;n=5;h=(b-a)/n;yl=f(a:h:b);y2=f(a+h/2:h:b-h/2);S=(yl(l)-Fyl(n+l)+2*sum(yl(2:n))+4*sum(y2))*h/6例6.3.2f^inline(,exp(l./x),；x,);a=l;b=2;n=5;h=(b-a)/n;y=f(a:h:b);T=(y(l)+2*sum(y(2:n))+y(n+l))*h/23、復(fù)合柯特斯公式在積分區(qū)間［ab］取點a=XoVX］V…<X4n/vx4n小將［ab］分為4n個子區(qū)間，令x4k-4-x4k=(b-a)/n=4h,在每個子區(qū)間Dugx4k］上應(yīng)用柯特斯公式得：n)(k=l,2,…，n)(k=l,2,…，n)=面(7丫4k-4+32y4k一3+12y4kT+32y4k_1+7y4k) (k=l,2,9(Jx4k-4得到復(fù)合柯特斯公式b 4h J, _n, j, ?一!ff(x)dx?—[7y0+32^y4k_3+12^y4k_2+32^y4k_,+14gy4k+7y4n]a 9U k=l k=l k=l k=l4、步長h的自動選擇（1）復(fù)合梯形公式|T2n-Tn|<3e （e表示允許誤差）其中I,和T2n分別表示取n和2n時用復(fù)合梯形公式計算得到的積分近似值。（2）夏合辛卜生公式|S2「Sn|<15e （e表示允許誤差）其中Sn和s2n分別表示取n和2n時用復(fù)合辛卜生公式計算得到的積分近似值。（3）復(fù)合柯特斯公式|C2「Cn|<63e （e表示允許誤差）其中Cn和C2n分別表示取n和2n時用復(fù)合柯特斯公式計算得到的積分近似值。TOC\o"1-5"\h\z1、梯形公式與辛卜生公式之間的關(guān)系4 1 4 1 4 1= /S2= T4 12,…，S?k= T9k+i T)k4-124-11 24-144-12 2 4-12 4-122、辛卜生公式與柯特斯公式之間的關(guān)系42 1 42 1 42 1C.=—: S2z S.C2=-r S4zS2,---,CK=-r SK+I；SK142-1242-11 242-1442-12 2 42-12 42-123、龍貝格積分公式§6.5高斯求積公式高斯積分問題的提出適當選擇插值求積公式的節(jié)點和系數(shù)可以提高公式的代數(shù)精確度。例如構(gòu)造兩點公式1Jf(x)dx=CO0f(x0)+ ) (6-23)-1如果限定節(jié)點刖=-1/=1，那么得到的積分公式是1j/(x)dr?/(-1)+/(1)-1代數(shù)精確度僅為k如果對(6-23)節(jié)點和系數(shù)都不加限制,則可以使其代數(shù)精確度大于1.事實上，要使(6-23)對函數(shù)1,x,f,/,都精確成立，只要(6一23)中的節(jié)點和系數(shù)滿足%+q=2COqXq+ =0一，再二一

313=>%=4=lx一，再二一

313+四x：=0所以，有,(x)dx?/(-事)+/(^)，這個積分公式至少有3次代數(shù)精確度。TOC\o"1-5"\h\zb, n定理6.5.1形如J/(x)dx=Z3*/(x*)的插值型積分公式的代數(shù)精確度最高不超過2n-l次。a &=。6.5.2高斯積分公式b, n定義651使插值求積公式J/(x)必0j(xj的代數(shù)精確度為2n-l的節(jié)點片,孫…冬稱為高斯點,a &=。對應(yīng)的插值求積公式稱為高斯求積公式。定理6.5.2節(jié)點Xi,X2,??.dn是高斯點的充分必要條件是在區(qū)間乂ab]匕n次多項式ban(x)=(x-x\)(x?X2)??.,Q-Xn)與所有次數(shù)不超過n-1的多項式。(五)都正交。即J〃〃(x)Q(x)dx=O。a定理6.5.3:若X1,X2,…,Xn是高斯點，則以(x)=2(〃,)/“(x) 4(X)=']_ .n(2〃)！" " 2"〃！dx2Ln(x)稱為勒讓德多項式。=2=2*=|定理6.5.4高斯求積公式的系數(shù)恒為正，且有以=以上討論實際上給出了構(gòu)造高斯積分公式的方法，即只要找出區(qū)間［ab］上n次正交多項式的n個零點，對于區(qū)間［ab］可做線性變換，x= +"■，使積分區(qū)間變?yōu)椋?11］,所以對于高斯積分公式，只考慮積分2 2區(qū)間［-11］?下表給出了1-5個節(jié)點的高斯積分公式的節(jié)點和系數(shù)，可以看出節(jié)點和系數(shù)具有對稱性。n節(jié)點Xk⑺系數(shù)3k⑺余項（誤差）積分區(qū)間[-11]102f⑵（n）/3 -i<n<i2-0.5773503+0.577350311f（4）（n）/i353-0.57735030+0,57735035/98/95/9f<6）（Ti）/157504-0.8611363-0.3399810+0.3399810+0.86113630.34785480.65214520.65214520.3478548f⑹01）/34728755-0.9061799-0.53846930+0.5384693+0.90617990.23692690.47862870.56888890.47862870.23692691⑼01）/1237732650二重積分的數(shù)值積分法數(shù)值微分給出表格函數(shù)如下表所示，求f(X)在節(jié)點f(Xj)處的導(dǎo)數(shù)的近似值。XiXoX1X2Xn-lXnf(Xi)yoyiyzyn-1Yn方法：用插值多項式的導(dǎo)數(shù)近似fi[x)的導(dǎo)數(shù)。f'(Xo)=4(-3yo+4y|一丫2)2hf*(xo)=p-(yo-2y!+y2)三點公式：”f'(X|)=!(一y()+y2) -2hf*(x1)=p-(y0-2