現(xiàn)代數(shù)值計(jì)算（第3版）課件 第8、9章 矩陣特征值與特征向量的計(jì)算、常微分方程初邊值問(wèn)題數(shù)值解

第八章AdvancedNumericalComputing矩陣特征值與特征向量的計(jì)算現(xiàn)代數(shù)值計(jì)算同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院目錄/Contents第八章矩陣特征值和特征向量的計(jì)算第一節(jié)特征值問(wèn)題的基本理論第二節(jié)冪方法第三節(jié)計(jì)算矩陣所有特征值的QR方法前言矩陣特征值問(wèn)題:給定

階方陣 ,尋找常數(shù)

和非零向量

,使得其中

和

分別表示復(fù)數(shù)域和

維復(fù)向量空間。前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法【定義8.1】對(duì)于

階方陣

,稱(chēng)

為矩陣

的特征多項(xiàng)式?！径x8.2】

是關(guān)于

的

次多項(xiàng)式,且

的零點(diǎn)即為矩陣

的全

部特征值。【定義8.3】設(shè)

為矩陣

的特征值,由于 ,故方程

必有非零解 ,稱(chēng)

為矩陣

對(duì)應(yīng)于特征值

的特征向量。前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法特征多項(xiàng)式與特征方程線(xiàn)性代數(shù)中介紹方法：第一步

求解特征方程，得到矩陣

的所有特征值；第二步

求解齊次方程

，

得到

所對(duì)應(yīng)的特征向量。矩陣階數(shù)較大時(shí),特征多項(xiàng)式的零點(diǎn)沒(méi)有簡(jiǎn)單的解析表達(dá)式,只能通過(guò)近似計(jì)算得到。而高次方程近似求根的穩(wěn)定性差,求得的近似解會(huì)有較大誤差!【例子8.1】計(jì)算矩陣

的全部特征值。解：矩陣

的特征多項(xiàng)式為矩陣

的特征值為

的解,即

,

。前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法(2) ?！径ɡ?.1】若

是矩陣

的特征值,則有(1) ，其中

表示矩陣

的跡；【定義8.4】設(shè)

, 都是

階矩陣,若存在可逆矩陣

,使得

,則稱(chēng)矩陣

與矩陣

相似?！径ɡ?.2】若矩陣

與矩陣

相似,則

與

具有相同的特征值。前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法【定理8.3】Gerschgorin圓盤(pán)定理設(shè), ，則

的每一個(gè)特征值必屬于下述某個(gè)圓盤(pán)之中其中上式表示復(fù)平面上以

為中心,以

為半的圓盤(pán)。證：

是任意特征值,

是對(duì)應(yīng)的特征向量,則記 ,第

個(gè)方程為前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法【例8.2】對(duì)于矩陣

，圓盤(pán)定理中的圓盤(pán)為：注：前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法【定義8.5】設(shè)

為

階實(shí)對(duì)稱(chēng)矩陣,對(duì)于任意非零向量

,稱(chēng)

為關(guān)于向量

的瑞雷(Rayleigh)商,其

中

為

中兩向量

, 的內(nèi)積。【定理8.4】設(shè)

為

階實(shí)對(duì)稱(chēng)矩陣, 為其全部特征值,則(1)(2)前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法證：因?yàn)?/p>

實(shí)對(duì)稱(chēng)矩陣,故有完全的特征向量系。設(shè)

,

,

,是

的規(guī)范化正交特征向量組,即

。任一非零向量

都能

表示為 ,且故前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法故同理可證在瑞雷商中,分別取

和 ,可達(dá)到瑞雷商的最大值

和最小值,故前言----特征值問(wèn)題的線(xiàn)性代數(shù)解法目錄/Contents第八章矩陣特征值和特征向量的計(jì)算第一節(jié)特征值問(wèn)題的基本理論第二節(jié)冪方法第三節(jié)計(jì)算矩陣所有特征值的QR方法8.1冪方法8.1.1乘冪法乘冪法是一種計(jì)算矩陣按模最大特征值的方法,該方法簡(jiǎn)單有效并能給出對(duì)應(yīng)的特征向量。乘冪法需對(duì)實(shí)矩陣

做如下假設(shè)：(1)具有

個(gè)線(xiàn)性無(wú)關(guān)的特征向量 , ;(2) ,其中

是實(shí)特征值且滿(mǎn)足算法描述（未考慮停止迭代條件）選定初始向量

；利用關(guān)系式

迭代計(jì)算。(1)給定初始非零向量

及控制參數(shù);(2)(3)(4)對(duì) ,做

【算法8.1】原始乘冪法當(dāng)

時(shí),停止迭代;(5)返回近似按模最大特征值

和近似特征向量

.

為向量

中非零分量所對(duì)應(yīng)的指標(biāo)集,

為向量

中非零分量的個(gè)數(shù).8.2冪方法8.2.1乘冪法【例8.3】用乘冪法計(jì)算矩陣

的按模最大特征值及其對(duì)應(yīng)的

特征向量。解：取初始向量

和控制參數(shù)

。利用迭代關(guān)系式

進(jìn)行計(jì)算。8.2冪方法8.2.1乘冪法【例8.3】用乘冪法計(jì)算矩陣

的按模最大特征值及其對(duì)應(yīng)的

特征向量。解：8.2冪方法8.2.1乘冪法0(1.0,0.5,-0.5)----1(18,13,-13)23.3333--2(372,290,-290)21.76071.57263(7944,6292,-6292)21.58260.17804(170832,135752,-135752)21.55170.03095(3679008,2925520,-2925520)21.54560.00616(79254336,63031328,-63031328)21.54430.00137(1707427968,1357964608,-1357964608)21.54410.00038(36784710912,2925607232,-2925607232)21.5440<0.0001解：取按模最大特征值 ,對(duì)應(yīng)的特征向量為注：特征向量

的分量隨著

的增大而不斷增大！8.2冪方法8.2.1乘冪法【例8.3】用乘冪法計(jì)算矩陣

