任課教師：楊超電子郵箱:計(jì)算方法與程序設(shè)計(jì)2第8章常微分方程數(shù)值解法實(shí)際問(wèn)題中，大量的常微分方程很難用初等解法求出其精確解，因此促使人們?nèi)ヌ接懬蠼獬Ｎ⒎址匠痰臄?shù)值方法。

計(jì)算機(jī)的發(fā)展，為常微分方程的數(shù)值求解提供了便利。

本章主要介紹如下的一階常微分方程初值問(wèn)題的數(shù)值解法。3第8章常微分方程數(shù)值解法定義1：設(shè)函數(shù)f(t,u)定義在區(qū)域D內(nèi)，若存在正整數(shù)L，使得對(duì)任意，成立下列不等式稱(chēng)函數(shù)f(t,u)在D內(nèi)關(guān)于u滿足Lipschitz條件(L-條件)，L稱(chēng)為L(zhǎng)ipschitz常數(shù)。定理1：若函數(shù)f(t,u)在區(qū)域D內(nèi)滿足:(1)連續(xù)；(2)關(guān)于u滿足Lipschitz條件，則存在唯一的連續(xù)可微函數(shù)u=u(t)滿足及初始條件u(0)=u0，且連續(xù)地依賴于初始值u0。4第8章常微分方程數(shù)值解法所謂常微分方程數(shù)值解法就是求初值問(wèn)題(1)的解u(t)關(guān)于變量t在指定的離散點(diǎn)列(或稱(chēng)為節(jié)點(diǎn))上的近似值相鄰兩節(jié)點(diǎn)間的距離稱(chēng)為步長(zhǎng)，一般地，采取等距節(jié)點(diǎn)，即步長(zhǎng)hi=h為常數(shù)，有或5第8章常微分方程數(shù)值解法由于初值問(wèn)題中的初始條件u(0)=u0為已知，故可利用已知的來(lái)求下一節(jié)點(diǎn)t1處u(t1)的近似值u1，再利用u1來(lái)求u2，……，直到求出un為止。這種按節(jié)點(diǎn)的排列順序一步一步地向前推進(jìn)的方式求解的算法稱(chēng)為步進(jìn)式或遞推式算法。它是初值問(wèn)題數(shù)值解法的各種計(jì)算格式具有的共同特點(diǎn)。因此，只要能寫(xiě)出由前幾步已知信息u0，u1，…，un來(lái)計(jì)算un+1的遞推公式(通常為差分格式)，就可以完全表達(dá)該種算法。計(jì)算格式的構(gòu)造通常采用微商化為差商、Taylor級(jí)數(shù)展開(kāi)和數(shù)值積分等方法。6第8章常常微分方程數(shù)值解法Euler方法預(yù)估-校正法Runge-Kutta法線性多步法單步法的收斂性與穩(wěn)定性主要內(nèi)容：7第8章第1節(jié)Euler方法1.Euler格式

(1)數(shù)值積分法。將微分方程在區(qū)間[ti,ti+1]上積分，可得若在區(qū)間[ti,ti+1]上將f(t,u)近似看作常數(shù)f(ti,ui)，則有上式稱(chēng)為Euler格式。若u0已知，可依次求出u1,u2,…。8第8章第1節(jié)Euler方法

(2)數(shù)值微分法。設(shè)在區(qū)間[ti,ti+1]上，u

(t)可用線性插值函數(shù)近似，即將上式代入下式的左邊令，得Euler格式。9第8章第1節(jié)Euler方法

(3)級(jí)數(shù)展開(kāi)法。設(shè)微分方程的初值問(wèn)題的解將u(ti+1)在t=ti處展開(kāi)，可得忽略h的二次項(xiàng)，注意到du/dt=f(t,u)，得Euler格式Euler格式的幾何解釋?zhuān)河梦⑿^(qū)間[ti,ti+1]上的直線段作為u(t)在t=ti處切線的近似，即用直線段的斜率作為切線斜率的近似。10第8章第1節(jié)Euler方法從(t0,u0)開(kāi)始，依次應(yīng)用Euler格式，可得到(t1,u1)、(t2,u2)，……。這類(lèi)方法叫單步法。算法實(shí)現(xiàn)過(guò)程：步驟1：輸入必要的初始數(shù)據(jù)t0,u0,T,h步驟2：對(duì)t0≤T循環(huán)求的u1值

