常微分方程數(shù)值解法

本文格式為Word版，下載可任意編輯——常微分方程數(shù)值解法第七章常微分方程數(shù)值解法

常微分方程中只有一些典型方程能求出初等解(用初等函數(shù)表示的解)，大部分的方程是求不出初等解的。另外，有些初值問(wèn)題雖然有初等解，但由于形式太繁雜不便于應(yīng)用。因此，有必要探討常微分方程初值問(wèn)題的數(shù)值解法。本章主要介紹一階常微分方程初值問(wèn)題的歐拉法、龍格-庫(kù)塔法、阿達(dá)姆斯方法，在此基礎(chǔ)上推出一階微分方程組與高階方程初值問(wèn)題的數(shù)值解法；此外，還將簡(jiǎn)要介紹求解二階常微分方程值問(wèn)題的差分方法、試射法。

第一節(jié)歐拉法

求解常微分方程初值問(wèn)題

?dy??f(x,y)?dx??y(x0)?y0（1）

的數(shù)值解，就是尋求確鑿解y(x)在一系列離散節(jié)點(diǎn)x0?x1?x2???xn??上的近似值y0,y1,y2,?,yn,?

?yn?稱(chēng)為問(wèn)題的數(shù)值解，數(shù)值解所滿足的離散方程統(tǒng)稱(chēng)為差分格式，hi稱(chēng)為步長(zhǎng)，實(shí)用中常取定步長(zhǎng)。

?xi?xi?1顯然，只有當(dāng)時(shí)值問(wèn)題(1)的解存在且唯一時(shí)，使用數(shù)值解法才有意義，這一前提條件由下面定理保證。

定理設(shè)函數(shù)f?x,y?在區(qū)域D:a?x?b,???y???

上連續(xù)，且在區(qū)域D內(nèi)滿足李普希茲(Lipschitz)條件，即存在正數(shù)L，使得對(duì)于R內(nèi)任意兩點(diǎn)?x,y1?與?x,y2?,恒有的解y?x?存在并且唯一。一、歐拉法(歐拉折線法)

若將函數(shù)

yxf?x,y1??f?x,y2??Ly1?y2則初值問(wèn)題(1)

?在點(diǎn)xhn?處的導(dǎo)數(shù)y?xn?用兩點(diǎn)式代替，即

y??xn??y(xn?1)?y(xn)，再用yn近似地代替y?xn?,則初值問(wèn)題(1)變?yōu)?/p>

?yn?1?yn?hf(xn?yn)(2)?y?y(x),n?0,1,2,?0?0

1

(2)式就是著名的歐拉(Euler)公式。以上方法稱(chēng)為歐拉法或歐拉折線法。歐拉公式有明顯的幾何意義。從幾何上看，求解初值問(wèn)題(1)就是xy平面上求一條通過(guò)點(diǎn)?x0,y0?的曲線y?y?x?，并使曲線上任意一點(diǎn)?x,y?處的切線斜率為

f?x,y?。歐拉公式的幾何意義就是從點(diǎn)

P0x0,y0?出發(fā)作一斜率為f?x0,y0?的

直線交直線x?x1于點(diǎn)P1?x1,y1?，P1點(diǎn)的縱坐標(biāo)y1就是y?x1?的近似值；再?gòu)狞c(diǎn)P1作一斜率為f?x1,y1?的直線x?x2交直線于點(diǎn)P2?x2,y2?,P2點(diǎn)的縱坐標(biāo)y2就是的近似值y?x2?；如此繼續(xù)進(jìn)行，得一條折線P0P1P2?。該折線就是解y?y?x?的近似圖形，如圖7-1。

圖7-1

歐拉法的其它幾種解釋?zhuān)?/p>

1．

假設(shè)y?x?在xn附近展開(kāi)成泰勒級(jí)數(shù)

y?xn?1??y?xn??hy??xn???y?xn??hf?xn,y?xn???hh222y???xn???2y???xn???

取h的線性部分，并用yn作為y?xn?的近似值，得yn?1?yn?hf?xn,yn?

dy2．對(duì)方程

dx?f?x,y?兩邊從xn到xn?1積分，得

y(xn?1)?y?xn???xn?1xnf(x,y(x))dx(3)

2

用矩形公式計(jì)算上式右側(cè)積分，即?xn?1xnf(x,y(x))dx??xn?1xnf?x,y?x??dx

并用yn作為的近似值y?xn?，得yn?1?yn?hf?xn,yn?故歐拉法也稱(chēng)為矩形法。二、改進(jìn)的歐拉法(梯形法)

歐拉法形式簡(jiǎn)單，低，特別當(dāng)曲線y=y(x)計(jì)算便利，但精度比較的曲率較大時(shí)，歐拉法的效果更差。為了達(dá)到較高精度的計(jì)算公式，對(duì)歐拉法進(jìn)行改進(jìn)，用梯形公式計(jì)算(3)式右側(cè)積分，即