的按模最大特征值及其對(duì)應(yīng)的

特征向量。取初始向量

線(xiàn)性無(wú)關(guān),故利用

計(jì)算得向量序列對(duì) ,因?yàn)?/p>

,所以歸納法得8.2冪方法8.2.2算法分析由,當(dāng)

充分大時(shí),所以于是注意到則有這表明

可以近似地作為矩陣

按模最大特征值

所對(duì)應(yīng)的特征向量。若將向量

的第

個(gè)分量記為

,則或8.2冪方法8.2.2算法分析“上溢”和“下溢”現(xiàn)象從上例可以看出,隨著迭代步數(shù)

的增加,向量的分量逐漸增大。事實(shí)上,當(dāng)

或

時(shí),向量序列

中,各個(gè)不為零的分量將會(huì)隨著

的增大,而或無(wú)限增大或趨于零。規(guī)范化規(guī)范化的方法為用

的按模最大分量,記為

,去除所有分量,得到向量 .利用規(guī)范化的方法可將向量

的各個(gè)分量控制在

中.例如,若,

,則

,

.8.2冪方法8.2.2算法分析【算法8.2】改進(jìn)乘冪法8.2冪方法8.2.3改進(jìn)乘冪法(1)給定初始非零向量

及控制參數(shù)

;(2) 為向量

的按模最大分量;(3) ;(4)對(duì)

,做為向量

的按模最大分量;當(dāng)

時(shí),停止迭代;(5)返回近似按模最大特征值

和近似特征向量

。解：取初始向量

和控制參數(shù)

.

則

和

.【例8.4】用改進(jìn)乘冪法計(jì)算例2.1中矩陣

的按模最大特征值及其對(duì)應(yīng)的特征向量。8.2冪方法8.2.3改進(jìn)乘冪法解：【例8.4】用改進(jìn)乘冪法計(jì)算例2.1中矩陣

的按模最大特征值及其對(duì)應(yīng)的特征向量。8.2冪方法8.2.3改進(jìn)乘冪法0(1.0,0.5,-0.5)1.0000(1.0000,0.5000,-0.5000)1(18,13,-13)18.0000(1.0000,0.7222,-0.7222)2(20.6667,16.1111,-16.1111)20.6667(1.0000,0.7796,-0.7796)3(21.3548,16.9140,-16.914021.3548(1.0000,0.7920,-0.7920)4(21.5045,17.0886,-17.0886)21.5045(1.0000,0.7947,-0.7947)5(21.5358,17.1251,-17.1251)21.5358(1.0000,0.7952,-0.7952)6(21.5423,17.1327,-17.1327)21.5423(1.0000,0.7953,-0.7953)7(21.5437,17.1343,-17.1343)21.5437(1.0000,0.7953,-0.7953)8(21.543917.1346,-17.1346)21.5439(1.0000,0.7953,-0.7953)9(21.5440,17.1347,-17.1347)21.5440(1.0000,0.7953,-0.7953)解：當(dāng)

時(shí) ,迭代停止。

矩陣

的按模最大特征值 ,

對(duì)應(yīng)的特征向量 .8.2冪方法8.2.3改進(jìn)乘冪法【例8.4】用改進(jìn)乘冪法計(jì)算例2.1中矩陣

的按模最大特征值及其對(duì)應(yīng)的特征向量。取初始向量 ,由得

和則8.2冪方法8.2.4算法分析所以,由

,當(dāng)

充分大時(shí),8.2冪方法8.2.4算法分析從而則

可作為按模最大特征值

的近似, 可作為對(duì)應(yīng)的近似特征向量。8.2冪方法8.2.4算法分析function[t,y]=eigIPower(a,xinit,ep)

v0=xinit;

[tv,ti]=max(abs(v0));

lam0=v0(ti);

u0=v0/lam0;

flag=0;

while(flag==0)

v1=a*u0;

[tv,ti]=max(abs(v1));8.2冪方法8.2.5源程序

lam1=v1(ti);

u0=v1/lam1;

err=abs(lam0-lam1);

if(err<=ep)

flag=1;

end

lam0=lam1;

end

t=lam1;

y=u0;【例8.5】用改進(jìn)的乘冪法求下面矩陣的按模最大的特征值和對(duì)應(yīng)的特征向量,精度為

.8.2冪方法8.2冪方法【例8.5】用改進(jìn)的乘冪法求下面矩陣的按模最大的特征值和對(duì)應(yīng)的特征向量,精度為

.8.2冪方法【例8.5】用改進(jìn)的乘冪法求下面矩陣的按模最大的特征值和對(duì)應(yīng)的特征向量,精度為

.8.2冪方法【例8.5】用改進(jìn)的乘冪法求下面矩陣的按模最大的特征值和對(duì)應(yīng)的特征向量,精度為

.8.2冪方法【例8.5】用改進(jìn)的乘冪法求下面矩陣的按模最大的特征值和對(duì)應(yīng)的特征向量,精度為

.8.2冪方法【例8.5】用改進(jìn)的乘冪法求下面矩陣的按模最大的特征值和對(duì)應(yīng)的特征向量,精度為

.8.2.6反冪法反冪法可用來(lái)計(jì)算非奇異矩陣的按模最小特征值及其對(duì)應(yīng)的特征向量。(1)具有

個(gè)線(xiàn)性無(wú)關(guān)的特征向量 ;(2) ,其中

是實(shí)特征值且滿(mǎn)足

.反冪法需對(duì)實(shí)矩陣

作如下假設(shè):因?yàn)樗?/p>

的特征值滿(mǎn)足乘冪法求得矩陣

的按摸最大特征值及特征向量。得矩陣

的按摸最小特征值

及特征向量。8.2冪方法【算法8.3】反冪法(1)給定初始非零向量及控制參數(shù);(2) 為向量

的按模最大分量;(3) ;(4)對(duì)

,做為向量

的按模最大分量;當(dāng)

時(shí),停止迭代;(5)]返回近似按模最小特征值

和近似特征向量

