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

《常微分方程數(shù)值解法》由會員分享，可在線閱讀，更多相關(guān)《常微分方程數(shù)值解法（37頁珍藏版）》請在裝配圖網(wǎng)上搜索。

1、數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS第8章常微分方程 實(shí)際中，很多問題的數(shù)學(xué)模型都是微分方程。我們可以研究它們的一些性質(zhì)。但是，只有極少數(shù)特殊的方程有解析解。對于絕大部分的微分方程是沒有解析解的。常微分方程作為微分方程的基本類型之一，在自然界與工程界有很廣泛的應(yīng)用。很多問題的數(shù)學(xué)表述都可以歸結(jié)為常微分方程的定解問題。很多偏微分方程問題，也可以化為常微分方程問題來近似求解。本章討論常微分方程的數(shù)值解法數(shù) 學(xué) 系University of Science and Technology

2、of ChinaDEPARTMENT OF MATHEMATICS對于一個常微分方程：,),(baxyxfdxdyy通常會有無窮個解。如：Raaxyxdxdy,)sin()cos(因此，我們要加入一個限定條件。通常會在端點(diǎn)出給出，如下面的初值問題：0)(,),(yaybaxyxfdxdy為了使解存在唯一，一般，要加限制條件在f上，要求f對y滿足Lipschitz條件：2121),(),(yyLyxfyxf數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS 常微分方程的解是一個函數(shù)，但是，計算機(jī)沒

3、有辦法對函數(shù)進(jìn)行運(yùn)算。因此，常微分方程的數(shù)值解并不是求函數(shù)的近似，而是求解函數(shù)在某些節(jié)點(diǎn)的近似值。例：我們對區(qū)間做等距分割：,()/ixh ihbam設(shè)解函數(shù)在節(jié)點(diǎn)的近似為iy由數(shù)值微分公式，我們有iixxxxyxfdxdy),(，則：),(1iiiiyxfhyy向前差商公式),(1iiiiyxfhyy可以看到，給出初值，就可以用上式求出所有的iy數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS基本步驟如下：解差分方程，求出格點(diǎn)函數(shù) 對區(qū)間作分割：bxxxanI10:求 在 上的近似值 。)(

4、xyixiyiy稱為分割I(lǐng)上的格點(diǎn)函數(shù) 由微分方程出發(fā)，建立求格點(diǎn)函數(shù)的差分方程。這個方程應(yīng)該滿足：A、解存在唯一；B、穩(wěn)定，收斂；C、相容數(shù)值方法，主要研究步驟，即如何建立差分方程，并研究差分方程的性質(zhì)。這種方法，稱為數(shù)值離散方法。求的是在一系列離散點(diǎn)列上，求未知函數(shù)y在這些點(diǎn)上的值的近似。我們的目的，就是求這個格點(diǎn)函數(shù)數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS為了考察數(shù)值方法提供的數(shù)值解，是否有實(shí)用價值，需要知道如下幾個結(jié)論：步長充分小時，所得到的數(shù)值解能否逼近問題得真解；即收斂性問

5、題 誤差估計 產(chǎn)生得舍入誤差，在以后得各步計算中，是否會無限制擴(kuò)大；穩(wěn)定性問題數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS8.1 Euler公式公式做等距分割，利用數(shù)值微分代替導(dǎo)數(shù)項，建立差分方程。imabxiI:1、向前差商公式)(2)()()(1nnnnyhxyhxyxy)(2)(,()()(1nnnnnyhxyxfhxyxy)(2)(,()()(21nnnnnyhxyxhfxyxy所以，可以構(gòu)造差分方程),(1nnnnyxhfyy稱為局部截斷誤差。顯然，這個誤差在逐步計算過程中會傳播

6、，積累。因此還要估計這種積累數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS定義定義在假設(shè)在假設(shè) yi=y(xi)，即第，即第 i 步計算是精確的前提下，考慮步計算是精確的前提下，考慮的截斷誤差的截斷誤差 Ri=y(xi+1)yi+1 稱為稱為局部截斷誤差局部截斷誤差/*local truncation error*/。定義定義若某算法的局部截斷誤差為若某算法的局部截斷誤差為O(hp+1)，則稱該算法有，則稱該算法有p 階精度。階精度。記為2、收斂性)(2)(,()()(21nnnnnyhxy

7、xhfxyxy),(1nnnnyxhfyy1nhT考察局部誤差的傳播和積累數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS111)(nnnyxye1),()(,()(nnnnnnnThyxfxyxfhyxy1)(nnnnThyxyhLejjnTThTehLmax ,)1(hThTehLhLn1)1()1(hThLehLn1)1()1(12hThLhTehLhLn1)1()1()1(22hThLhLehLn1)1()1()1(223數(shù) 學(xué) 系University of Science and T

8、echnology of ChinaDEPARTMENT OF MATHEMATICShThLhLehLnn1)1()1()1(01hThLhLehLnn)1(1)1(1)1(101hThLhLehLnn101)1()1(LTehLn01)1(LTehLn)1()(1hOen)()1(00hOTexenxn稱為整體截斷誤差是1階方法數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS3、穩(wěn)定性誤差在以后各步的計算中不會無限制擴(kuò)大。是格式對舍入誤差的抑止作用我們考慮一種簡單情況，即僅初值有誤差，而其

9、他計算步驟無誤差。設(shè)iz是初值有誤差后的計算值，則),(),(11nnnnnnnnzxhfzzyxhfyy所以，我們有：),(),(111nnnnnnnnzxfyxfhezye)1(hLezyhLennnnhLnneehLe)1(010)1(可以看出，向前差商公式關(guān)于初值是穩(wěn)定的。當(dāng)初始誤差充分小，以后各步的誤差也充分小數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS4、向后差商公式)(2)()()(11nnnnyhxyhxyxy)(2)(,()()(111nnnnnyhxyxfhxyxy)(

10、2)(,()()(2111nnnnnyhxyxhfxyxy),(111nnnnyxhfyy是隱格式，要迭代求解)0(1)(11)1(1),(nknnnknyyxhfyy可以由向前差商公式求出數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS5、中心差商公式)(2)()()(111nnnnyhxyhxyxy),(1111nnnnyxhfyy是多步，2階格式，該格式不穩(wěn)定6、梯形法基于數(shù)值積分的公式對微分方程,),(baxyxfdxdyy做積分，則：1),(nnxxyxfdxdy數(shù) 學(xué) 系Unive

11、rsity of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS)(21hOen類似，可以算出其誤差估計式：2階的方法所以，有格式為：),(),(2111nnnnnnyxfyxfhyy是個隱式的方法，要用迭代法求解11()()(,()nnxnnxy xy xf x y x dx)(12)(,()(,(2)()(2111fhxyxfxyxfhxyxynnnnnn局部截斷誤差數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS8.2 R

12、ungeKutta法法由Taylor展開)()!1()(!)()()()1(1)(1nkknkknnnykhxykhxhyxyxy記為1nhT),()(yxfxy),(),()(yyxfyxfxyyx)(xy所以，可以構(gòu)造格式),(),(),(!2),(21nnnnynnxnnnnyxfyxfyxfhyxhfyy這種格式使用到了各階偏導(dǎo)數(shù)，使用不便。數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS從另一個角度看，11),(,()()(nnnnnhTfxyxhhFxyxy取(x,y)及其附近的點(diǎn)

13、做線性組合，表示F，問題就好辦了。當(dāng)然，要求此時的展開精度相同。這種方法稱為RungeKutta法數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS在(x,y)處展開，)(),(),(),(),(),(),(221221hOyxfyxhfbyxhfayxfcyxfcfyxhFyx比較11),(),(),(!2),()()(nnnnnynnxnnnnhTyxfyxfyxfhyxfhxyxy以2階為例，設(shè)),(,(),(),(21221yxhfbyhaxfcyxfcfyxhF數(shù) 學(xué) 系Univers

14、ity of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS有：2/12/112122221bcaccc112/121221bacc1、改進(jìn)的Euler公式),(),(2/121211hKyhxfKyxfKKKhyynnnnnn3/23/24/3,4/121221bacc2、Heun公式數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS一般的RungeKutta法構(gòu)造111212212211,(),(),(),(miimimmmm

15、KbhyhaxfKhKbyhaxfKyxfKKcKcKcfyxhF常見的為3階，4階公式數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS8.3 線性多步法線性多步法用用若干若干節(jié)點(diǎn)處的節(jié)點(diǎn)處的 y 及及 y 值的值的線性組合線性組合來近似來近似y(xn+1)。).(.110111101knknnnknknnnffffhyyyy b bb bb bb ba aa aa a其通式可寫為：其通式可寫為：),(nnnyxff 當(dāng)當(dāng) b b 1 0 時，為時，為隱式公隱式公式式;b b 1=0 則為則為

16、顯式公式顯式公式。數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS 基于數(shù)值積分的構(gòu)造法基于數(shù)值積分的構(gòu)造法將將 在在 上積分，得到上積分，得到),(yxfy,1npnxx1)(,()()(1ipnxxpnndxxyxfxyxy只要只要近似地算出右邊的積分近似地算出右邊的積分 ，則可通，則可通過過 近似近似y(xn+1)。而。而選用不同近似式選用不同近似式 Ik，可得到不，可得到不同的計算公式同的計算公式。1)(,(npnxxkdxxyxfIkpnnIyy1數(shù) 學(xué) 系University of

17、 Science and Technology of ChinaDEPARTMENT OF MATHEMATICS1)()!1()()2(npnxxqqdxxqy1)(npnxxdxxy若積分用節(jié)點(diǎn)qnnnxxx,1作為積分點(diǎn)，則有1110)()()()(1nqnqnnxxhTxyaxyaxyahdxxynpn)(,()(nnnxyxfxy101)(,()()(nqjjnjnjpnnhTxyxfahxyxy積分系數(shù)1)(npnxxjjdxxlha這是顯格式，q+1階r+1步格式。r=maxp,q11,qnnnxxx為積分節(jié)點(diǎn)，可以構(gòu)造r+1步q+1階隱格式局部截斷誤差同樣，若以數(shù) 學(xué) 系Uni

18、versity of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS例：建立p=1,q=2的顯格式p=1，q=2，顯格式，11)(nnxxdxxy積分區(qū)間為積分節(jié)點(diǎn)為21,nnnxxxhdxxxxxxxxxhannxxnnnnnn37)()(1121210所以hdxxxxxxxxxhannxxnnnnnn32)()(1121121hdxxxxxxxxxhannxxnnnnnn31)()(1112212數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF M

19、ATHEMATICS)(31)()()!3()()4(421)4(111yhdxxxxxxxyTnnxxnnnn例：建立p=2,q=2的隱格式p=2，q=2，隱格式，12)(nnxxdxxy積分區(qū)間為積分節(jié)點(diǎn)為11,nnnxxxhdxxxxxxxxxhannxxnnnnnn43)()(1211110所以0)()(1211111nnxxnnnnnndxxxxxxxxxha數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS)(83)()()!3()()4(411)4(111yhdxxxxxxxyTn

20、nxxnnnnhdxxxxxxxxxhannxxnnnnnn49)()(1211110它的截斷誤差較 顯格式 小，通常也具有更好的穩(wěn)定性。Adams公式公式 p=0 時候的多步法時候的多步法參見書數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS8.4 方程組和高階方程的數(shù)值解法方程組和高階方程的數(shù)值解法bxaayayyyxfdxdyyyxfdxdymmmmmm ,)()(),(),(111111寫成向量的形式：)(),(aYYxFdxdY數(shù) 學(xué) 系University of Science a

21、nd Technology of ChinaDEPARTMENT OF MATHEMATICS各種方法都可以直接運(yùn)用過來。bxazazyayzyxgdxdzzyxfdxdy ,)()(),(),(00Euler公式),(),(11nnnnnnnnnnzyxhgzzzyxhfyy以兩個方程的方程組為例數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICSRunge-Kutta公式112341(22)6nnnnyyhKKKKzz)2,2(12KhYhxFK1(1)(2)112(1)(2)11(,)(,)

22、(,)222(,)222nnnnnnnnnnnnf xyzKg xyzhhhf xyKzKKhhhg xyKzK 數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS(1)(1)223(1)(1)22(1)(1)334(1)(1)33(,)222(,)222(,)(,)nnnnnnnnnnnnhhhf xyKzKKhhhg xyKzKf xh yhKzhKKg xh yhKzhK數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT

23、 OF MATHEMATICS0.05(1)0.002200.09(1)0.1515(0)0.193(0)0.083duuuuvdtdvvvuvdtuv1、(,)0.05(1/20)0.002(,)0.09(1/15)0.15f u v tuuuvFg u v tvvuv2、確定方法，然后求解（0.20276 0.0881157）（0.213007 0.0934037）（0.223763 0.0988499）（0.235052 0.104437）（0.246902 0.110146）4階Runge-Kutta法，h=1數(shù) 學(xué) 系University of Science and Technol

24、ogy of ChinaDEPARTMENT OF MATHEMATICS高階方程bxaayayayyyyxfdxydmmmmm ,)()()(),()1(21)1(1數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS則有：bxaayayyyxfdxdyydxdymmmm ,)()(),(11121令mmydxdyydxdyyy1211數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS例：例：考察

25、初值問題考察初值問題 在區(qū)間在區(qū)間0,0.5上的解。上的解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計算數(shù)值解。分別用歐拉顯、隱式格式和改進(jìn)的歐拉格式計算數(shù)值解。1)0()(30)(yxyxy0.00.10.20.30.40.5精確解精確解改進(jìn)歐拉法改進(jìn)歐拉法 歐歐拉隱式拉隱式歐拉顯式歐拉顯式 節(jié)點(diǎn)節(jié)點(diǎn) xixey30 1.0000 2.0000 4.0000 8.0000 1.6000 101 3.2000 101 1.00002.5000 10 1 6.2500 10 21.5625 10 23.9063 10 39.7656 10 41.00002.50006.25001.5626 101

26、3.9063 1019.7656 1011.00004.9787 10 22.4788 10 31.2341 10 46.1442 10 63.0590 10 7What is wrong?!數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS8.5 差分方程的絕對穩(wěn)定性差分方程的絕對穩(wěn)定性對于一般的差分方程kjjnjnjkjjnjyxfbhya00),(由初始誤差產(chǎn)生了差分解的誤差，實(shí)際上是同一差分方程，取不同初值所得到的2組差分解之間的差。這個差不僅于差分方程本身有關(guān)，而且與微分方程本身有關(guān)。

27、如果微分方程本身是不穩(wěn)定，那就沒理由要求這2組解充分接近。因此，差分方程的穩(wěn)定性概念是建立在微分方程穩(wěn)定的基礎(chǔ)上的。把這個典型微分方程規(guī)定為：仍然考慮最簡單的模型，即只有初值產(chǎn)生誤差，看看這個誤差的傳播。0)()0(Re yayydxdy數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS差分方程運(yùn)用到如上的微分方程后，可以得到0)Re(,00kjjnjkjjnjybhya對于給定的初始誤差110,keee，誤差方程具有一樣的形式0)Re(,00kjjnjkjjnjebhea數(shù) 學(xué) 系Univer

28、sity of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS定義：差分方程稱為絕對穩(wěn)定的，若差分方程作用到微分方程)0(Re ydxdy時，對任意的初值，總存在左半復(fù)平面上的一個區(qū)域，當(dāng) 在這個區(qū)域時，差分方程的解趨于0。這個區(qū)域稱為穩(wěn)定區(qū)域h例：向后Euler公式的穩(wěn)定性11nnnhyyy誤差方程：11nnnheeeheenn1111210ReImg數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS考察隱式歐拉法考察隱式歐拉法

29、11 iiiyhyy 111iiyyh11011iih可見絕對穩(wěn)定區(qū)域?yàn)椋嚎梢娊^對穩(wěn)定區(qū)域?yàn)椋簗1|1h210ReImg注：注：一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式一般來說，隱式歐拉法的絕對穩(wěn)定性比同階的顯式法的好。法的好。數(shù) 學(xué) 系University of Science and Technology of ChinaDEPARTMENT OF MATHEMATICS3階RungeKutta211132111)211(46hhyKhyKyKKKKhyynnnnn32161211hhhyynn顯式顯式 1 4 階方法的絕對穩(wěn)定區(qū)域?yàn)殡A方法的絕對穩(wěn)定區(qū)域?yàn)閗=1k=2k=3k=4-1-2-3-123ReImg