?xn?1xnf(x,y(x))dx?h2?f?xn,y?xn???f?xn?1,y?xn?1???

并用yn作為y(xn)的近似值，得到改進(jìn)的歐拉公式

yn?1?yn?h2?f?xn,yn??f?xn?1,yn?1??(4)

上述方法稱(chēng)為改進(jìn)的歐拉法，也稱(chēng)為梯形法。

不難發(fā)現(xiàn)，歐拉公式是關(guān)于yn?1的顯式，即只要已知yn,經(jīng)過(guò)一次計(jì)算便可得yn?1的值，而改進(jìn)的歐拉公式是以yn?1的隱式方程給出，不能直接得到y(tǒng)n?1。隱式方程(4)尋常用迭代法求解，而迭代過(guò)程的實(shí)質(zhì)是逐步顯式化。

先用歐拉公式

yn?1?yn?hf?xn,yn??0?

給出yn?1的迭代初值，然后再用改進(jìn)的歐拉公式(4)進(jìn)行迭代，即有

(0)?yn?1?yn?hf(xn,yn)?h?(k?1)(k)y?y?f(xn,yn)?f(xn?1,yn?1)?n?1n(5)2?k?0,1,2,?????迭代過(guò)程進(jìn)行到連續(xù)兩次迭代結(jié)果之差的絕對(duì)值小于給定的精度?即

yn?1?yn?1???k?1?k

為止，這時(shí)取

yn?1?yn?1?k?1?

然后再轉(zhuǎn)入下一步計(jì)算。

3

???y?下面探討是否收斂；若收斂，它的極限是否滿足(4)式。

kn?1假設(shè)f(x,y)滿足李普希茲條件f?x,y1??f?x,y2??L?y1?y2?則

yn?1?yn?1???k?1??k?h2f?xn?1,yn?1??f?xn?1,yn?1?k??k?1??hL2yn?1?yn?12?k??k?1??hL?????2?yn?1?yn?12?k?1?k?2?hL??????2??yn?1?yn?1?1??0?

?hL??k→∞時(shí)，有??2?k由此可見(jiàn)，只要

kn?1hL2?1(這里只要步長(zhǎng)h足夠小即可)，當(dāng)

?0，

???y?所以收斂。又由于f(x,y)對(duì)y連續(xù)，當(dāng)k→∞時(shí)，對(duì)等式

yn?1?k?1??yn?h2?f?xn,yn??fxn?1,yn?1??k???

兩端取極限，得

yn?1?yn?h2?f?xn,yn??f?xn?1,yn?1??kn?1

n?1???y?因此，只要步長(zhǎng)h足夠小，就可保證收斂到滿足(4)式的y。

三、預(yù)估一校正法

改進(jìn)的歐拉公式在實(shí)際計(jì)算時(shí)要進(jìn)行屢屢迭代，因而計(jì)算量較大。在實(shí)用上，對(duì)于改進(jìn)的歐拉公式(5)只迭代一次，即先用歐拉公式算出yn?1的預(yù)估值y

(0)n?1，再用改進(jìn)的歐拉公式(4)進(jìn)行一次迭代得到校正值yn?1，即

(0)?yn?1?yn?hf(xn,yn)?h?yn?1?yn??f?xn,yn??f?xn?1,yn?1??,n?0,1,2,??2?(6)

預(yù)估—校正公式也常寫(xiě)成以下形式：

4

11?y?y?k?k2n1?n?122?k1?hf(xn,yn),n?0,1,2,???k?hf?x?h,y?k?nn1?2?(7)

四、公式的截?cái)嗾`差

定義若某種微分方程數(shù)值解公式的截?cái)嗾`差是O(h

k?1)，則稱(chēng)這種方法是

k階方法。為了簡(jiǎn)化分析，在進(jìn)行誤差分析時(shí)，我們假設(shè)前一步的結(jié)果是確鑿的，即在yn=y(xn)的前提下，考慮用yn?1作為y(xn?1)的近似值而產(chǎn)生的截?cái)嗾`差，這種誤差稱(chēng)為局部截?cái)嗾`差。由泰勒公式

y?xn?1??y?xn?h??y?xn??by??xn??h22!y???xn???

對(duì)于歐拉公式，有

yn?1?yn?hf?xn,yn??y?xn??hy??xn?于是

y?xn?1??yn?1?Oh2??

2則歐拉公式的截?cái)嗾`差為O(h),所以歐拉法是一階方法。

對(duì)于預(yù)估—校正公式，有

k1?hf?xn,yn??hy??xn?k2?hf?xn?h,yn?k1??hf?xn?h,y?xn??k1??hf?xn,y?xn???hf?hf?xn,y?xn???h2?x?xn,y?xn???k1fy?xn,y?xn?????xn?f?x,y?xn???y??xn?fy?xn,y?xn?????

而

y??x??f?x,y?x??y???x??fx?x,y?x