t1=t0+h,

u1=u0+h*f(t0,u0),輸出t1,u1,

t0=t1,u0=u1,循環(huán)結(jié)束。11第8章第1節(jié)Euler方法例1：在區(qū)間[0,1.6]上，用Euler方法求解初值問(wèn)題解：t0=0,u0=1,T=1.6,可分別取h=0.1,0.01,0.001進(jìn)行計(jì)算，將計(jì)算所得的數(shù)據(jù)點(diǎn)(t0,u0),(t1,u1),……,連成曲線，并與曲線進(jìn)行比較。結(jié)果表明，h越小，計(jì)算結(jié)果與精確解越接近，但總的來(lái)說(shuō)，Euler格式的精度不高。Euler格式是穩(wěn)定的算法。12第8章第2節(jié)預(yù)估-校正法1.改進(jìn)的Euler方法對(duì)于初值問(wèn)題改進(jìn)的Euler格式(又稱(chēng)為梯形公式)為式(2)只給出了ui和ui+1之間的關(guān)系，在已知的條件下，由需要解方程(2)，式(2)稱(chēng)隱式格式。非線性方程(2)采用迭代法進(jìn)行求解，迭代初值由Euler格式求出。13第8章第2節(jié)預(yù)估-校正法取迭代初值為采用迭代公式求出迭代序列：若f(t,u)在所討論的區(qū)域上關(guān)于變?cè)猽滿足Lipschitz條件(Lipschitz常數(shù)為L(zhǎng))，且步長(zhǎng)h滿足式hL<2則迭代序列收斂于方程(2)的解ui+1。14第8章第2節(jié)預(yù)估-校正法2.預(yù)估-校正法Euler方法的局部截?cái)嗾`差的階為O(h2)，改進(jìn)的Euler方法的局部截?cái)嗾`差的階為O(h3)，即改進(jìn)的Euler方法比Euler方法精確。但改進(jìn)的Euler方法的每一步都要解方程，計(jì)算量較大。如果改進(jìn)的Euler方法只使用迭代公式一次，就作為方程（2）解的近似值，計(jì)算量減小。

(*)和(**)式稱(chēng)為預(yù)估-校正格式。15第8章第2節(jié)預(yù)估-校正法預(yù)估-校正格式可寫(xiě)成如下形式：預(yù)估-校正格式的局部截?cái)嗾`差為O(h3)，與改進(jìn)的Euler方法的局部截?cái)嗾`差階數(shù)相同，但預(yù)估-校正格式的計(jì)算量要小。16第8章第2節(jié)預(yù)估-校正法預(yù)估-校正格式的實(shí)現(xiàn)過(guò)程：步驟1：輸入必要的初始數(shù)據(jù)：t0,u0,T,h步驟2：對(duì)t0≤T循環(huán)求u1值

t1

=t0+h

k1

=f(t0,u0)

k2

=f(t1,u0+h*k1)

u1

=u0+h(k1+k2

)/2

輸出t1,u1t0=t1,u0=u1循環(huán)結(jié)束。17第8章第3節(jié)Runge-Kutta法1.二階Runge-Kutta法對(duì)于初值問(wèn)題采用如下形式的數(shù)值計(jì)算公式如果取c1=c2=1/2,a2=1，得預(yù)估-校正法計(jì)算公式。如果取c1=0，c2=1,a2=1/2，遞推公式變?yōu)?8第8章第3節(jié)Runge-Kutta法上式稱(chēng)為二階Runge-Kutta格式。2.三階Runge-Kutta法求初值問(wèn)題(*)的數(shù)值解采用如下形式計(jì)算公式19第8章第3節(jié)Runge-Kutta法取c1=c3=1/6,c2=2/3,a2=1/2,a3=1,b31=-1,b32=2,得三階Runge-Kutta格式：20第8章第3節(jié)Runge-Kutta法3.四階Runge-Kutta法標(biāo)準(zhǔn)四階Runge-Kutta格式：

