數(shù)值計算教程

數(shù)值計算方法的對象與特點(diǎn)什么是數(shù)值計算方法現(xiàn)代的科學(xué)技術(shù)發(fā)展十分迅速，他們有一個共同的特點(diǎn)，就是都有大量的數(shù)據(jù)問題。比如，發(fā)射一顆探測宇宙奧秘的衛(wèi)星，從衛(wèi)星設(shè)計開始到發(fā)射、回收為止，科學(xué)家和工程技術(shù)人員、工人就要對衛(wèi)星的總體、部件進(jìn)行全面的設(shè)計和生產(chǎn)，要對選用的火箭進(jìn)行設(shè)計和生產(chǎn)，這里面就有許許多多的數(shù)據(jù)要進(jìn)行準(zhǔn)確的計算。發(fā)射和回收的時候，又有關(guān)于發(fā)射角度、軌道、遙控、回收下落角度等等需要進(jìn)行精確的計算。有如，在高能加速器里進(jìn)行高能物理試驗，研究具有很高能量的基本粒子的性質(zhì)、它們之間的相互作用和轉(zhuǎn)化規(guī)律，這里面也有大量的數(shù)據(jù)計算問題。計算問題可以數(shù)是現(xiàn)代社會各個領(lǐng)域普遍存在的共同問題，工業(yè)、農(nóng)業(yè)、交通運(yùn)輸、醫(yī)療衛(wèi)生、文化教育等等，各行各業(yè)都有許多數(shù)據(jù)需要計算，通過數(shù)據(jù)分析，以便掌握事物發(fā)展的規(guī)律。研究計算問題的解決方法和有關(guān)數(shù)學(xué)理論問題的一門學(xué)科就叫做計算方法。計算方法屬于應(yīng)用數(shù)學(xué)的范疇，它主要研究有關(guān)的數(shù)學(xué)和邏輯問題怎樣由計算機(jī)加以有效解決。數(shù)值計算方法的內(nèi)容數(shù)值計算方法也叫做計算數(shù)學(xué)或數(shù)值分析。數(shù)值計算方法主要內(nèi)容包括非線性方程求根、線性代數(shù)方程組解法、微分方程的數(shù)值解法、插值問題、函數(shù)的數(shù)值逼近問題、概率統(tǒng)計計算問題等等，還要研究解的存在性、惟一性、收斂性和誤差分析等理論問題。我們知道五次及五次以上的代數(shù)方程不存在求根公式，因此，要求出五次以上的高次代數(shù)方程的解，一般只能求它的近似解，求近似解的方法就是數(shù)值分析的方法。對于一般的超越方程，如對數(shù)方程、三角方程等等也只能采用數(shù)值分析的辦法。怎樣找出比較簡潔、誤差比較小、花費(fèi)時間比較少的計算方法是數(shù)值分析的主要課題。在求解方程的辦法中，常用的辦法之一是迭代法，也叫做逐次逼近法。迭代法的計算是比較簡單的，是比較容易進(jìn)行的。迭代法還可以用來求解線性方程組的解。求方程組的近似解也要選擇適當(dāng)?shù)牡?使得收斂速度快，近似誤差小。在線性代數(shù)方程組的解法中，常用的有塞德爾迭代法、共鈍斜量法、超松弛迭代法等等。此外，一些比較古老的普通消去法，如高斯法、追趕法等等，在利用計算機(jī)的條件下也可以得到廣泛的應(yīng)用。在數(shù)值計算方法中，數(shù)值逼近也是常用的基本方法。數(shù)值逼近也叫近似代替，就是用簡單的函數(shù)去代替比較復(fù)雜的函數(shù)，或者代替不能用解析表達(dá)式表示的函數(shù)。數(shù)值逼近的基本方法是插值法。初等數(shù)學(xué)里的三角函數(shù)表，對數(shù)表中的修正值，就是根據(jù)插值法制成的。在遇到求微分和積分的時候，如何利用簡單的函數(shù)去近似代替所給的函數(shù)，以便容易求微分和求積分,也是計算方法的個主要內(nèi)容。微分方程的數(shù)值解法也是近似解法。常微分方程的數(shù)值解法有歐拉法、預(yù)測校正法等。偏微分方程的初值問題或邊值問題，目前常用的是有限差分法、有限元素法等。有限差分法的基本思想是用離散的、只含有限個未知數(shù)的差分方程去代替連續(xù)變量的微分方程和定解條件。求出差分方程的解法作為求偏微分方程的近似解。有限元素法是近代才發(fā)展起來的，它是以變分原理和剖分差值作為基礎(chǔ)的方法。在解決橢圓形方程邊值問題上得到了廣泛的應(yīng)用。目前，有許多人正在研究用有限元素法來解雙曲形和拋物形的方程。數(shù)值計算方法的內(nèi)容十分豐富，它在科學(xué)技術(shù)中正發(fā)揮著越來越大的作用。數(shù)值計算方法的特點(diǎn)第一，面向計算機(jī)，要根據(jù)計算機(jī)特點(diǎn)提供實際可行的有效算法。即算法只能包括加、減、乘、除和邏輯運(yùn)算，是計算機(jī)能直接處理的。第二，有可靠的理論分析，能任意逼近并達(dá)到精度要求，對近似算法要保證收斂性和數(shù)值穩(wěn)定性，還要對誤差進(jìn)行分析。第三，要有好的計算復(fù)雜性，時間復(fù)雜性好是指節(jié)省時間，空間復(fù)雜性好是指節(jié)省存儲量，這也是建立算法要研究的問題，它關(guān)系到算法能否在計算機(jī)上實現(xiàn)。第四，要有數(shù)值實驗，即任何一個算法除了從理論上要滿足上述三點(diǎn)外，還要通過數(shù)值實驗證明是行之有效的。根據(jù)“數(shù)值計算方法''的特點(diǎn)，學(xué)習(xí)時我們要注意掌握方法的基本原理和思想，要注意方法處理的技巧及其與計算機(jī)結(jié)合，要重視誤差分析、收斂性及穩(wěn)定性的基本理論：其次，要通過例子，學(xué)習(xí)使用各種數(shù)值方法解決實際計算問題；最后，為了掌握本課的內(nèi)容，還應(yīng)做一定數(shù)量的理論分析與計算問題練習(xí)。由于本課內(nèi)容包括了微積分、代數(shù)、常微分方程的數(shù)值方法，以及高級程序設(shè)計的方法，讀者必須掌握這幾門課的基本內(nèi)容才能學(xué)好這一課程。誤差來源與誤差分析誤差的來源早在中學(xué)我們就接觸過誤差的概念，如在做熱力學(xué)實驗中，從溫度計上讀出的溫度是23.4C就不是一個精確的值，而是含有誤差的近似值。事實上，誤差的在我們的生活中無處不在，無處不有。如量體裁衣,量與裁的結(jié)果都不是精確無誤的，都有誤差。人們可能會問：如果使用計算機(jī)來解決問題，結(jié)果還會有誤差嗎？下面我們通過考察數(shù)學(xué)方法解決實際問題的主要過程來思考這個問題。用數(shù)學(xué)方法解決一個具體的實際問題，首先要建立數(shù)學(xué)模型，這就要對實際問題進(jìn)行抽象、簡化，因而數(shù)學(xué)模型本身總含有誤差，這種誤差叫做模型誤差。在數(shù)學(xué)模型中通常包含各種各樣的參變量，如溫度、長度、電壓等，這些參數(shù)往往都是通過觀測得到的，因此也帶來了誤差，這種誤差叫做觀測誤差。當(dāng)數(shù)學(xué)模型不能得到精確解時，通常要建立一套行之有效的數(shù)值方法求它的近似解，近似解與準(zhǔn)確解之間的誤差就稱為截斷誤差或方法誤差.由于在計算機(jī)中浮點(diǎn)數(shù)只能表示實數(shù)的近似值，因此用計算機(jī)進(jìn)行實際計算時每一步都可能有誤差，這種誤差稱為舍入誤差。