.8.2.6反冪法8.2冪方法解：取初始向量

和控制參數(shù)

.

則

和

.為避免矩陣求逆,線(xiàn)性方程組

可以采用三角分解法求解。由于矩陣三角分解僅需計(jì)算一次,所以計(jì)算量可大為減少?！纠?.6】用反冪法計(jì)算例2.1中矩陣

的按模最小特征值及其對(duì)應(yīng)的特征向量。8.2.6反冪法8.2冪方法因此,矩陣

的按模最小特征值約為

,對(duì)應(yīng)的特征向量約為8.2.6反冪法8.2冪方法【例8.6】用反冪法計(jì)算例2.1中矩陣

的按模最小特征值及其對(duì)應(yīng)的特征向量。0(1.0,0.5,-0.5)1.0000(1.0000,-0.5000,0.5000)1(0.2083,-0.1250,0.1250)4.8000(1.0000,-0.6000,0.6000)2(0.2208,-0.1375,0.1375)4.5283(1.0000,-0.62260.6226)3(0.2237,-0.1403,0.1403)4.4710(1.0000,-0.6274,0.6274)4(0.2243,-0.1409,0.1409)4.4591(1.0000,-0.6284,0.6284)5(0.2244,-0.1411,0.1411)4.4566(1.0000,-0.6286,0.6286)6(0.2244,-0.1411,0.1411)4.4561(1.0000,-0.6287,0.6287)7(0.2244,-0.1411,0.1411)4.4560(1.0000,-0.6287,0.6287)8(0.2244,-0.1411,0.1411)4.4560(1.0000,-0.6287,0.6287)若已知矩陣

的某特征值的粗略近似值

,則可利用結(jié)合原點(diǎn)平移技巧的反冪法來(lái)計(jì)算該矩陣與數(shù)

最接近的特征值及其對(duì)應(yīng)的特征向量。首先作如下假設(shè):(1)數(shù)

是矩陣

的第

個(gè)特征值的近似,且滿(mǎn)足

由算法假設(shè)可知:矩陣

滿(mǎn)足反冪法所要求的假設(shè),且

是矩陣

的按摸最小特征值。對(duì)

應(yīng)用反冪法求出及其對(duì)應(yīng)特征向量的近似值,而矩陣

與數(shù)

最接近的特征值,其對(duì)應(yīng)的特征向量為.8.2.7結(jié)合原點(diǎn)平移的反冪法8.2冪方法【算法8.4】帶原點(diǎn)平移的反冪法(1)給定初始非零向量 ,控制參數(shù)以及數(shù)

;(2)為中按模最大的那個(gè)分量;(3) ;(4)對(duì) ,做為向量的按模最大分量;當(dāng)

時(shí),停止迭代;(5)返回近似特征值

和近似特征向量.8.2.7結(jié)合原點(diǎn)平移的反冪法8.2冪方法解：對(duì)矩陣

進(jìn)行LU分解,有取初始向量 ,利用結(jié)合原點(diǎn)平移的反冪法進(jìn)行計(jì)算。8.2.7結(jié)合原點(diǎn)平移的反冪法8.2冪方法【例8.7】求矩陣

與數(shù)

最接近的特征值及其對(duì)應(yīng)的特征向量。因此，矩陣與

最接近的特征值約為

,對(duì)應(yīng)的特征向量約為

.8.2.7結(jié)合原點(diǎn)平移的反冪法8.2冪方法【例8.7】求矩陣

與數(shù)

最接近的特征值及其對(duì)應(yīng)的特征向量。0(1,1,1)(1,1,1)-121(-2.4545,0.6667,04848)(1.0000,-0.2716,-0.1975)-13.40742(-4.5971,1.0782,0.7875)(1.0000,-0.2345,-0.1713)-13.21753(-4.5409,1.0676,0.7793)(1.0000,-0.2351,-0.1716)-13.22024(-4.5418,1.0678,0.7795)(1.0000,-0.2351,-0.1716)-13.2202目錄/Contents第八章矩陣特征值和特征向量的計(jì)算第一節(jié)特征值問(wèn)題的基本理論第二節(jié)冪方法第三節(jié)計(jì)算矩陣所有特征值的QR方法8.3

QR方法QR方法以矩陣的QR分解為基礎(chǔ),是計(jì)算中小型矩陣全部特征值的最有效方法之一,且具有收斂速度快,算法穩(wěn)定等特點(diǎn)。QR方法的算法描述如下:【算法8.5】QR方法(1)令

,并給定控制參數(shù).(2)對(duì)

,做這里

表示由矩陣

對(duì)角線(xiàn)元素構(gòu)成的列向量。當(dāng)

時(shí),停止迭代并返回

。計(jì)算矩陣的QR分解

;計(jì)算

;解：令

,利用QR分解算法對(duì)矩陣

進(jìn)行QR分解。8.3

QR方法【例8.8】用QR方法計(jì)算矩陣

的全部特征值.解：令

,利用QR分解算法對(duì)矩陣進(jìn)行QR分解。8.3

QR方法【例8.8】用QR方法計(jì)算矩陣

的全部特征值.解：迭代20次后得到8.3

QR方法【例8.8】用QR方法計(jì)算矩陣

的全部特征值.解：所以

的兩個(gè)特征值為

和 .其它兩個(gè)特征值滿(mǎn)足方程即

具有兩個(gè)共軛復(fù)特征值8.3

QR方法事實(shí)上,矩陣

的四個(gè)特征值分別為

.【例8.8】用QR方法計(jì)算矩陣

的全部特征值.Matlab簡(jiǎn)介：eig函數(shù)對(duì)于矩陣特征值問(wèn)題,Matlab軟件提供了eig函數(shù)。其格式如下:其中對(duì)角陣D表示矩陣A的全部特征值,矩陣V的第

列為D中第

列對(duì)角線(xiàn)元素所對(duì)應(yīng)的特征向量。示例：利用eig函數(shù)計(jì)算矩陣