21第8章第3節(jié)Runge-Kutta法四階Runge-Kutta格式的實(shí)現(xiàn)過(guò)程：步驟1：輸入必要的初始數(shù)據(jù)t0,u0,T,h步驟2：對(duì)t0≤T循環(huán)求u1值

t1=t0+h,k1=f(t0,u0)k2=f(t0+h/2,u0+h*k1/2)

k3=f(t0+h/2,u0+h*k2/2)

k4=f(t0+h,u0+h*k3)

u1=u0+h

(k1+2k2+2k3+k4)/6

輸出t1,u1

t0=t1,u0=u1,

循環(huán)結(jié)束。22第8章第3節(jié)Runge-Kutta法例2：在區(qū)間[0,1.2]上分別用Euler法、預(yù)估-校正法和四階R-K求下列初值問(wèn)題：解：真實(shí)解為23第8章第4節(jié)線性多步法前面介紹的初值問(wèn)題的數(shù)值算法，都是在求ui+1時(shí)，僅僅用到前一步求出的ui的值，這類(lèi)方法統(tǒng)稱(chēng)為單步法。為了導(dǎo)出精度比較高的數(shù)值方法，在求ui+1時(shí)用到前面若干步的值ui，ui-1，…，此類(lèi)方法統(tǒng)稱(chēng)為多步法。假定下面的初值問(wèn)題在區(qū)間[0,T]上有充分光滑的解u(t);f(t,u)在區(qū)域上關(guān)于足夠多次連續(xù)可微；步長(zhǎng)h取定值，ti=t0+ih；記24第8章第4節(jié)線性多步法1.線性二步法借助于冪級(jí)數(shù)展開(kāi)式并使用待定系數(shù)法，可導(dǎo)出一些多步算法。求ui+2時(shí)，用到ui+1，ui的值，這類(lèi)方法稱(chēng)為二步法。中點(diǎn)公式：中點(diǎn)公式的局部截?cái)嗾`差的階為O(h3)。Milne(米爾恩)格式：Milne格式是隱式格式。25第8章第4節(jié)線性多步法2.Adams外推法設(shè)已知求u(t)在ti+1=t0+(i+1)h處的近似值。將初值問(wèn)題中的微分方程寫(xiě)成如下等價(jià)形式在(1)式已知的情況下，可以求出通過(guò)k+1個(gè)數(shù)據(jù)點(diǎn)做插值多項(xiàng)式得26第8章第4節(jié)線性多步法代入下式得社區(qū)余項(xiàng)得u(ti+1)的近似值ui+1的表達(dá)式插值點(diǎn)在ti和ti+1之間，設(shè)t=ti+th，用Newton后插公式得27第8章第4節(jié)線性多步法部分bm的值如下：m012345bm1?5/123/8251/72095/288上述方法，將插值點(diǎn)取在包含k+1個(gè)節(jié)點(diǎn)的最小區(qū)間[ti-k,ti]之外，故稱(chēng)之為外推法，上面的(*)式稱(chēng)為Adams外推公式。當(dāng)k=0時(shí)，即為Euler公式。28第8章第4節(jié)線性多步法Adams外推公式的局部截?cái)嗾`差的階為O(hk+2)。將高階差商Δmfi-m化為fi-m，…，fi的線性組合得29第8章第4節(jié)線性多步法3.Adams內(nèi)插法通過(guò)k+1個(gè)數(shù)據(jù)點(diǎn)做插值多項(xiàng)式得舍去余項(xiàng)，得u(ti+1)的近似值的表達(dá)式30第8章第4節(jié)線性多步法部分cm的值如下：m0123