例如，函數(shù)#幻用泰勒(Taylor)展開式"必平"家一鏟近似代替，則數(shù)值方法的截斷誤差是(:介于o(:介于o與，之間)又如，在計算時用3.14159近似替代：產(chǎn)生的誤差*■<-3.14159=0.0000026 就是舍入誤差。上述種種意味著都會影響計算結(jié)果的準(zhǔn)確性，因此需要了解與研究誤差。在數(shù)值計算方法中將著重研究截斷誤差、舍入誤差，并對它們的傳播與積累作出分析。絕對誤差、相對誤差與有效數(shù)字.絕對誤差與絕對誤差限定義1.1設(shè)乂為準(zhǔn)確值，：為x的一個近似值，稱/?/一玄為近似值的絕對誤差，或誤差。通常我們無法知道準(zhǔn)確值2,也不能算出誤差的準(zhǔn)確值一，只能根據(jù)測量或計算估計出誤差的絕對值不超過某個正數(shù)鼠，則/稱為絕對誤差限。有了絕對誤差限就可知道x的范圍,--4*4即△落在區(qū)間［*?一『?內(nèi)。例如用毫米測度尺測量一長度？，讀出的長度為23mm,則有^一^“06mm。由此例也可以看到絕對誤差是有量綱和單位的。.相對誤差與相對誤差限只用絕對誤差還不能說明數(shù)的近似程度，例如甲打字時每百個字錯一個，乙打字時平均每千個字錯一個，他們的誤差都是錯一個，但顯然乙要準(zhǔn)確些。這就啟發(fā)我們除了要看絕對誤差大小外，還必須顧及量的本身。定義1.2把近似值的誤差￡與準(zhǔn)確值工的比值稱為近似值T的相對誤差，記作實際計算時，由于真值工通常是不知道的，通常取'7。相對誤差也可正可負(fù)，它的絕對值的上■?■8 1* '|t*| 4—,1%界叫做相對誤差限。記作即II。根據(jù)定義，甲打字時的相對誤差 100 ,乙打字時的r；5-L-01%相對誤差 wooo易知相對誤差是個無量綱量。.有效數(shù)字當(dāng)下有多位時，常常接四舍五入的原則得到*的前幾位近似值例如x-x-314159265*-取3位k；?3.?440002；取5位它們的誤差都不超過末位數(shù)字的半個單位，即卜一生H|￡Lx10、,卜-3J4iq4lxVT*2 2現(xiàn)在我們將四舍五入抽象成數(shù)學(xué)語言，并引入一個新名詞“有效數(shù)字”來描述它。定義1.3若近似數(shù)T的誤差限是某一位的半個單位，該位到丁的第一位非零數(shù)字共有q位，我們就說，?有七位有效數(shù)字。如取x?“殳14作丁的近似值，工?就有3位有效數(shù)字；取x'=3.M]6作：的近似值，丁就有5位有效數(shù)字。有效數(shù)字也可采用以下定義:X的近似數(shù)X*寫成標(biāo)準(zhǔn)形式:其中,是0到9的一個數(shù)字，且為不為0,全為整數(shù)，若(1.2)則稱，有。位有效數(shù)字。例1.1依四舍五入原則寫出下列各數(shù)具有5位有效數(shù)字的近似數(shù)，913.95872,39.1882,0.0143254,8.000033解：按定義，上述的5位有效數(shù)字的近似數(shù)分別為913,96,39.188,0.014325,8.0000注意，8.000033的5位有效數(shù)字的近似值是&0000而不是8,8只有一位有效數(shù)字。.有效數(shù)字與誤差關(guān)系不難看出，若由（1.1）給出某近似數(shù)有V位有效數(shù)字，則可以從（1.2）求得這個數(shù)的絕對誤差限2因此，在色相同的情況下，':越大則IL就越小，故有效數(shù)字越多，絕對誤差限越小。定理1.1用（1.1）表示的近似數(shù)/，若工”具有,：位有效數(shù)字，則其相對誤差限為陰. —■—*IOT*** .反之，若f的相對誤差限 則T■至少具有工位有效數(shù)字。證明由式（1.1）得所當(dāng)？‘有暗位有效數(shù)字時所當(dāng)？‘有暗位有效數(shù)字時IH卜】戶12?冒,球產(chǎn)?如。3<1--—心2，反之，若一的相對誤差限 2gl+D有：卜"卜kka><(Wl卜"卜kka><(Wl>x1oWX2^由式（1.2）式知道，寸至少有三位有效數(shù)字，證畢這說明有效數(shù)字越多，相對誤差限越小。例1.2要使亞的近似值的相對誤差限小于om,要取兒個有效數(shù)字？解：因為屈.4.47…，所以q7,令—XI0*^01%

2"故取K-4即可滿足。5.函數(shù)計算的誤差估計一般情況，當(dāng)自變量有誤差時計算相應(yīng)的函數(shù)值也會產(chǎn)生誤差，其誤差限可利用函數(shù)的泰勒展開式估計。設(shè)“X）是一元函數(shù)，工的近似值為味以/（,）近似,㈤，其誤差界記作式/6?可用泰勒展開式/W 力-八戲”-力+嬰Qr尸:介于了，丁之間，取絕對值得假定義電與/*（吟的比值不太大，可忽略貳”■）的高階項，于是可得計算函數(shù)值的誤差限加**力W4解|aT%#a.>例1.3設(shè)/的相對誤差限為2%,則/"）.（*?『的相對誤差限是多少？解：陽m”戶浮Cx2M?“g?所以/（0的相對誤差限為0.02no當(dāng)J為多元函數(shù)時，如計算?■/（%卬…，4）。如果、4…'*"的近似值為年*…H,則看的近似值尸心，于是函數(shù)值三*的誤差心"）山泰勒展開式得/瓜,冷….4)-/(。卬….Q于是誤差限為如疆*力例1.4已測得某場地長’的值為寬丁的值為已知卜?。?。**?,試求面積的絕對誤差限馬相對誤差限。s-4/.— -/解：因也㈤，那么3牌卜）俘I"）其中令，■，-80w,琮），?r-110*1于是絕對誤差限?（『）-?*02*110x01?2X6相對誤差限1.3誤差傳播與若干防治辦法由前述可知，在數(shù)值計算中每步都可能產(chǎn)生誤差。而一個問題的解決，往往要經(jīng)過成千上萬次運(yùn)算，我們不可能每步都有加以分析，下面，通過對誤差的某些傳播規(guī)律的簡單分析，指出在數(shù)值計算中應(yīng)注意的幾個原則，它有助于鑒別計算結(jié)果的可靠性并防止誤差危害現(xiàn)象的產(chǎn)生。.要避免兩相近數(shù)相減在數(shù)值運(yùn)算中兩相近數(shù)相減會使有效數(shù)字嚴(yán)而損失，例如*'5父85.尸?53252都具有五位有效數(shù)字，但只有兩位有效數(shù)字，所以要盡量避免這類運(yùn)算。1g%-1g嗎-1g—通常采用的方法是改變計算公式，例如當(dāng)天與勺很接近時，由于 巧那么可用右端的公式代替左端的公式計算，有效數(shù)字就不會損失。當(dāng)工很大時， 必也可用右端來代替左端。一般情況，當(dāng)時，可用泰勒展開，）+，）+q：）（x-/+…取右端的有限項近似左端。如果計算公式不能改變，則可采用增加有效位數(shù)的方法。.要防止大數(shù)“吃掉”小數(shù)若參加運(yùn)算的數(shù)的數(shù)量級相差很大，而計算機(jī)的位數(shù)有限，如不注意運(yùn)算次序，就可能出現(xiàn)大數(shù)“吃掉“小數(shù)的現(xiàn)象，影響計算結(jié)果。例如在五位十進(jìn)制計算機(jī)上，計算寫成規(guī)格化形式MM $A-0.52492?10,+^0.000001x10由于計算時要對階，在計算機(jī)中表示為o,計算出來A?0524S=1Q‘，結(jié)果嚴(yán)重失真!m小如果計算時，先將H計算出來，再與52492相加，就不會出現(xiàn)大數(shù)“吃掉”小數(shù)的現(xiàn)象了。.注意簡化計算步驟，減少運(yùn)算次數(shù)同樣一個計算題，如果能減少運(yùn)算次數(shù)，不但可節(jié)省計算機(jī)的計算時間，還能減少舍入誤差。例如，計算x255的值，如果逐個相乘要用254次乘法，但若寫成產(chǎn)-X//只要做14次乘法運(yùn)算即可。.絕對值太小的數(shù)不宜作除數(shù)?X ?X?十X*—設(shè)T與；分別有近似值 ，的近似值，，其絕對誤差F?卜卜廿聲，M