的全部特征值及其對(duì)應(yīng)的特征向量.>>[V,D]=eig(A)可在Matlab軟件的命令窗口輸入如下命令：>>[V,D]=eig([5,-2,-5,-1;1,0,-3,2;0,2,2,-3;0,0,1,-2])Matlab簡(jiǎn)介：eig函數(shù)V=

0.9787

0.6602

0.6602

0.5774

0.1596

0.2330+0.3884i

0.2330-0.3884i

-0.0000

0.1277

0.4272-0.3884i

0.4272+0.3884i

0.5774

0.0213

0.0388-0.1553i

0.0388+0.1553i

0.5774D=

4.0000

0

0

0

0

1.0000+2.0000i

0

0

0

0

1.0000-2.0000i

0

0

0

0

-1.0000>>在命令窗口可看到計(jì)算結(jié)果：Matlab簡(jiǎn)介：eig函數(shù)>>[V,D]=eig([5,-2,-5,-1;1,0,-3,2;0,2,2,-3;0,0,1,-2])AdvancedNumericalComputing現(xiàn)代數(shù)值計(jì)算同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院學(xué)海無(wú)涯，祝你成功！第九章AdvancedNumericalComputing常微分方程初邊值問(wèn)題數(shù)值解現(xiàn)代數(shù)值計(jì)算同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院目錄/Contents第九章常微分方程初邊值問(wèn)題數(shù)值解第一節(jié)歐拉公式及其改進(jìn)第二節(jié)龍格-庫(kù)塔公式第三節(jié)收斂性與穩(wěn)定性第四節(jié)微分方程組和剛性問(wèn)題第五節(jié)有限差分法前言【定義9.1】含有未知函數(shù)及其導(dǎo)數(shù)的導(dǎo)數(shù)或微分的等式,常稱(chēng)為微分方程。一階常微分方程(組)的初值問(wèn)題:其中

是定義在上的

維函數(shù)向量,是定義在維區(qū)域上的

維已知函數(shù)向量.前言?xún)呻A常微分方程邊值問(wèn)題:其中

和

在區(qū)間

上連續(xù),.前言一階常微分方程(組)初值問(wèn)題:如果函數(shù)在區(qū)域中連續(xù),且關(guān)于滿(mǎn)足李普希茨(Lipschitz)條件,即對(duì)所有及

,

,總存在常數(shù)

,使得那么一階常微分方程(組)的初值問(wèn)題存在唯一解,且解連續(xù)依賴(lài)于初始條件和右端項(xiàng).兩階常微分方程邊值問(wèn)題：假設(shè)兩階常微分方程邊值問(wèn)題存在唯一解,且解連續(xù)依賴(lài)于邊界條件和右端項(xiàng).前言無(wú)論是初值問(wèn)題還是邊值問(wèn)題,其解

都是區(qū)間

上關(guān)于變量

的函數(shù)或函數(shù)向量.記

為求解區(qū)域中的一系列節(jié)點(diǎn).本章所討論的數(shù)值解指的是精確解

在這些節(jié)點(diǎn)

處的近似值

.簡(jiǎn)單起見(jiàn),假設(shè)點(diǎn)列滿(mǎn)足

其中

稱(chēng)為步長(zhǎng).9.1歐拉公式及其改進(jìn)9.1.1歐拉公式基本思想把方程中的導(dǎo)數(shù)項(xiàng)

用差商逼近,從而把一個(gè)微分方程轉(zhuǎn)化成為一個(gè)代數(shù)方程求解。在節(jié)點(diǎn)

處,

由于故將

的近似值

, (真解

是未知)代入上式就得到歐拉公式【算法9.1】歐拉方法(1)給定整數(shù)

;(2)(3)(4)對(duì)

做(5)返回向量

.9.1歐拉公式及其改進(jìn)9.1.1歐拉公式function[x,y]=odeeuler(f,y0,a,b,n)

y(1)=y0;

h=(b-a)/n;

x=a:h:b;

fori=1:n,

y(i+1)=y(i)+h*feval(f,x(i),y(i));

end歐拉方法的Matlab程序如下:9.1歐拉公式及其改進(jìn)9.1.1歐拉公式【例9.1】取步長(zhǎng)

，用歐拉方法求解如下初值問(wèn)題在區(qū)間

上的數(shù)值解.

(解析解為 )解：令. 由歐拉公式得到具體計(jì)算公式9.1歐拉公式及其改進(jìn)9.1.1歐拉公式9.1歐拉公式及其改進(jìn)9.1.1歐拉公式【例9.1】取步長(zhǎng)

，用歐拉方法求解如下初值問(wèn)題在區(qū)間

上的數(shù)值解.

(解析解為 )9.1歐拉公式及其改進(jìn)9.1.1歐拉公式【例9.1】取步長(zhǎng)

，用歐拉方法求解如下初值問(wèn)題在區(qū)間

上的數(shù)值解.

(解析解為 )0.00.10.20.30.40.51.00001.10001.19181.27741.35821.43511.00001.09541.18321.26491.34161.41420.00000.00460.00860.01250.01660.02090.60.70.80.91.01.50901.58031.64981.71781.78481.48321.54921.61251.67331.73210.02570.03110.03730.04450.0527

一階常微分方程初值問(wèn)題的真解

是

平面上過(guò)點(diǎn)

的一條特殊積分曲線(xiàn).利用歐拉公式得到數(shù)值解是從

出發(fā)的折線(xiàn)

,

如圖所示.所以歐拉方法又稱(chēng)為歐拉折線(xiàn)法.注意：誤差在不斷的積累和傳播！9.1歐拉公式及其改進(jìn)9.1.2歐拉方法的幾何解釋常微分方程

可改寫(xiě)為兩邊關(guān)于

在區(qū)間

上積分,得要得到

的值,就必須計(jì)算右端的積分選用不同的數(shù)值積分公式,就會(huì)得到不同的計(jì)算公式.這種將初值問(wèn)題離散化的方法稱(chēng)為求解初值問(wèn)題的數(shù)值積分方法.9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法左矩形公式代入得用

代替上式右端的

,得

的近似值 ,

這就是歐拉公式.9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法右矩形公式得后退歐拉公式9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法梯形求積公式得梯形公式注：由于梯形求積公式比矩形求積公式的代數(shù)精度高,因此一般情況下梯形公式的精度比歐拉公式高.9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法【算法9.2】梯形方法(1)給定整數(shù);(2)(3)(4)對(duì)

做求解方程(5)返回向量

.9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法【例9.2】用梯形公式求解例2.1中的初值問(wèn)題.解：梯形公式的具體形式為利用非線(xiàn)性方程的簡(jiǎn)單的迭代法求解由壓縮映像定理可以證明:當(dāng)

充分小時(shí),該迭代法收斂于唯一不動(dòng)點(diǎn).9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法【例9.2】用梯形公式求解例2.1中的初值問(wèn)題.0.01.00001.00000.00000.11.09571.09540.00020.21.18361.18320.00040.31.26541.26490.00050.41.34231.34160.00070.51.41511.41420.00080.61.48431.48320.00100.71.55041.54920.00120.81.61391.61250.00150.91.67511.67330.00181.01.73411.73210.0021利用高階數(shù)值積分公式,還可構(gòu)造其他離散格式.例如,阿達(dá)姆斯-巴士福德(Adams-Bashforth)公式阿達(dá)姆斯-莫爾頓(Adams-Moulton)公式9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法單步法與多步法，顯格式與隱格式歐拉公式與梯形公式僅含有

和

,只要知道

即可求出

，這類(lèi)公式稱(chēng)為單步法.阿達(dá)姆斯-巴士福德公式與阿達(dá)姆斯-莫爾頓公式計(jì)算時(shí),除了要知道

外,還需知道

,

等,這類(lèi)公式稱(chēng)為多步法.歐拉公式和阿達(dá)姆斯-巴士福德公式右端不含有

,把右端的已知量算出來(lái)即可求出

,這類(lèi)方法稱(chēng)為顯式格式.梯形公式與阿達(dá)姆斯-莫爾頓公式的右端隱含,求時(shí)通常需要解方程,此類(lèi)公式稱(chēng)為隱式格式.特別當(dāng)

關(guān)于

為非線(xiàn)性時(shí),就需要解非線(xiàn)性方程.9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法局部截?cái)嗾`差可用于表征求解初值問(wèn)題的數(shù)值方法的計(jì)算精度.【定義9.2】一般地,常微分方程初值問(wèn)題的數(shù)值解滿(mǎn)足形如的等式,其中為在

個(gè)節(jié)點(diǎn)處的數(shù)值解.假設(shè)上式中等號(hào)右端的數(shù)值解都是準(zhǔn)確的,則利用該式得到的數(shù)值解與真解之間的差異被稱(chēng)為數(shù)值方法的局部截?cái)嗾`差.9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法【定義9.3】若初值問(wèn)題離散格式的局部截?cái)嗾`差為

,則稱(chēng)該離散格式公式所代表的數(shù)值方法具有

階精度或稱(chēng)為

階方法.【例9.3】歐拉公式是一階方法.證：設(shè)

在

上充分光滑,假設(shè) ,則9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法證：設(shè)

在

上充分光滑,由表明歐拉公式的局部截?cái)嗾`差是

的同階無(wú)窮小,即歐拉公式是一階方法.9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法【定義9.3】若初值問(wèn)題離散格式的局部截?cái)嗾`差為

,則稱(chēng)該離散格式公式所代表的數(shù)值方法具有

階精度或稱(chēng)為

階方法.【例9.3】歐拉公式是一階方法.類(lèi)似可證梯形公式是二階方法,而阿達(dá)姆斯-巴士福德公式與阿達(dá)姆斯-莫爾頓公式均為四階方法.9.1歐拉公式及其改進(jìn)9.1.3數(shù)值積分與多步法【定義9.3】若初值問(wèn)題離散格式的局部截?cái)嗾`差為

,則稱(chēng)該離散格式公式所代表的數(shù)值方法具有

階精度或稱(chēng)為

階方法.【例9.3】歐拉公式是一階方法.對(duì)比例2.1和例2.2的計(jì)算結(jié)果:【定義9.4】預(yù)估校正方法為了避免隱格式的迭代求解,一種有效簡(jiǎn)化計(jì)算的方法是:當(dāng)

較小時(shí),先用顯格式計(jì)算一個(gè)合適的預(yù)估值,然后利用隱格式迭代一、二次得到校正值.預(yù)估校正這就是一個(gè)預(yù)估校正公式，上式通常稱(chēng)為改進(jìn)的歐拉公式.梯形方法相對(duì)于歐拉方法雖然精度提高了,但是每一步都要用迭代法解方程.每迭代一次,都要重新計(jì)算函數(shù)

的值,所以計(jì)算量也增加很多.9.1歐拉公式及其改進(jìn)9.1.4預(yù)估校正格式以下是改進(jìn)歐拉方法的算法描述:【算法9.4】改進(jìn)歐拉方法(1)給定整數(shù) ;(2)

;(3);(4)對(duì)

做(5)返回向量9.1歐拉公式及其改進(jìn)9.1.4預(yù)估校正格式function[x,y]=odeIEuler(f,y0,a,b,n)

y(1)=y0;

h=(b-a)/n;

x=a:h:b;

fori=1:n,

yp=y(i)+h*feval(f,x(i),y(i));

yc=y(i)+h*feval(f,x(i+1),yp);

y(i+1)=0.5*(yp+yc);