i/r1?嚴(yán)?顯然,當(dāng)g很小時，近似值工”的絕對誤差,仁）有可能很大。因此，不宜把絕對值太小的數(shù)作除數(shù)。認(rèn).，■爺? ■631%KrI。W第2章非線性方程求根與線性方程相比，非線性方程問題無論是從理論上還是從計算公式上，都要復(fù)雜得多。對于一般的非線性方程/（力?0,計算方程的根既無一定章程可尋也無直接法可言。例如，求解高次方程組7/?x-L5=0的根，求解含有指數(shù)和正弦函數(shù)的超越方程/■幻的零點(diǎn)。解非線性方程或非線性方程組也是計算方法中的?個主題。一般地，我們用符號/（*）來表示方程左端的函數(shù)，方程--般形式表示為〃工）?°,方程的解稱為方程的根或函數(shù)的零點(diǎn)。通常，非線性方程的根不止一個，對于非線性方程，一般用迭代法求解。因此，在求解非線性方程時,要給定初始值或求解范圍。實根的對分法使用對分法的條件對分法或稱二分法是求方程近似解的一種簡單直觀的方法。設(shè)函數(shù)人力在〔。力1上連續(xù)，且/㈤則在［6網(wǎng)上至少有-零點(diǎn)，這是微積分中的介值定理，也是使用對分法的前題條件。計算中通過對分區(qū)間、縮小區(qū)間范圍的步驟搜索零點(diǎn)的位置。例2.1用對分法求解,（x）*37.7/+19.2x-153在區(qū)間［詞之間的根.解：/G）--2-8/o-Q3由介值定理可得有根區(qū)間1層四?［閉。（2）計算均-。+期2-1.“。9?《45有根區(qū)間［?］叩.切（3）計算三=（16+2）/2=1.75,/Q4=0.078125,有根區(qū)間—。工工了可。（4）?直做到y(tǒng)（/）K.（計算前給定的精度）或I。一坪￡時停止。2.1.2對分法求根算法計算/（動■。的一般計算步驟如下：.輸入求根區(qū)間SQ1和誤差控制量凡定義函數(shù)〃力。IF THEN退出選用其他求求根方法.WH/LE?時（1）計算中點(diǎn)*.g*劫,2以及/（目的值；（2）分情況處理l/WK?:停止計算，?x,轉(zhuǎn)向步驟4/（?VW（Q修正區(qū)間141卜>（??*!/W/W<0修正區(qū)間［3卜》［@*1ENDWHILE.輸出近似根在算法中，常用血貳/3）噸底/@>°代替八視/⑸>°的判斷，以避免/1（?）?/（力數(shù)值溢出。對分法的算法簡單，然而，若/（*）在S/1是有幾個零點(diǎn)時，只能算出其中一個零點(diǎn)；另一方面，即使〃力在上有零點(diǎn)，也未必有/觸怨）〈匕這就限制了對分法的使用范圍。對分法只能計算方程的實根。2.2迭代法對給定的方程將它轉(zhuǎn)達(dá)換成等價形式：給定初值外，由此來構(gòu)造迭代序列尸，如果迭代法收斂，即出3,此武與）”，尸，如果迭代法收斂，即出3,此武與）”，有a?a°）,則現(xiàn)就是方程(2)--5(2)--52r?P-5,迭代格式 2人才■。的根。在計算中當(dāng)I&*~匕1小于給定的精度控制量時，取@為方程的根。例如，代數(shù)方程？-2x-5?0的三種等價形式及其迭代格式如下:(1)/■&f+5.x■胱〃，5.迭代格式%.4匕*+5(1)(3)迭代格式(3)對于方程構(gòu)造的多種迭代格式怎樣判斷構(gòu)造的迭代格式是否收斂？收斂是否與迭代的初值有關(guān)？定理2.1 若定義在［。為|上，如果滿足（D當(dāng)有2］有《配（2）KG在上可導(dǎo)，并且存在正數(shù)"1,使對任意的問6村，有IrlWl"；則在［。.句上有惟一的點(diǎn)片滿足,"?<**>,稱寸為K?的不動點(diǎn)。而且迭代格式“_對任意的初值均收斂于?<*）的不動點(diǎn)黑■,并有誤差估計式證明：（1）先證明存在性： ）?x-a*）,則有/何…故有a 使得（2）再證明惟一性：設(shè)都是?<*）的不動點(diǎn)，且父**?，則有與假設(shè)矛盾，這表明*~xn,即不動點(diǎn)是惟一的。(3)當(dāng)*0—同時，由于?x)w[a向可用歸納法證明，迭代序列QJu(?網(wǎng)，于是由微分中值定理加和得"|--，卜4川|勺-，| (2.2)因為4<1,所以當(dāng)北~>9時，""70%"*’即迭代格式加?^^收斂。(4)誤差估計式：設(shè)A固定，對于任意的正整數(shù)3有，1%,-qH7V-9.1Hl+1+-￡(*4+tf*Z 卜注：定理2.1是判斷迭代法收斂的充分條件，而非必要條件。要構(gòu)造滿足定理條件的等價形式一般是難于做到的。事實上，如果寸為JQ)的零點(diǎn)，若能構(gòu)造等價形式“女工而"2KI,由。口)的連續(xù)性，一定存在￡的鄰域仁?艮八「)，其上有I"（才KL〈i,這時若初值迭代也就收斂了。由此構(gòu)造收斂迭代式有兩個要素，其一，等價形式*■?*）應(yīng)滿足i.s'ki；其二，初值必須取自寸的充分小鄰域，這個鄰域大小決定于函數(shù)/a）,及做出的等價形式“■??。例2.2求代數(shù)方程？-2x-5-0,在4?2附近的實根。解：1）2?2/5―廣取%+5、2 1 T’（2x+/.（杯】，當(dāng)x€[l.*25]，構(gòu)造的迭代序列收斂。取勺.？則=2.08008, 4=2.09235, =2.094217"*=2.094494,馬=2.094543,三=2.094550準(zhǔn)確的解是工=2.094551481502）將迭代格式寫為亡-5,%a?-5

54--^-0任）?一^-力。；。）卜容卜產(chǎn)口」2S]’?迭代格式.仍（Q不能保證收斂。對于收斂的迭代過程，只要迭代足夠多次，就可以使結(jié)果達(dá)到任意的精度，但有時迭代過程收斂緩慢,從而使計算量變得很大，因此迭代過程的加速是個重要的課題。以下介紹一種埃特金(Aitken)方法。對X.山)方程，構(gòu)造加速過程，算法如下：(1)預(yù)測：(2)校正：(3)改進(jìn)：有些不收斂的迭代法，經(jīng)過埃特金方法處理后，變得收斂了。例2.3求方程J? 在勺附近的根一。解：若采用迭代公式："=匕7迭代法是發(fā)散的，我們現(xiàn)在以這種迭代公式為基礎(chǔ)形成埃特金算法:取&■15計算結(jié)果如表2.1所示:表2.1計算結(jié)果t%舉?】01.512.375001339651.4162921.840925.238881.3556531.491402.317281.3289541.347101.444351.3248051.325181.327141.32472我們看到，將發(fā)散的迭代公式通過埃特金方法處理后，競獲得了相當(dāng)好的收斂性。2.4牛頓迭代法對方程」(*)‘°可構(gòu)造多種迭代格式，牛頓迭代法是借助于對函數(shù),Q)■°作泰勒展開而構(gòu)造的一種迭代格式。將八外?。在初始值工作泰勒展開：/w- a*彎.--專尸*??4■取展開式的線性部分作為的近似值，則有設(shè)fa/o,則勺/W令飛■/6)類似地，再將/(*)-0在均作泰勒展開并取其線性部分得到：一宜做下去得到牛頓迭代格式：■?4? 上?L2…/⑷ , (2.3)牛頓迭代格式對應(yīng)于/?)?°的等價方程是