end改進(jìn)歐拉方法的Matlab程序如下:9.1歐拉公式及其改進(jìn)9.1.4預(yù)估校正格式解：取步長(zhǎng)由改進(jìn)的歐拉公式得.9.1歐拉公式及其改進(jìn)9.1.4預(yù)估校正格式【例9.4】用改進(jìn)的歐拉公式求解例2.1中的初值問(wèn)題.【例9.4】用改進(jìn)的歐拉公式求解例2.1中的初值問(wèn)題.9.1歐拉公式及其改進(jìn)9.1.4預(yù)估校正格式0.01.00001.00000.00000.11.09591.09540.00050.21.18411.18320.00090.31.26621.26490.00130.41.34341.34160.00170.51.41641.41420.00220.61.48601.48320.00270.71.55251.54920.00330.81.61651.61250.00400.91.67821.67330.00481.01.73791.73210.0058簡(jiǎn)單的預(yù)估校正格式都包含兩個(gè)計(jì)算公式：

一個(gè)是顯式公式,作為預(yù)估公式；

另一個(gè)是隱式公式作為校正公式.當(dāng)然也可以構(gòu)造包含多個(gè)計(jì)算公式的預(yù)估校正格式.構(gòu)造預(yù)估校正格式時(shí),應(yīng)注意階數(shù)的匹配.例如在改進(jìn)的歐拉公式中,校正公式具有二階精度,而預(yù)估公式僅具有一階精度.因?yàn)樘峁┑念A(yù)估精度較差,且僅經(jīng)一次校正,校正值的精度不會(huì)太高.可改為下面的預(yù)估校正格式預(yù)估校正9.1歐拉公式及其改進(jìn)9.1.4預(yù)估校正格式四階阿達(dá)姆斯預(yù)估校正系統(tǒng)9.1歐拉公式及其改進(jìn)9.1.4預(yù)估校正格式目錄/Contents第九章常微分方程初邊值問(wèn)題數(shù)值解第一節(jié)歐拉公式及其改進(jìn)第二節(jié)龍格-庫(kù)塔公式第三節(jié)收斂性與穩(wěn)定性第四節(jié)微分方程組和剛性問(wèn)題第五節(jié)有限差分法9.2龍格-庫(kù)塔公式9.2.1基本思想Taylor展開(kāi)則最后有9.2龍格-庫(kù)塔公式9.2.1基本思想（

）先對(duì)歐拉公式和預(yù)估-校正公式作進(jìn)一步的分析:歐拉公式可改寫(xiě)成(1)計(jì)算

,需計(jì)算一次

的值;(2)若設(shè)

,

則

的表達(dá)式與的泰勒展開(kāi)式的前兩項(xiàng)

完全相同,即局部截?cái)嗾`差為

.9.2龍格-庫(kù)塔公式9.2.1基本思想預(yù)估校正公式可改寫(xiě)為(1)計(jì)算

,需計(jì)算兩次

的值;(2)若設(shè)

,則

的展開(kāi)式與的泰勒展開(kāi)式的前三項(xiàng)

完全相同,即局部截?cái)嗾`差為

.9.2龍格-庫(kù)塔公式9.2.1基本思想龍格-庫(kù)塔公式的基本思想：設(shè)法計(jì)算

在某些點(diǎn)上的函數(shù)值,然后對(duì)這些函數(shù)值作線(xiàn)性組合,構(gòu)造近似計(jì)算公式,再把近似公式和解的泰勒展開(kāi)式相比較,使前面的若干項(xiàng)吻合,從而獲得達(dá)到一定精度的數(shù)值計(jì)算公式.9.2龍格-庫(kù)塔公式9.2.1基本思想【定義9.5】一般的,顯式龍格-庫(kù)塔公式的形式為其中參數(shù)

,,是與初值問(wèn)題右端

和步長(zhǎng)均無(wú)關(guān)的常數(shù).上式稱(chēng)為

段的龍格-庫(kù)塔公式.特別地,若上式與

的泰勒展開(kāi)式的前

項(xiàng)完全一致,即局部截?cái)嗾`差達(dá)到

,則稱(chēng)為階

段龍格-庫(kù)塔公式.9.2龍格-庫(kù)塔公式9.2.1基本思想【例9.5】二階二段R-K方法9.2龍格-庫(kù)塔公式9.2.2二階龍格-庫(kù)塔公式預(yù)估校正公式中點(diǎn)公式9.2龍格-庫(kù)塔公式9.2.2二階龍格-庫(kù)塔公式二階二段龍格-庫(kù)塔算法：(1)給定整數(shù) ;(2)

;(3) ;(4)對(duì)

做

(5)返回向量

.9.2龍格-庫(kù)塔公式9.2.2二階龍格-庫(kù)塔公式二階二段龍格-庫(kù)塔Matlab程序：function[x,y]=rk2(dfun,a,b,y0,h)x=a:h:b;n=length(x);y(1)=y0;fork=2:n,

k1=h*feval(dfun,x(k-1),y(k-1));

k2=h*feval(dfun,x(k-1)+h/2,y(k-1)+k1/2);

y(k)=y(k-1)+1/2*(k1+k2);end9.2龍格-庫(kù)塔公式9.2.2二階龍格-庫(kù)塔公式高階龍格-庫(kù)塔公式的推導(dǎo)方法與上面討論的二階龍格-庫(kù)塔公式推導(dǎo)類(lèi)似.這里僅列舉幾個(gè)常用的公式:(1)三階三段龍格-庫(kù)塔公式(2)三階三段Heun公式9.2龍格-庫(kù)塔公式9.2.3高階龍格-庫(kù)塔公式(3)標(biāo)準(zhǔn)四階四段龍格-庫(kù)塔公式(4)四階四段Gill公式9.2龍格-庫(kù)塔公式9.2.3高階龍格-庫(kù)塔公式理論上,可構(gòu)造任意高階的計(jì)算方法,但方法的階數(shù)與計(jì)算函數(shù)值的次數(shù)之間的關(guān)系并非等量增加的,關(guān)系如下表:由于計(jì)算量較大,很少使用更高階的龍格-庫(kù)塔公式.事實(shí)上,對(duì)于大量的實(shí)際問(wèn)題,四階的龍格-庫(kù)塔公式已可滿(mǎn)足對(duì)精度要求.9.2龍格-庫(kù)塔公式9.2.3高階龍格-庫(kù)塔公式計(jì)算函數(shù)值次數(shù)r1234567階數(shù)p1234456【算法9.5】標(biāo)準(zhǔn)四階四段龍格-庫(kù)塔方法(1)給定整數(shù)

;(2) ;(3)

;(4)對(duì)

做(5)返回向量

.9.2龍格-庫(kù)塔公式9.2.3高階龍格-庫(kù)塔公式【例9.6】用標(biāo)準(zhǔn)四階四段龍格-庫(kù)塔方法求解例2.1中的初值問(wèn)題.解：已知 , ,并取

.

將

,代入標(biāo)準(zhǔn)的四階

四段龍格-庫(kù)塔格式,得具

體計(jì)算公式9.2龍格-庫(kù)塔公式9.2.3高階龍格-庫(kù)塔公式9.2龍格-庫(kù)塔公式9.2.3高階龍格-庫(kù)塔公式【例9.6】用標(biāo)準(zhǔn)四階四段龍格-庫(kù)塔方法求解例2.1中的初值問(wèn)題.0.01.000000001.000000000.0000e-50.11.095445531.095445120.0417e-50.21.183216741.183215960.0789e-50.31.264912231.264911060.1164e-50.41.341642541.341640790.1567e-50.51.414215581.414213560.2016e-50.61.483242221.483239690.2525e-50.71.591964521.549193340.3114e-50.81.612455351.612551550.3800e-50.91.673324661.673320050.4606e-51.01.732056371.732150810.5558e-5目錄/Contents第九章常微分方程初邊值問(wèn)題數(shù)值解第一節(jié)歐拉公式及其改進(jìn)第二節(jié)龍格-庫(kù)塔公式第三節(jié)收斂性與穩(wěn)定性第四節(jié)微分方程組和剛性問(wèn)題第五節(jié)有限差分法9.3收斂性與穩(wěn)定性9.3.1單步法的收斂性【定義9.6】若某數(shù)值方法對(duì)任意固定節(jié)點(diǎn)

,當(dāng)

時(shí),有

,則稱(chēng)該方法是收斂的.一般來(lái)說(shuō),顯式單步法總可寫(xiě)成如下形式:其中

稱(chēng)為增量函數(shù).增量函數(shù)依賴(lài)于

且僅是關(guān)于

,和

的連續(xù)函數(shù).例如歐拉公式的增量函數(shù)為改進(jìn)的歐拉公式的增量函數(shù)為下面的定理給出了一個(gè)判斷

階方法收斂性的方法,并指明了局部截?cái)嗾`差與收斂性的關(guān)系.【定理9.1】若

階單步法

中的增量函數(shù)

在區(qū)域

, , 上連續(xù),并且關(guān)于滿(mǎn)足李普希茨條件,即其中

為正常數(shù),則由該單步法得到的

滿(mǎn)足如下誤差估計(jì):9.3收斂性與穩(wěn)定性9.3.1單步法的收斂性證：令

,則有其中

表示局部截?cái)嗾`差.9.3收斂性與穩(wěn)定性9.3.1單步法的收斂性利用三角不等式,可得

階單步法的整體截?cái)嗾`差滿(mǎn)足如下估計(jì)式:9.3收斂性與穩(wěn)定性9.3.1單步法的收斂性注意到故定理得證.9.3收斂性與穩(wěn)定性9.3.1單步法的收斂性應(yīng)用上述定理可知:當(dāng)

時(shí),單步法是收斂的;當(dāng)局部截?cái)嗾`差

時(shí),其整體截?cái)嗾`差

.由于整體截?cái)嗾`差比局部截?cái)嗾`差低一階,因此在構(gòu)造高精度的計(jì)算方法時(shí),只要設(shè)法提高方法的局部截?cái)嗾`差.對(duì)于歐拉方法,其增量函數(shù)

,所以當(dāng)初值問(wèn)題的右端函數(shù)

滿(mǎn)足李普希茨條件時(shí),它是收斂的.9.3收斂性與穩(wěn)定性9.3.1單步法的收斂性對(duì)于標(biāo)準(zhǔn)四階四段龍格-庫(kù)塔公式,若

滿(mǎn)足李普希茨條件,則龍格-庫(kù)塔公式的增量函數(shù)可寫(xiě)成設(shè)李普希茨常數(shù)為

.因?yàn)? ,所以9.3收斂性與穩(wěn)定性9.3.1單步法的收斂性故因此,標(biāo)準(zhǔn)四階四段龍格-庫(kù)塔方法收斂.9.3收斂性與穩(wěn)定性9.3.1單步法的收斂性對(duì)于標(biāo)準(zhǔn)四階四段龍格-庫(kù)塔公式,若

滿(mǎn)足李普希茨條件,則龍格-庫(kù)塔公式的增量函數(shù)可寫(xiě)成設(shè)某數(shù)值方法在節(jié)點(diǎn)

處的數(shù)值解應(yīng)為,但實(shí)際計(jì)算時(shí)因?yàn)橛姓`差,算得的值為.記

,并稱(chēng)之為第

步數(shù)值解的擾動(dòng).9.3收斂性與穩(wěn)定性9.3.2單步法的穩(wěn)定性【定義9.7】若一種數(shù)值方法在節(jié)點(diǎn)處的數(shù)值解

伴隨有大小為

的擾動(dòng),由的傳播而導(dǎo)致以后各節(jié)點(diǎn)上的值

產(chǎn)生的擾動(dòng)值 ,其絕對(duì)值均不超過(guò),即則稱(chēng)該數(shù)值方法是絕對(duì)穩(wěn)定的.特別地,一個(gè)約定俗成的方法是對(duì)如下的“模型方程”來(lái)討論數(shù)值方法的絕對(duì)穩(wěn)定性.數(shù)值穩(wěn)定性的分析是相當(dāng)復(fù)雜的,且總跟微分方程的右端