(2.4)cr?)a(2.4)若M是/㈤的單根時，/8）?0./9）"°則仃*4°,只要初值三充分接近工，所以牛頓迭代收斂。當(dāng)二為的五重根時，取下面迭代格式：比72“牛頓法的幾何意義以?。ㄉ祝樾甭首鬟^點(diǎn)的切線，即作了（外在士點(diǎn)的切線方程令*-Q,則得此切線與工軸的交點(diǎn)電，即再作在《》/（砧）處的切線，得交點(diǎn)與，逐步逼近方程的二根.如圖2.2所示圖2.2牛頓切線法示意圖在區(qū)域［勺，4 的局部“以直代曲'’是處理非線性問題的常用手法。在泰勒展開中，截取函數(shù)展開的線性部分替代/㈤。例2.4用牛頓迭代法求方程/（x）-/-7.7/+19一2X-I5.3,在、附近的根。-7.7^*19.2^-15.3解：XA―乂754—192計算結(jié)果列于表2.2中。表2.2計算結(jié)果Kxkf(x)01.00-2.811.41176-0.72707121.62424-0.14549331.6923-0.013168241.69991-0.000151551.70比較表2.1和表2.2的數(shù)值，可以看到牛頓迭代法的收斂速度明顯快于對分法。牛頓迭代法也有局限性。在牛頓迭代法中，選取適當(dāng)?shù)跏贾等乔蠼獾那邦}，當(dāng)?shù)某跏贾等稍谀掣母浇鼤r迭代才能收斂到這個根，有時會發(fā)生從一個根附近跳向另一個根附近的情況，尤其在導(dǎo)數(shù)數(shù)值很小時，如圖2.3所示。圖2.3失效的牛頓迭代法牛頓迭代算法.定義函數(shù)/(外，?(?),輸入或定義迭代初始值和控制精度epsilon。.FORSOtoMAXREPTx_k\-g(x_JW) "g(才是牛慎迭代公式IF(\x_*1-x_JtO|<apsdoM)oRa/(x_JtQK^oK)THEN｛輸出滿足給定精度的近似解結(jié)束｝ENDFOR.輸出：在J*0附近無根。注：在偽碼中，"/r'表示注釋語句。偽碼又稱過程設(shè)計語言，主要用于描述算法。它是某種高級語言和自然語言的混雜語言，它取某種高級語言中的一些關(guān)鍵字，用于描述算法的結(jié)構(gòu)化構(gòu)造和數(shù)據(jù)說明等。偽碼的語句中嵌有自然語言的敘述，偽碼易于2.5弦截法弦截法迭代格式j(luò)-X.2--TOC\o"1-5"\h\z在牛頓迭代格式中： 『g\o"CurrentDocument",一㈤?"32 〃、用差商 q-Ri代替導(dǎo)數(shù)并給定兩個初始值七和工】，那么迭代格式可寫成如下形式：〃必-〃2 （2.5）上式稱為弦截法。用弦截法迭代求根，每次只需計算一次函數(shù)值，而用牛頓迭代法每次要計算一次函數(shù)值和一次導(dǎo)數(shù)值。但弦截法收斂速度稍慢于牛頓迭代法。弦截法的幾何意義做過兩點(diǎn)（、?/6》和的一條直線（弦），該直線與％軸的交點(diǎn)就是生成的迭代點(diǎn)為，再做過（%/（??和（馬的一條直線，力是該直線與工軸的交點(diǎn)，繼續(xù)做下去得到方程的根/依）-0,如圖2.4所示圖2.4弦截法示意圖例2.5用弦截法求方程〃X）■/-7.7?*19.2x-15.3的根，取勺.1.5.-4.0。解：“‘

計算結(jié)果列于表2.3中。表2.3計算結(jié)果上二/W01.5-0.45142.321.909090.24883531.65543-0.080569241.717480.028745651.701160.0019590261.69997-0.000053924671.79.45910-8弦截法算法.定義函數(shù)輸入或定義控制精度值epsilon和迭代初始值計算/匕?/。_初 表示％.勺.FOR:=0,1…MAXREPT/21?/?_北為2->x_A:-x_Jt2-/2(x_A2〃/_上表示4MTHEN｛輸出滿足給定精度的近似解結(jié)束｝2.4x_kl-x_k2x_k2:~x_k〃為下一次迭代準(zhǔn)備數(shù)值ENDFOR.輸出：在初始值、?劫*?上2附近無根。2.6非線性方程組的牛頓方法為了敘述簡單，我們以解二階非線性方程組為例演示解題的方法和步步驟，類似地可以得到解更高階非線性方程組的方法和步驟。設(shè)二階方程組以5?。 (2.6)其中X?，為自變量。為了方便起見，將方程組寫成向量形式：“您做fl舊仁川，其中二?J將附近作二元泰勒展開，并取其線性部分，得到方程組4(勺㈤+(r)當(dāng)型?。-%)絲-oTOC\o"1-5"\h\zox dy〃勺㈤WR好必”卬逐衿?0(M .令X-(.&.y一“則有dx a&當(dāng)也唱出rg川w qy

弧也生生如果dx分，解出&A,得再列出方程組：曄小巡守32?_〃”)n dy‘^2q --工的比M 夕解出,&7r 修?k田得出繼續(xù)做下去，每一次迭代都是解一個方程組…腳￡圖力(2.7)&卬Ar.九*-%即“?".&①”?為.2,直到為止。例2.6求解非線性方程組4GM?4---y-0/式八力-y-0取初始值

畫漁，〃、班打 -2y\解:玄笠-lJ5Sy解:"43(；：sa：：)(-2&+32--OLll-27l328Ar-^>-001828產(chǎn).00042%]a(qJI-0.029849J解方程得'■國曾…斗尸力5J l-,7J(-0-029849J(-1.729849J繼續(xù)做下去，直到mW&U》D<l°"時停止。將兩個變量的非線性方程組推廣到多變量的非線性方程組：人風(fēng),盯?…。-0((X*,-0 (28)記*…的向量形式:F(X)-0

的向量形式:F(X)-0用牛頓迭代法求方程組的解為:就⑺彌X) %⑺% I %跳⑶里(箱...秘⑶

血 與 改乜⑺孤就⑺彌X) %⑺% I %跳⑶里(箱...秘⑶

血 與 改乜⑺孤X)..^(X)"X）?其中：稱J（x）為雅可比（jacobi）矩陣。計算中采用4獷乂/川-義與?化專NjMZ'-一或評 （29）一直做到小于給定精度為止。在上的領(lǐng)域中若”《也"AT,則迭代收斂。.7程序示例程序2.1用牛頓迭代法求/S）在工:附近的根。算法描述給定/（*）,從三開始，根據(jù)牛頓迭代公式計算」(x)在L附近的根。程序源碼IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII// Purpose:牛頓迭代法求根 〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h#include<math.h>//iterate。)為"#defineiterate(x)(x-((x*x-3)*x-1)/(3*x*x-3))#define/(力 ((x*x-3)*x-l)#definexO1.5 〃迭代初值xO#defineMAXREPT1000//最大迭代次數(shù)#defineepsilon0.0001 //精度voidmain(){inti;doublex_k=x0,x_kl=x0for(i=i<MAXREPT；i++)(printfCtGot...%f7n,\x_kl);x_k1=iterate(x_k); 〃迭代iffabs(x_kl-x_k)<epsilon||fabs(fi[x_kl)<epsilonprintfC'!Root:%f\n”,x_kl);〃滿足精度，輸出return;{x_k=x_kI{printfi(tfcAfter%drepeate,nosolvedAm'\MAXREPT)；// EndofFile 計算實例已知取4 用牛頓迭代法計算/㈤的根。程序輸入輸出由本程序預(yù)設(shè)的及小.L5,得到輸出!Root:1.879385程序2.2用弦截法求/(“)在內(nèi)通附近的根。算法描述給定/Q),從、?再開始，根據(jù)弦截法迭代公式求得/(*)在其附近的根。程序源碼〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃〃///Purpose：弦截法求根〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h>#include<math.h>#definef(x)(x*x*x-7.7*x*x+19.2*x-15.3)//f函數(shù)#definexO0.0〃初值xO,x1#definexl1.0#defineMAXREPT1000 //最大迭代次數(shù)#defineepsilon0.00001〃求解精度voidmain(){inti;doublex_k=x0,x_k1=x1,x_k2=x1;fbr(I=O;i<MAXREPT;i++={printf("!Root:%f\n”,x_k2);x_k2=x_kl-(fi[x_kl)*x_k1?x_k))/(f(x_k1)?Rx_k));〃弦截求新x_niffabs(x_k2-x_k1)<epsilon||fabs(f(s_k2))<epsilon{printf("!Root:%f\n,,,x_k2);〃滿足精度,輸出return;(x_k=x_kl;x_kl=x_k2; 〃反復(fù){printfi(t4After%drepdate,nosolvde.'n",MAXREPT);// EndofFile 計算實例用弦截法求方程〃x)?J-7.7/+19.2x-15.3的根，取.-I'q程序輸入輸出由本程序的,Q)及4?'得到輸出!Root:1.700000第3章解線性方程組的直接法在近代數(shù)學(xué)數(shù)值計算和工程應(yīng)用中，求解線性方程組是雨:要的課題。例如，樣條插值中形成的“關(guān)系式，曲線擬合形成的法方程等，都落實到解一個v元線性方程組，尤其是大型方程組的求解，即求線性方程組(3.1)的未知量的數(shù)值?！鲈?aw+—>a.x.-4%4*31*-***-A與*""，-4 (31)其中Gjbi為常數(shù)。上式可寫成矩陣形式4x=b,即其中，AF哂為系數(shù)矩陣("吃"，產(chǎn)為解向量…為常數(shù)向量。當(dāng)deS=￡)0時，由線性代數(shù)中的克萊姆法則，方程組及的解存在且惟?,且有q為系數(shù)矩陣三的第7列元素以◎代替的矩陣的行列式的值。克萊姆法則在建立線性方程組解的理論基礎(chǔ)中功不可沒，但是在實際計算中，我們難以承受它的計算量。例如，解一個100階的線性方程組，乘除法次數(shù)約為(101T00!-99),即使以每秒】?的運(yùn)算速度，也需要近產(chǎn)年的時間。在石油勘探、天氣預(yù)報等問題中常常出現(xiàn)成百上千階的方程組，也就產(chǎn)生了各種形式方程組數(shù)值解法的需求。研究大型方程組的解是H前計算數(shù)學(xué)中的一個重要方向和課題。解方程組的方法可歸納為直接解法和迭代解法。從理論上來說，直接法經(jīng)過有限次四則運(yùn)算，假定每一步運(yùn)算過程中沒有舍入誤差，那么，最后得到方程組的解就是精確解。但是，這只是理想化的假定，在計算過程中，完全杜絕舍入誤差是不可能的，只能控制和約束由有限位算術(shù)運(yùn)算帶來的舍入誤差的增長和危害，這樣直接法得到的解也不一定是絕對精確的。迭代法是將方程組的解看作某種極限過程的向量極限的值，像第2章中非線性方程求解一樣，計算極限過程是用迭代過程完成的，只不過將迭代式中單變量X**】換成向量而已。在用迭代算法時，我們不可能將極限過程算到底，只能將迭代進(jìn)行有限多次，得到滿足一定精度要求的方程組的近似解。在數(shù)值計算歷史上，直接解法和迭代解法交替生輝。一種解法的興旺與計算機(jī)的硬件環(huán)境和問題規(guī)模是密切相關(guān)的。一般說來，對同等規(guī)模的線性方程組，直接法對計算機(jī)的要求高于迭代法。對于中等規(guī)模的線性方程組伊^^，由于直接法的準(zhǔn)確性和可靠性高，一般都用直接法求解。對于高階方程組和稀疏方程組（非零元素較少），一般用迭代法求解。消元法三角形方程組的解形如下面三種形式的線性方程組較容易求解。對角形方程組TOC\o"1-5"\h\z%^勺 ■與■■' %勺-X （3.3）設(shè)對每一個方程，/,取,TN…』。顯然，求解〃階對角方程的運(yùn)算量為下三角方程組，馬 -4* ■巧 "與+ 4勺*…? ? 4 （3.4）按照方程組的順序，從第一個方程至第K個方程，逐個解出由方程得看將』的值代入到第二個方程?'gK,■幻得 勺-眄-’11砧〃必將的值代入到第？個方程。凝+…一x「孰■IN”-m得 P計算不需要,次乘法或除法運(yùn)算，i-t2-,■**.因此，求解過程中的運(yùn)算量為g?罕?*)上三角方程組%-*x* … *F.4* - *%。_%<…433* -如=% (3.5)與計算下三角方程組的次序相反，從第七個方程至第一個方程，逐個解出2j.從"-1?…2?L由第t個方程%。?.待。將三的值代入到第4-1個方程將、?/士…Eu的值代入到第;個方程F0-4得解的通式—部W.，…V.A1計算號需要“T+1次乘法或除法運(yùn)算1**?"-L…?2?L因此求解過程中的運(yùn)算量為7?I?￡/?心；D?吁）1 /-I 2消元法的基本思想就是通過對方程組做初等變換，把?般形式的方程組化為等價的具有上述形式的易解方程組。3.1.2高斯消元法與列主元消元法高斯消元法高斯消元法是我們熟悉的古老、簡單而有效的解方程組的方法。下面是中學(xué)階段解二元方程組（高斯消元法）的步驟：fx,-X,-1 （3.6）13X（*2x）-8 （3.7）方程（3.6）乘以-3加到第（3.7）個方程中得03-57代入（3.6）得“IN其方法相當(dāng)于對方程組的增廣矩陣做行的初等變換:

（3。：：）3』（；；；卜史）W已是上三角矩陣，而為原方程組的等價方程組，已化成易解的方程組形式。再用回代方法求解，得到：-2這就是高斯消元法解方程組的消元和回代過程。一般地，可對線性方程組（3.1）施行以下一系列變換:（1）對換某兩個方程的次序;（2）對其中某個方程的兩邊同乘?個不為零的數(shù);（3）把某?個方程兩邊同乘?個常數(shù)后加到另一個方程的兩邊。記變換后的方程組為:才“%*-?!。-E(38)顯然方程組（3.1）與（3.8）是等價方程組,顯然方程組（3.1）與（3.8）是等價方程組,的增廣矩陣為:可以看出，J實際上是由《4與按-系列初等換后得到的(1)對換(/楊某兩行元素；(2)(九力中的某行乘一個不為零的數(shù)；(3)把(4年的某一行乘一個常數(shù)后加到另一行。高斯消元法就是通過以上(3)的變換，把《49化為等價的上三角形式。下面我們以"?彳為例演示消元過程。設(shè)方程組：ZlA -瓦0M巧*0^X2+?占+0Mx.-6,?Mi*<V3+.5F玉f%。+4用+3,-4 (39)其增廣矩陣為：(1)若,則將第一行乘以-%/%i加到第二行上；將第乘以-加到第三行上；將第一行乘以一加到第四行上得到a*-呼其中:a,..-egj.2%*f°叱-5%/,/-2,X4⑵若詞“a則將第二行乘以一造,右加到第三行上；將第一行乘以-馬‘啰加到第四行上,得到all

00(3.11)0(3.11)其中:譚■碎-啰-3.4

中-if9-岬噌3,4⑶鏟ff"則將第三行乘以一點(diǎn)/蟾*加到第四行上，得到(3.12)其中喑?虐-噌喈礎(chǔ)＞--已是上三角矩陣，于是得到了與原方程等價的易解形式的方程組：再對方程組（3.13）依次回代解出

從式(3.12)可以得到系數(shù)矩陣＜的行列式的值為彳的對角元素的乘積。即這也正是在計算機(jī)上計算方陣W的行列式的常規(guī)方法。要將上述解方程步驟推廣到三階方程組，只需將對控制量“4”的操作改成對控制量V的操作即可。t元方程組高斯消元法的步驟如下：對于士-L2”.〃-l約定《■。有6?鏟，鏟/邛)中6?鏟，鏟/邛)中**l<i<m經(jīng)過以上消元，我們已得到與出等價的方程組后.A,其中之已是一個上三角矩陣。為簡單起見，仍記2的元素為“宕的兀素力當(dāng)月歸即己得到原方程組的解。高斯消元法算法1.輸入：方程組階數(shù)〃，方程組系數(shù)矩陣4和常數(shù)向量入2.FORk:=lTOn-1〃消元過程2.FORk:=lTOn-1〃消元過程FORi:=k+lTOn[1■/aFORj:=k+lTOn

{■產(chǎn)}〃endforj//ENDFORI//ENDFORK3.FORi:=nTOI〃回代求解3.FORi:=nTOI〃回代求解{s:=biFORj:=i+lTOnDO4.輸出方程組的解勺J-12…。高斯消元法的運(yùn)算量整個消元過程即式（3.14）的乘法和除法的運(yùn)算量為JU4 Al 4回代過程即式（3.15）的乘除運(yùn)算量為網(wǎng)”.1）-2~故高斯消元法的運(yùn)算量為(3.16)3+- ■g,）

(3.16)高斯消元法的可行性在上面的消元法中，未知量是按照在方程組中的自然順序消去的，也叫順序消元法。在消元過程中假定對然元素小 消元步驟才能順利進(jìn)行，由于順序消元不改變￡的主子式值,故高斯消元法可行的充分必要條件為W的各階主子式不為零。但是，實際上只要修才~“方程組就有解。故高斯消元法本身具有局限性。另一方面，即使高斯消元法可行，如果簿很小，在運(yùn)算中用它作為除法的分母會導(dǎo)致其它元素數(shù)量級的嚴(yán)重增長和舍入誤差的擴(kuò)散。這是高斯消元法的另一缺陷。(0.0003X.*10000x.-20001 (3J7)L000011*10OOOx,-1.0000 (3.18)的精確解為:2-3在高斯消元法計算中取5的精確解為:2-3在高斯消元法計算中取5位有效數(shù)字。解：方程(3.17)x(-1)/0.0003+方程(3.18)得:09003不?39000*,-200019999.0?s-6666.0,-.Xj-0.M67代入方程(3.17)得“】?°。由此得到的解完全失真，如果交換兩個方程的順序,得到等價方程組I.OOOOAj*1.0000^-1.00W

10.0003jh+39000啊-20001經(jīng)高斯消元后有(10000、+10000x2-100002勞的》?1詡8.-.Xj-06667,?-0.3333由此可看到，在有些情況下，調(diào)換方程組的次序?qū)Ψ匠探M的解是有影響的，對消元法中抑制舍入誤差的增長是有效的。如果不調(diào)換方程組的次序，取6位有效數(shù)字計算方程組的解，得到取9位有效數(shù)字計算方程組的解，得到力?06666￡7.X1-0333333由此可見有效數(shù)字在數(shù)值計算中的作用。列主元消元法如果在一列中選取按模最大的元素，將其調(diào)到主干方程位置再做消元，則稱為列主元消元法。調(diào)換方程組的次序是為了使運(yùn)算中做分母量的絕對值盡量地大，減少舍入誤差的影響。用列主元方法可以克服高斯消元法的額外限制，只要方程組有解，列主元消元法就能暢通無阻地順利求解，同時又提高了解的精確度。max{kD?kl更具體地，第一步在第一列元素中選出絕對值最大的元素,交換第一行和第噌行的所有元素，再做化簡嗎為零的操作。對于每個七狂12…-T在做消元前，選出閭…中絕對值最大的元素對既行和航行交換后，再做消元操作，這就是列主元消元法的操作步驟。由于修//0,可證<?射，蜀震…中至少有一個元素不為零，從理論上保證了列主元消元法的可行性。列主元消元法與高斯消元法相比，只增加了選列主元和交換兩個方程組（即兩行元素）的過程。列主元消元法算法.輸入：方程組階數(shù)用，方程組系數(shù)矩陣W和常數(shù)向量項￡。.FORt=lTOd-1 〃選主元的消元過程{〃選擇. 小K*lFOR?=k+!TOnIFk|〉6THBH{胸?《；6?卜JFORk^TOq //交換第七行和第z行9f；@a：?7IFORr=k*lTOn（匕?。/小FORj:=kHTOn{,?％-?*/}ffENDPORJW4} //ENDFORI} //ENDFORK3.FORi:=nTO1 //回代求解4.輸出方程組的解如果對于第克步，從全行至七行和從裊列至y列中選取按模最大的wW5記月卜蚓1aH對第—行交換,對第i列交換,這就是全主元消元法。在《列和第V列交換時，還要記錄下v的序號，以便恢復(fù)未知量M和XV的位置。3.1.3高斯一若爾當(dāng)(Gauss—Jordan)消元法高斯消元法將系數(shù)矩陣化為上三角矩陣，再進(jìn)行回代求解；高斯一若爾當(dāng)消元法是將系數(shù)矩陣化為對角矩陣，再進(jìn)行求解，即對高斯消元法或列主元消元法得到的等價增廣矩陣：alt%flU *10圖4第b?0。電審Ap)00 0a曾即用第4行乘右加到第3行上,用第4行乘如曾加到第2行上,用第4行乘布里加到第1行上，得到?0°裁0aSe0守0。嫉。守用第3行乘號'卡加到第2行上，用第3行乘-媾加到第1行上，再用第2行乘一//姆加到第1行上，得到'.0 0 06fo0喈0Q中0。點(diǎn)。理0 0 。眼,用 (3,19)為方便起見，仍記式(3.19)系數(shù)矩陣元素為氣,常數(shù)項元素為當(dāng)。那么用初等變換化系數(shù)矩陣為對角矩陣的方法稱為高斯一若爾當(dāng)消元法。從形式上看對角矩陣（3.15）比上三角矩陣（3.11）更為簡單，易于計算三但是將上三角矩陣（3.11）化為對角矩陣（3.15）的工作量略大于上三角矩陣回代的工作量。高斯一若爾當(dāng)消元法求逆矩陣設(shè)a為非奇異矩陣，方程組/**/?的增廣矩陣為,-（工〃）。如果對。應(yīng)用高斯-若爾當(dāng)消元法化為則小，123,/-245例3.2用高斯-若爾當(dāng)消元法求P56J的逆矩陣TOC\o"1-5"\h\z1 0 -1/2fo -5/2 2->0 1 3/20 3/2 -I0 0 1/21 -172 03.2直接三角分解法

相當(dāng)于用了三個初等矩陣左乘山和土。記仍以為例，在高斯消元法中，從對方程組進(jìn)行初等變換的角度觀察方程組系數(shù)矩陣的演變過程。相當(dāng)于用了三個初等矩陣左乘山和土。記第一次消元步驟將方程組（3.9）化為方程組（3.10）,〃.. .2*4,容易驗證,1000W10 1000W10 100 0由卒”■Q?得到■6嗎其中將方程組（3.10）化為方程組（3.11）,相當(dāng)于用了兩個初等矩陣左乘士和E。記Q?咽，渡八切.有I0 I0 00,0 1000 0 100-3014 0 0 0、0 10 00 -Zn I 00 0 0 10 0同理，將方程組(3.11)化為方程組(0 0同理，將方程組(3.11)化為方程組(3.12),相當(dāng)于用一個初等矩陣左乘完成了消元過程，有亦有所有消元步驟表示為：（4與左乘一系列下三角初等矩陣。容易驗證，這些下三角矩陣的乘積工仍為下三角矩陣，并有「I000'*11 *12 ■ °/?IGiLq彳彳/■a于是有X-7X或這里『‘仍為下三角矩陣，其對角元素為1,稱為單位下三角矩陣，而三已是上三角矩陣。記L?￡,]■",則有A-LU結(jié)果表明若消元過程可行，可以將三分解為單位下三角矩陣上與上三角矩陣力的乘積。由此派生出解方程組的宜接分解法。由高斯消元法得到啟發(fā)，對E消元的過程相當(dāng)于將工分解為一個上三角矩陣和一個下三角矩陣的過程。如果直接分解W得到上和二'，A-LUQ這時方程癡化為,令由。解出了；再由改解出工。這就是直接分解法。將方陣工分解為4.AU,當(dāng)A是單位下三角矩陣，二.是上三角矩陣時，稱為多利特爾（Doclittle）分解；當(dāng)上是下三角矩陣，二是單位上三角矩陣時，稱為庫朗（Courant）分解。分解的條件是若方陣三的各階主子式不為零，則多利特爾分解或庫朗分解是可行和惟一的。多利特爾分解多利特爾分解步驟若W的各階主子式不為零，三可分解為A其中心為單位下三角矩陣，二為上三角矩陣，即矩陣上和下共有『個未知元素，按照二的行元素上的列元素的順序，對每個,列出式（3.16）兩邊對應(yīng)的元素關(guān)系式，一個關(guān)系式解出一個工或二’的元素。下面是計算上和二’的元素的步驟：（1）計算二的第一行元素要計算則列出式（3.20）等號兩邊的第1行第1列元素的關(guān)系式:.1?各也＜■(10 0故一般地，由了的第一行元素的關(guān)系式Jav?■(10 0)得到

Vv-OqJ-1,2.--.*(2)計算上的第一列元素?Ju要計算‘二】，則列出式(3.20)等號兩邊的第2行第1列元素的關(guān)系式:.1?！?T故一般地，由￡的第1列元素的關(guān)系式得到ifl-Ofl/Ujpt-N….“(3)計算-的第2行元素(4)計算上的第2列元素若已算出二”的前￡一】行，上的前左-1列的元素，則(5)計算二的第2二行元素*里*由二的第2行元素的關(guān)系式:

、■江，7皿41…,3是上三角矩陣，"列標(biāo)二二行標(biāo)上。* A 1-4%N以0 ?2以f?-1 ，■】 ，」得到口(3.21)氣■、?L??9m

z(3.21)(6)計算L的第C列元素j匕尸‘電由￡的第2列元素的關(guān)系式：、?金5■4…41Tla…⑼藍(lán)o,:L是下三角矩陣，二行標(biāo)標(biāo)：-行標(biāo):二o得到一直做到上的第,7列，二的第工行為止。用4U直接分解方法求解方程組所需要的計算量仍為3 ,和用高斯消元法的計算量基本相同?？梢钥吹皆诜纸庵械拿總€元素只在式(3.21)或(3.22)中做而且僅做一次貢獻(xiàn)，如果需要節(jié)省空間，可將二以及上的元素直接放在矩陣占相應(yīng)元素的位置上。用直接分解法解方程4*首先作出分解/■幺“令解方程組“由于上是單位下三角矩陣，容易得到Mk (3.23)再解方程組及■，，其中時工進(jìn)行4U分解時，并不涉及常數(shù)項因此，若需要解具有相同系數(shù)的?系列線性方程組時，使用直接分解法可以達(dá)到事半功倍的效果。多利特爾直接分解算法1.輸入：方程組階數(shù)3系數(shù)矩陣-二和常數(shù)項2FORk=lTOn{FORj4TOn//計算:的第七行元素FORi=kMTO?〃計算￡的第士列元素完成分解////解方程組0-6////解方程組0-6//解方程組3.FORI-ITOn4FORi-ITOn5.輸出方程組的解巧/例3.2用多利特爾分解求解方程組:對*-1,有?v-?vJ-U3?n-Z?R-l?u-14??，2?23.-1/2-0.S^m-1/2-0.5對比-2,有?? -3-05-15??"??-的口-2-05-1.5外-圓-WiI)g-(2-05)/25-0.6對*-3,有?fT—?Q6.』?4力?，力?06解皿?，，即-Lq-l3.2.2追趕法很多科學(xué)計算問題中，常常所要求解的方程組為三對角方程組其中(3.26)并且滿足條件:Ml*#。lARk1*1*11^，。八1.2.…l4H，po稱a為對角占優(yōu)的三對角線矩陣。其特征是非零元素僅分布在對角線及對角線兩側(cè)的位置。在樣條插值函數(shù)的“關(guān)系式中就出現(xiàn)過這類矩陣。事實上，許多連續(xù)問題經(jīng)離散化后得到的線性方程組，其系數(shù)矩陣就是三對角或五對角形式的對角矩陣。利用矩陣直接分解法，求解具有三對角線系數(shù)矩陣的線性方程組十分簡單而有效。現(xiàn)將工分解為下三角矩陣￡,及單位上三角矩陣）的乘積。即‘其中利用矩陣乘法公式，比較人.工?“兩邊元素，可得到u,G=N"=2.3「避-1于是有!4.巧.嗎.令，嗎吟一q,_1/-2,W…,“(3.28)H叫/-L2,…/-I

(3.28)由些可見將三分解為￡及了，只需計算值｝及｛呼｝兩組數(shù)，然后解37,計算公式為:(3.29)M.工"Of-<W.i)/-2.3.--.M(3.29)再解“■’則得x.-XA-x.-XA-Xf4**m-L…21(3.30)整個求解過程是先由(3.28)及(3.29)求MLM)及E),這時i=12“工是，，追，，的過程，再由(3.24)求出｛&L這時1=?，…/是往回趕，故求解方程組(3.25)的整個過程稱為追趕法。3.2.3對稱矩陣的皿f分解對稱正定矩陣也是很多物理問題產(chǎn)生的?類矩陣，正定矩陣的各階主子式大于零?？梢宰C明，若三正定，則存在下三角矩陣工，使/?￡￡「，直接分解/的分解方法，稱為平方根法。由于在平方根法中含有計算平方根的步驟，因此很少采用平方根法的分解形式。對于對稱矩陣，常用的是人?￡。上1'分解。

由不對稱正定，可得% 令由d的對?稱性和分解的惟一性可證￡是對角元素為1的單位下三角矩陣。對矩陣工作多利特爾或庫朗分解，共計算爐個矩陣元素；對稱矩陣的UE分解，只需計算2個元素，減少了近一半的工作量。借助于多利特爾或庫朗分解計算公式，容易得到分解計算公式。設(shè)t有多利特爾分解形式：]Pi也-小?山-?；?.心；呼其中4在分解可套用多利特爾分解公式，只要計算下三角矩陣上和總的對角元素d*。計算中只需保存,"3)的元素，三的「行/列的元素用工的L表示。由于對稱正定矩陣的各階主子式大于零，直接調(diào)用多利特爾或庫朗分解公式可完3d分解計算，而不必借助于列主元的分解算法。對于1.1?2?一?4.有_J_J4?%-2初F-IM由UDlT^x~b，解方程組AxF可分三步完成：(1)解方程組&.，其中z…-加―(2)解方程組。其中,.dx。尤?〃4, 1-1,2(3)解方程局再**_習(xí)j ….1(3.32)JL4…,〃(3.33)(3.34)…” (3.35)(3.36)對稱矩陣的LDU分解算法.輸入：方程組階數(shù)七，系數(shù)矩陣上和常數(shù)項X.FORk=1TOnd*：=oSdd：r-1FORi.HcMTOn.略去解方程組步驟從矩陣分解角度看，直接分解法與消元法本質(zhì)上沒有多大區(qū)別，但實際計算時它們各有所長。一般來說，如果僅用單字長進(jìn)行計算，列主元消元法具有運(yùn)算量較少、精度高的優(yōu)點(diǎn)，故是常用的方法。但是，為了提高精度往往采取單字長數(shù)雙倍內(nèi)積的辦法（即作向量內(nèi)積計算時，采用雙倍加法，最終結(jié)果再舍入成單字長數(shù)），這時直接分解法能獲得較高的精度。3.4用LDC分解求解方程組：J28t.r1 o o|p'TOC\o"1-5"\h\z-1 I 0.0- 21-0.5I 39 \由有Lx.仇z-?-4冉?1,xmytx?-\2Li3.3范數(shù)向量范數(shù)在一維空間中，實軸上任意兩點(diǎn)41*距離用兩點(diǎn)差的絕對值|*b|表示。絕對值是一種度量形式的定義。范數(shù)是對函數(shù)、向量和矩陣定義的一種度量形式。任何對象的范數(shù)值都是一個非負(fù)實數(shù)。使用范數(shù)可以測量兩個函數(shù)、向量或矩陣之間的距離。向量范數(shù)是度量向量長度的一種定義形式。范數(shù)有多種定義形式，只要滿足下面的三個條件即可定義為一個范數(shù)。同一向量，采用不同的范數(shù)定義，可得到不同的范數(shù)值。定義3.1對任一向量/按照一個規(guī)則確定一個實數(shù)與它對應(yīng)，記該實數(shù)記為帆1,若滿足卜面上個性質(zhì)：⑴有P1W當(dāng)且僅當(dāng)*?。時，（非負(fù)性）⑵Uek,aw虐，有（齊次性）」（3.37）⑶“心肥，有產(chǎn)?平陽中（三角不等式）那么稱該頭數(shù)L■為向量上的氾數(shù)。幾個常用向量范數(shù)向量的J范數(shù)定義為巴?序if其中，經(jīng)常使用的是P?L2.三種向量范數(shù)。雙■力訃附+NI++M?】wl,標(biāo)+x：或?qū)懗蒻?國WL?靖{Id?g{kL同'>kO例3.5計算向量*?a?M_5f.P.L28的三種范數(shù)。設(shè)卜1*，*5-979161|X|L-max{l.X|-5|-5向量范數(shù)的等價性上兩種不同的范數(shù)定有限維線性空間片中任意向量范數(shù)的定義都是等價的。若4⑶出⑶是上兩種不同的范數(shù)定義，則必存在0<*<*《。,使￥2€肥均有

（證明略）向量的極限有了向量范數(shù)的定義，也就有了度量向量距離的標(biāo)準(zhǔn)，即可定義向量的極限和收斂概念了。設(shè)為X上向量序列，若存在向量使!m-十°,則稱向量列心是收斂的（H是某種向量范數(shù)），儀稱為該向量序列的極限。由向量范數(shù)的等價知，向量序列是否收斂與選取哪種范數(shù)無關(guān)。向量序列m.（物‘守""產(chǎn)），m=12…收斂的充分必要條件為其序列的每個分量收斂,即出者B存在。若媽/■角，則” …J就是向量序列（榕卜｛（產(chǎn)，近—y）的極限。解：算出每個向量分量的極限后得取格,為極限在計算方法中，計算的向量序列都是數(shù)據(jù)序列，當(dāng)pL-ei小于給定精度時,向量*二取格,為極限3.3.2矩陣范數(shù)矩陣范數(shù)定義定義3.2如果矩陣,.（哂電廣的某個非負(fù)實函數(shù)*3,記作風(fēng)，滿足條件：（1）Ml?。,當(dāng)且僅當(dāng)人?。時，°（非負(fù)性）⑵M_HM?tfe*（齊次性）⑶對于任意兩個階數(shù)相同的矩陣a?b有（三角不等式性）（4）46矩陣為同階矩陣（相容性）則稱M為矩陣范數(shù)。矩陣的算子范數(shù)常用的矩陣范數(shù)是矩陣的算子范數(shù)，可用向量范數(shù)定義：設(shè)AcL.Xe**,記方陣上的范數(shù)為阿，那么皆或 3（3.38）稱為矩陣的算子范數(shù)或從屬范數(shù)。這樣定義的矩陣范數(shù)滿足矩陣范數(shù)的所有性質(zhì)外，還滿足相容性:W為華階矩陣，丫￡wR■.恒有卬BMai （3.39）根據(jù)定義，對任一種從屬范數(shù)有即單位矩陣的范數(shù)是1。常用矩陣范數(shù)向量有三種常用范數(shù)，相對應(yīng)的矩陣范數(shù)的三種形式為：同,陶加《的行范數(shù)）（3.40）

（正的列范數(shù)）(3.41)（正的列范數(shù)）(3.41)（為是"4的最大特征值）（4的2范數(shù)）（3.42）證明：既然矩陣的算子范數(shù)是配上滿足牘I.〕向量范數(shù)wi的上確界，那么，找到這個上確界也就找到了矩陣的范數(shù)。⑴任取萬名必.設(shè)|?1,則血卜:MI/-4血卜:MI/-4喀快胴《卜加）副小踹加即14,踹加另一方面設(shè)極大值在上列達(dá)到，即端UI歲喇取入7。,0「?.0」.0「-，09,￥除第2分量為1外，其余分量均為0,于是有W卜隆.冬引r[?gjoj■能卻J由于H",故因此有

⑵任取二…則皿■二皿工2卜播熱舊卜即ML，福第J另一方面設(shè)極大值在上行達(dá)到，取”??…產(chǎn)叫了這里{l.akO-l.a<0于是M辦I故心㈱割(3)，)為對稱非負(fù)矩陣，具有非負(fù)特征值，并具有v個相互正交的單位特征向量。設(shè)*'A的特征值為4 2— 相應(yīng)的特征向量為足，其中聘為相互正交的單位向量。設(shè)XrgfAT*并旦Ml即不"，則j^AX■ +乙巧當(dāng)?…+斗!4"I"(山幽?，找㈤■之常j為4’-A1 M即對任意均有卬L"JZ,故于是如果A是對稱矩陣，那么,4?44,設(shè)W的特征值是4,則有14?°?口-皿=卜1還有一種與向量范數(shù)1rL相容的矩陣范數(shù)，稱為歐幾里得（Euclid）范數(shù)或舒爾（Schur）范數(shù)，用Ml表示，其定義為因為歐幾里得范數(shù)易于計算，在實用中是一種十分有用的范數(shù)。但它不能從屬于任何一種范數(shù)，因為與向量范數(shù)的等價性質(zhì)類似，矩陣范數(shù)之間也是等價的。?KKKM例3.7I,I,分別有 。解：-mae(卜1|?3,2?17)■9pyL-皿吠卜1|+2?3*7)70的特征值4■依電4?281Pl-^019Ml-（1+4+9+49）%-屈矩陣范數(shù)等價性定理對l上的任兩種范數(shù)Hi及H,存在常數(shù)。了三使矩陣特征值與范數(shù)關(guān)系若W是矩陣三的特征值（即存在非零向量者使得：AK=OX）,對任一算子范數(shù)有（相容性）.,MNMII即矩陣特征值的模不大于矩陣的任一范數(shù)。譜半徑若3有特征值卬?記。（&-慵"1"則稱。（'）為上的譜半徑。有了譜半徑的定義，矩陣的2范數(shù)可記為：風(fēng)■山（，'）。譜半徑與矩陣范數(shù)關(guān)系由矩陣譜半徑定義，可得到矩陣范數(shù)的另一重要性質(zhì)，。（⑷“ML在解方程組時，我們總是假定系數(shù)矩陣片和常數(shù)項上是準(zhǔn)確的，而在實際問題中，系數(shù)矩陣A和常數(shù)項士往往是從前面的近似計算所得，元素的誤差是不可避免的。這些誤差會對方程組小.5的解2產(chǎn)生多大的影響？矩陣的條件數(shù)將給出一種粗略的衡量尺度。定義3.3若三非奇異稱為三的條件數(shù)。其中■表示矩陣的某種范數(shù)。用矩陣V及其逆矩陣的范數(shù)的乘積表示矩陣的條件數(shù)，由于矩陣范數(shù)的定義不同，因而其條件數(shù)也不同，但是由于矩陣范數(shù)的等價性，故在不同范數(shù)下的條件數(shù)也是等價的。矩陣條件數(shù)的大小是衡量矩陣“壞”或“好”的標(biāo)志。對于線性方程組若系數(shù)矩陣有小擾動W這時方程組的解也有擾動，于是\-Cottd"+3A）（X+6工）?8。34對，*的影響可表示為：\-Cottd(3.44)若常數(shù)上有小擾動.，其對$方的影響表示為:(3.45)(3.45)因此，C*（內(nèi)大的矩陣稱為“壞矩陣''或"病態(tài)矩陣對于式4）大的矩陣，小的誤差可能會引起解的失真。一般說來，若上的按模最大特征值與按模最小特征值之比值較大時，矩陣就會呈病態(tài)。特別地，當(dāng)修4很小時，/總是病態(tài)的。例3.8方程組2160^+01441X,-01440H2969jq+0耐呢-08M2方程組的解為對常數(shù)項上引入擾動熱9SOOEMOgT,則解為化…]串”1]司[-0.48701即取?(-1.009,1.5130);雖然.很小，但解已完全失真。detX--Idr,^-1^(°8648…][-1,29690.2161)14,=21617卜、L=L5l3xl(f故.&2TO02lx|”,4的病態(tài)是顯而易見的。.5病態(tài)方程組的解法如果系數(shù)矩陣W的條件數(shù)遠(yuǎn)大于1,則&T 為病態(tài)方程組，對病態(tài)方程組求解可采用以下措施：(1)采用高精度運(yùn)算，減輕病態(tài)影響，例如用雙倍字長運(yùn)算；⑵用預(yù)處理方法改善后的條件數(shù)，即選擇非奇異陣RQeL,使尸等價，而".儂的條件數(shù)比工改善，而求無■陽的解￡?0"x?即才為原方程組的解，計算時可選擇P?Q為對角陣或三角陣:

(3)平衡方法，當(dāng)三中元素的數(shù)量級相差很大，可采用行平衡或列平衡的方法改善■的條件數(shù)。設(shè)加廿l非奇異,計門移令。-蛔于是求…等價于求&檢.。6,或法■瓦這時孤-”的條件數(shù)可得到改善，這就是行均衡法。例3.9給定方程組為A的條件數(shù)Cg〃).7()f,若用行均衡法可取D"箍趴曠，1),則平衡后的方程九為3/而.*4,用三位有效數(shù)字的列主元消元法求解得nmimV3.6程序示例3.6程序示例程序3.1用列主元高斯消元法求解方程組:算法描述輸入：方程組階數(shù)n,矩陣a及列向量b。分解過程：對充?L2…選擇q*?g{l%U%uLT%iQ交換第上行和第產(chǎn);行方程對….熊jn….〃對f-A+l…回代過程：對…」?解得■夫?色-石仙)―程序源碼IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII〃RT/KW?則主元離裔箱元相隔口?〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h>#include<stdio.h>#defineMAX-N20〃(xi,yi)的最大維數(shù)intmain(){intn;intij,k;intmi;doublemx,tmp;staticoublea[MAXN],[MAXN],b[MAXN],x[MAX_N];printf(44\ninputnvalue(dimofAx=b):");〃輸入Ax=b的維數(shù)scanR"％d'',&n);ifi(n>MAXN){printf(4*TheinputnislargerthanMAXN,pleaseredefinetheMAXN.\n");return1;}if(n<=0){printfi(tfcPleaseinputanumberbetween1and%d.\n",MAX_N);return1;(〃輸入Ax=b的A矩陣printf(4fcNowinputthematrixa(i,j),i,j=0,...,%d:\n”,n-1);for(i=0,i<n;j++)fbr(j=Ou<n;j-H-scanC%lf,&a[i][j])〃輸入b向量printR"Nowinputthematrixb(i),i=0,…,%d:\n”,n-1)；fbr(i=O；i<n;i++)scanf(u%If”,&b[i]);fbr(i=O；i<n-l;i++)fbr(j=i+l,mi=i,mx=fabs(a[i][i]);j<n；j++)〃找主元素mi=j;mx=fabs(a[j][i]);{if(i<mi) 〃交換兩行{tmp=b[i];b[i]=b[mi];b[mi]=tmp;for(j=l;j<n;j++)(tmp=a[i][j];a[i][j]=a[mi][j]；a[mi][j]=tmp;))for(j=i+l;j<n;j++) 〃高斯消元{tmp=-a[j][i]/a[i][i];b[j]+=b[i]*tmp;for=(k=i,k<n;k-H-)a|j][k]+=a[i][k]*tmp;})x[n-l]=b[n-l]/a[n-l][n-l]; 〃求解方程fbr(i=n-2;i>0;i—){x[i]=b[i];

x[i]-=a[i][j]*xU]；x[i]/=a[i][i];〃輸出)〃輸出printfpSlove…xj=\n")；fbr(i=0;i<n++)printf("％f\n,x[i「)；return0;)// EndofFile 計算實例用列主元高斯消元法水求解方程組2、ff-4

<勺+34Zr*■6r*2j^*2xi-5程序輸入輸出Inputnvalue(dimofAx=b):3Nowinputthematrixa(i,j),i,j=0,...211132122Nowinputthematrixb(i),i=0, 2:465Solve...x_i=1.0000001.0000001.000000程序3.2用庫朗分解公式求解方程組:算法描述輸入方程組階數(shù)巴矩陣工及列向量入將矩陣工分解為￡口，其中Ai￡-勺?,4】。對土■對，■**7…上j，?%-2kJ%-?r￡3NJH對i=l2…4對i=nja-I,…AJ程序源碼IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII//Purpose:庫朗分解解Ax=b//IIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h>#include<stdio.h>#defineMAXN20 〃(xi,y_i)的最大維數(shù)intmain(){intn;inti,j,k;intmi;doublemx,tmp;statedoudlea[MAXN][MAXN],b[MAXN],x[MAXN],y[MAXN];statedoubleI[MAXN][MAXN],u[MAXN][MAXN];printfC'nlnputnvalue(dimofAx=b):");〃輸入加小的維數(shù)scanfT%d”,&n);ifi(n>MAXN){printf(fc*TheinputnislargerthenMAXN,pleadefinetheMAXN.\n");return1;if(n<=0){printfi[fcTleaseinputanumberbetween1and%d.\n",MAX_N);retun1;〃輸入AxW>的才矩陣printff'Nowinputthematrixa(ij),ij=O...，%d.\n,,,n-l);fbr(i=O;i<n;i-H-)for(j=O;j<n;j-H-)scanC4%If;&a[i][j]);〃輸入b矩陣printfit^Nowinputthematrixb(i),i=0,...d:\n,,,n-l);fdr(i=O;i<n;i++)scanfi[44%If;&b[i]);for(i=0;i<n;i++)u[i][i]=l;//U矩陣對角元素為1for(k=0;k<n;k++){for(i=k;i<n;i++)〃計算上矩陣的第三列元素(l[i][k]=a[i][k];For(j=0;j<=k-lJ-H-,l[i][k]-=(l[i][j]*u]j][k])；)for(j=k+l；j<n;j++) 〃計算二?矩陣的第匕行元素(u[k][j]=a[k][j];fbr(i=O;i<=k-l;i++)u[k]U]-=(l[k]U]*u[i]UD；for(i=0;i<n;i++)//Ly=b計算y{y[i]=b[i]；for(j=0;i<=i-l;j++)y[i]-=(i[i]U]*yU]);y[i]/=l[i][i];(fbr(i=n-1;i>=0;i-)//Ux=y計算x{x[i]=y[i];for(j=i+l;j<n;j++)x[i]-=(u[i][j]*x[j]);)printf("Solve...x」=\n")；〃輸出for(i=0;i<=n;I++)printf(11%return0;}計算實例用庫朗分解公式求解方程組：程序輸入輸出Inputnvalue(dimofAx=b):3Nowinputthematrixa(ij)ij=O,...,2：121-2-1-50-16Nowinputthematrixb(i),i=0,...,2：24-6350solve…x-i=7.0000004.0000009.000000程序3.3用Sf分解求解方程組:算法描述輸入矩陣/及列向量5。將矩陣工分解為人?皿,，其中解小?A,先解￡Z?凱再由0■”得I-12,-,M最后由/x?尸得a4.乂?就"%，'rx-1…21程序源碼IIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIPurpose:LD6分解解方程組〃IIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#include<stdio.h>#include<math.h>#defineMAX-N20〃(x?i,y?i)的最大維數(shù)intmain(){intn;inti,j,k;ilntmi;doublemx,tmp;staticdoublea[MAX-N][MAX-N],B[MAX-N],x[MAX-N],y[MAX-N]