有關(guān)，所以進(jìn)行一般性的討論常微分方程初值問(wèn)題的穩(wěn)定性很困難,一般只針對(duì)某一固定方程或某一類(lèi)方程來(lái)討論穩(wěn)定性.9.3收斂性與穩(wěn)定性9.3.2單步法的穩(wěn)定性歐拉公式的穩(wěn)定性:將歐拉公式用于模型方程可得其中, 設(shè)在計(jì)算時(shí)有擾動(dòng),由的傳播造成

產(chǎn)生擾動(dòng)值

,則所以,歐拉方法的絕對(duì)穩(wěn)定區(qū)域是復(fù)平面上以

為圓心,1為半徑的圓盤(pán).9.3收斂性與穩(wěn)定性9.3.2單步法的穩(wěn)定性故要?dú)W拉公式穩(wěn)定,只要后退歐拉公式的穩(wěn)定性:將后退歐拉公式用于模型方程可得即于是有為使后退歐拉公式穩(wěn)定,只要即

.所以,后退歐拉公式的絕對(duì)穩(wěn)定區(qū)域是復(fù)平面上以為圓心,1為半徑的圓外區(qū)域.顯然,它比歐拉方法的絕對(duì)穩(wěn)定區(qū)域要大很多.9.3收斂性與穩(wěn)定性9.3.2單步法的穩(wěn)定性標(biāo)準(zhǔn)四階四段龍格庫(kù)塔公式的絕對(duì)穩(wěn)定區(qū)域?yàn)槿?/p>

為實(shí)數(shù),則由上式可近似得到

.另外,當(dāng)是負(fù)實(shí)數(shù)時(shí),阿達(dá)姆斯-巴士福德公式的絕對(duì)穩(wěn)定區(qū)域?yàn)?/p>

,阿達(dá)姆斯-莫爾頓公式的絕對(duì)穩(wěn)定區(qū)域?yàn)?一般而言,隱式格式的穩(wěn)定性要比顯式格式的穩(wěn)定性要好!9.3收斂性與穩(wěn)定性9.3.2單步法的穩(wěn)定性目錄/Contents第九章常微分方程初邊值問(wèn)題數(shù)值解第一節(jié)歐拉公式及其改進(jìn)第二節(jié)龍格-庫(kù)塔公式第三節(jié)收斂性與穩(wěn)定性第四節(jié)微分方程組和剛性問(wèn)題第五節(jié)有限差分法9.4微分方程組和剛性問(wèn)題9.4.1一階常微分方程組前面介紹的方法大部分都可平行地推廣到方程組的情況.例如,考慮如下一階常微分方程組初值條件為將區(qū)間

進(jìn)行

等分,記

.將歐拉公式應(yīng)用于該問(wèn)題,則有且9.4微分方程組和剛性問(wèn)題9.4.1一階常微分方程組下面討論應(yīng)用歐拉方法求解上述常微分方程組的絕對(duì)穩(wěn)定性.簡(jiǎn)單計(jì)算可知其中

和

分別表示第

和步數(shù)值解的擾動(dòng).由于矩陣

的特征值為 ,故可以證明當(dāng)(或

)時(shí),即該方法是絕對(duì)穩(wěn)定的.9.4微分方程組和剛性問(wèn)題9.4.1一階常微分方程組對(duì)于高階常微分方程初值問(wèn)題若記則可將高階微分方程化為一階微分方程組,從而利用前面介紹的方法近似求解.9.4微分方程組和剛性問(wèn)題9.4.2高階常微分方程例如,考慮初值問(wèn)題令和

,則該初值問(wèn)題轉(zhuǎn)化成一階常微分方程組初值條件為9.4微分方程組和剛性問(wèn)題9.4.2高階常微分方程【定義9.8】常微分方程初值問(wèn)題（1.1）稱(chēng)為是剛性方程組,若對(duì)成立下列各式:(1)(2)其中

是雅可比矩陣

的特征值.

稱(chēng)為剛性比.9.4微分方程組和剛性問(wèn)題9.4.3剛性問(wèn)題例如,考慮常微分方程組初值條件為這就是一個(gè)剛性方程組,其雅可比矩陣為

,剛性比 .該問(wèn)題存在唯一解9.4微分方程組和剛性問(wèn)題9.4.3剛性問(wèn)題應(yīng)用標(biāo)準(zhǔn)四階龍格-庫(kù)塔公式求解該剛性問(wèn)題.由于受絕對(duì)穩(wěn)定性的限制,步長(zhǎng)

須滿(mǎn)足 ,其中為雅可比矩陣的特征值.這要求 .取和 ,用標(biāo)準(zhǔn)四階龍格-庫(kù)塔公式求解該剛性問(wèn)題.9.4微分方程組和剛性問(wèn)題9.4.3剛性問(wèn)題取,用標(biāo)準(zhǔn)四階龍格-庫(kù)塔公式求解該剛性問(wèn)題.9.4微分方程組和剛性問(wèn)題9.4.3剛性問(wèn)題0.11.793061-1.0320011.712219-0.87031520.21.423901-0.87468091.414070-0.85501480.31.131575-0.72499841.140523-0.72289100.40.9094086-0.60821410.9092763-0.60794750.50.7387877-0.51565750.7387506-0.51558100.60.6057094-0.44041080.6056833-0.44035580.70.4998603-0.37740380.4998361-0.37735400.80.4136714-0.32295350.4136490-0.32290780.90.3416143-0.27440880.3415939-0.27436731.00.2796748-0.22988770.2796568-0.2298511取 ,用標(biāo)準(zhǔn)四階龍格-庫(kù)塔公式求解該剛性問(wèn)題.9.4微分方程組和剛性問(wèn)題9.4.3剛性問(wèn)題0.11.793061-1.032001-2.6451697.8445270.21.423901-0.8746809-18.4515838.876310.31.131575-0.7249984-87.47221176.48280.40.9094086-0.6082141-934.0722789.35400.50.7387877-0.5156575-1760.0163520.000.60.6057094-0.4404108-7848.55015697.840.70.4998603-0.3774038-34989.6369979.870.80.4136714-0.3229535-155979.4311969.50.90.3416143-0.2744088-695332.013906641.00.2796748