《計(jì)算機(jī)數(shù)值方法教學(xué)課件》第二章-常微分方程數(shù)值解法

1、第二章 常微分方程數(shù)值解法Chapter 2 Numerical Solution of Ordinary Differential Equation(s)第二章 常微分方程數(shù)值解法Chapter 2 Numer22.1 引 言 Fm假設(shè)牽引力F為恒定值為了確定待定常數(shù)，可以給定初始條件：或者給定邊界條件：假設(shè)牽引力不恒定呢？求速度22.1 引 言 Fm假設(shè)牽引力F為恒定值為了確定待定常3雖然求解微分方程有許多解析方法，但解析方法只能夠求解一些特殊類型的方程。還有一類近似方法稱為數(shù)值方法，它可以給出解在一些離散點(diǎn)上的近似值。利用計(jì)算機(jī)解微分方程主要使用數(shù)值方法。 解析方法與數(shù)值方法3雖然求解微

2、分方程有許多解析方法，但解析方法只能夠求解一些特4主要研究對(duì)象：初值問題求數(shù)值解：求y(x)在離散數(shù)據(jù)點(diǎn)xk處的近似值yk 。 y=y(x)xyx0=ax1x2x3xk-1xn-1xn=b xk4主要研究對(duì)象：初值問題求數(shù)值解：求y(x)在離散數(shù)據(jù)點(diǎn)xk5則稱 f (x,y) 對(duì)y 滿足李普希茲條件，此時(shí)初值問題在a,b上存在唯一的連續(xù)可微的解。定理1：設(shè) f (x,y) 是定義在區(qū)域 G=(x,y)|axb, yR上的連續(xù)函數(shù)，若存在正的常數(shù) L 使：(Lipschitz)條件有解條件：5則稱 f (x,y) 對(duì)y 滿足李普希茲條件，此時(shí)初值問題6則Lipschitz條件成立：在 f (x,

3、y) 對(duì)y可微的情況下，若偏導(dǎo)數(shù)有界：有解條件的判斷：6則Lipschitz條件成立：在 f (x,y) 對(duì)y可微7定理2：如果f(x,y)在G=(x,y)|axb, yR上滿足Lipschitz條件，則初值問題是適定的。適定性：指初值問題中，初始值y0及微分方程的右端函數(shù)f (x,y) 有微小變化時(shí)，只能引起解的微小變化。適定性條件：7定理2：如果f(x,y)在G=(x,y)|axb, 8在解的存在區(qū)間 a, b上取n + 1個(gè)節(jié)點(diǎn) 這里把 稱為由xi到xi+1的步長 一般取成等間距的： 求解方法：步進(jìn)法（分為單步法和多步法） 數(shù)值方法的基本思想8在解的存在區(qū)間 a, b上取n + 1個(gè)節(jié)點(diǎn)

4、 這里把 9本章規(guī)定：在 處初值問題的理論解用 表示，數(shù)值解法的近似解用 表示。記 ，它和 是不同的，后者等于 。9本章規(guī)定：在 處初值問題的理論解用 102.2 幾種簡單的數(shù)值方法 （一） 歐拉（Euler）法102.2 幾種簡單的數(shù)值方法 （一） 歐拉（Euler）11、泰勒公式解釋 、求導(dǎo)的兩點(diǎn)公式解釋、積分公式解釋歐拉公式的的分析解釋11、泰勒公式解釋 歐拉公式的的分析解釋12泰勒公式解釋其中：可以得到：12泰勒公式解釋其中：可以得到：13求導(dǎo)的兩點(diǎn)公式解釋可以得到：13求導(dǎo)的兩點(diǎn)公式解釋可以得到：14對(duì)微分方程(1.1)兩端從進(jìn)行積分積分公式解釋14對(duì)微分方程(1.1)兩端從進(jìn)行積分

5、積分公式解釋15右端積分用左矩形數(shù)值求積公式：即：15右端積分用左矩形數(shù)值求積公式：即：16歐拉公式的的幾何描述yxx0 x1 x2 x3 x4y=y(x)16歐拉公式的的幾何描述yxx0 x1 x17例題1：（取步長h=0.1）用Euler方法求滿足條件 的y(t) 數(shù)值解。 解：17例題1：（取步長h=0.1）用Euler方法求滿足條件 18ntnyny(tn)y(tn)-yn01234101.01.11.21.31.43.00.00.271830.684761.276983.0935515.398240.00.345920.866641.607223.6203618.6830.00.07

6、4090.181880.330240.526813.28486數(shù)值解列表為 18ntnyny(tn)y(tn)-yn01.00.00.019歐拉方法的誤差估計(jì) 通過數(shù)值方法進(jìn)行計(jì)算時(shí)，考慮每一步產(chǎn)生的誤差，從x0開始一步步累積到xn，稱 為該數(shù)值方法在xn點(diǎn)處的整體截?cái)嗾`差，該誤差與xn 及之前的各步計(jì)算誤差都有關(guān)系。19歐拉方法的誤差估計(jì) 通過數(shù)值方法進(jìn)行計(jì)算時(shí)，考慮每一步產(chǎn)20歐拉方法的誤差估計(jì) 為了簡化分析，著重分析xn點(diǎn)單步計(jì)算產(chǎn)生的誤差，即把xn點(diǎn)之前的計(jì)算當(dāng)作無誤差： 稱該誤差為數(shù)值方法在xn+1點(diǎn)處的局部截?cái)嗾`差。局部截?cái)嗾`差的第一個(gè)非零項(xiàng)為局部截?cái)嗾`差主項(xiàng)。20歐拉方法的誤差估

7、計(jì) 為了簡化分析，著重分析xn點(diǎn)單步計(jì)算21歐拉方法的誤差估計(jì) 如果求解公式的局部截?cái)嗾`差為R (h)=O(hp+1)，則稱該求解公式具有p階精度，稱該方法為p階方法。定義：歐拉方法：具有1階精度。21歐拉方法的誤差估計(jì) 如果求解公式的局部截?cái)嗾`差為R (h22（二）向后歐拉法 （1）方法 其公式為：22（二）向后歐拉法 （1）方法 其公式為：23（2）局部截?cái)嗾`差23（2）局部截?cái)嗾`差24例題2：用向后Euler法解初值問題 向后Euler法的公式為 解：x0=0, y0=1, 取h=0.124例題2：用向后Euler法解初值問題 向后Euler法的25方法比較及推廣： Euler方法 顯式

8、公式向后Euler方法 隱式公式 解一個(gè)非線性方程 難求解 顯式和隱式相結(jié)合 隱式的顯化 25方法比較及推廣： Euler方法 顯式公式26計(jì)算公式為： 由顯式得到，稱為預(yù)估值；yn+1由隱式得到，稱為校正值。這種求解方法統(tǒng)稱為預(yù)估校正方法。其求解過程為：26計(jì)算公式為： 由顯式得到，稱為預(yù)估值；這種求27例3 用預(yù)估校正方法求解微分方程（取h=0.1）：解： 27例3 用預(yù)估校正方法求解微分方程（取h=0.1）：解28（三）梯形公式 28（三）梯形公式 29梯形公式局部截?cái)嗾`差29梯形公式局部截?cái)嗾`差30預(yù)估-校正方法：稱為改進(jìn)的Euler求解公式或改進(jìn)Euler法。30預(yù)估-校正方法：稱為

9、改進(jìn)的Euler求解公式或改進(jìn)Eul31為了表示方便，可以改寫為：31為了表示方便，可以改寫為：32（四）歐拉方法的收斂性分析由初值問題的單步法產(chǎn)生的近似解 ，如果對(duì)于任一固定的 均有 ，則稱該方法是收斂的。定義：局部截?cái)嗾`差：若初值問題的一個(gè)單步法的局部截?cái)嗾`差為： 定理：且增量函數(shù)關(guān)于y滿足Lipschitz條件，則整體截?cái)嗾`差： 整體截?cái)嗾`差比局部截?cái)嗾`差低1階32（四）歐拉方法的收斂性分析由初值問題的單步法產(chǎn)生的近似解33證明：存在常數(shù)c，使得33證明：存在常數(shù)c，使得34當(dāng) 固定時(shí)，所以34當(dāng) 固定時(shí)，所以35（五）歐拉方法的穩(wěn)定性分析問題：常微分方程初值問題數(shù)值解的每步計(jì)算都是在前

10、一步計(jì)算的結(jié)果上進(jìn)行的，所以必須考慮前面的誤差對(duì)以后計(jì)算結(jié)果的影響，誤差的積累會(huì)不會(huì)蓋過真解呢？選用代表性試驗(yàn)方程： y=y (Re()0，在計(jì)算yn時(shí)引入了誤差n 。若這個(gè)誤差在計(jì)算后面的yn+k(k=1,2,.)中所引的誤差n+k按絕對(duì)值均不增加，就說這個(gè)數(shù)值方法對(duì)于這個(gè)步長h和復(fù)數(shù)是絕對(duì)穩(wěn)定的。 若在區(qū)域R內(nèi)數(shù)值方法是絕對(duì)穩(wěn)定的，則R為該數(shù)值方法的絕對(duì)穩(wěn)定區(qū)域(區(qū)間)。 35（五）歐拉方法的穩(wěn)定性分析問題：常微分方程初值問題數(shù)值解36把歐拉方法用于試驗(yàn)方程： y=y誤差方程 ：要求誤差不增加： O-2-1Re(h)Im(h)36把歐拉方法用于試驗(yàn)方程： y=y誤差方程 ：要求誤差37把向

11、后歐拉方法用于試驗(yàn)方程： y=yO21Re(h)Im(h)要求誤差不增加： 可見隱式的向后Euler方法比顯式的Euler方法的絕對(duì)穩(wěn)定域要大得多。同階精度的數(shù)值方法，往往隱式方法比顯式方法的絕對(duì)穩(wěn)定域大。 37把向后歐拉方法用于試驗(yàn)方程： y=yO21Re(h38例4：以 y=y為例判斷梯形公式的穩(wěn)定性：解出：這種穩(wěn)定性稱為無條件穩(wěn)定！38例4：以 y=y為例判斷梯形公式的穩(wěn)定性：解出：這種39例5：用Euler法、向后Euler法、改進(jìn)的歐拉法（梯形公式）解初值問題 取步長h=0.2，小數(shù)點(diǎn)后至少保留4位。 解：Euler法向后Euler法39例5：用Euler法、向后Euler法、改進(jìn)的

12、歐拉法（梯40改進(jìn)的Euler法梯形公式法40改進(jìn)的Euler法梯形公式法41h=0.23.00003.30083.46593.55653.60623.63351.01.21.41.61.83.0真解改進(jìn)的歐拉法梯形公式向后EulerEulerxi 3.0000 3.0000 3.0000 3.0000 3.4000 3.2500 3.3077 3.2800 3.5600 3.4063 3.4734 3.4424 3.6240 3.5039 3.5626 3.5366 3.6496 3.5649 3.6106 3.5912 3.6598 3.6031 3.6365 3.6229如何得到高精度的

13、求解公式？41h=0.23.00001.0真解改進(jìn)的歐拉法梯形公式向后422.3 Runge-Kutta法 如何構(gòu)造更高精度的求解公式？一種思路是：對(duì)微分方程右端積分采用高次多項(xiàng)式近似，如用二次多項(xiàng)式近似可得到422.3 Runge-Kutta法 如何構(gòu)造更高精度43該公式稱為Simpson公式。該公式的局部截?cái)嗾`差為：43該公式稱為Simpson公式。44Taylor展開：另一種得到高階方法的想法是直接利用泰勒級(jí)數(shù)展開。如果能計(jì)算得到y(tǒng)的高階微商,則可寫出r階的計(jì)算公式44Taylor展開：另一種得到高階方法的想法是直接利用泰勒45例題：取h=0.1，求解初值問題 解： 一階公式：二階公式：

14、四階公式：45例題：取h=0.1，求解初值問題 解： 一階公式：46真解： p0.10.20.30.40.511.100001.221001.370081.557791.8004621.110001.246891.421741.662621.9208741.111101.249661.428481.666451.99942y(xn)1.111111.250021.428571.666672.00000 xn問題： 1、求微商麻煩； 2、計(jì)算量大。46真解： p0.10.20.30.40.511.100047R-K方法是通過對(duì)不同點(diǎn)上的函數(shù)值做線性組合，構(gòu)造近似公式，把近似公式和泰勒展開相比較，

15、使前面的若干項(xiàng)相吻合，從而使近似公式達(dá)到一定的階數(shù)。 問題：如何組合函數(shù)值？其中，47R-K方法是通過對(duì)不同點(diǎn)上的函數(shù)值做線性組合，構(gòu)造近似公48選擇參數(shù)ci, ai, bij 的原則是，要求的Runge-Kutta法的右端項(xiàng)在(xn,yn)處泰勒展開后按h的冪次重新整理得到的結(jié)果與微分方程的解y(xn+1)在xn處的Taylor展開式 有盡可能多的項(xiàng)重合。 48選擇參數(shù)ci, ai, bij 的原則是，要求的Rung49以計(jì)算兩個(gè)函數(shù)值為例說明選擇c1、c2、a2、b21 階最高為了計(jì)算49以計(jì)算兩個(gè)函數(shù)值為例說明選擇c1、c2、a2、b21 階50于是若要求局部截?cái)嗾`差達(dá)到,則要求有 選取

16、50于是若要求局部截?cái)嗾`差達(dá)到,則要求有 選取51二階R-K方法的公式 它的局部階段誤差為O(h3)。這是計(jì)算兩次函數(shù)值的情況下所能達(dá)到的最高階。 51二階R-K方法的公式 它的局部階段誤差為O(h3)。這是52中間點(diǎn)法： 二階休恩(Heun)法： 52中間點(diǎn)法： 二階休恩(Heun)法： 53經(jīng)典的R-K方法是一個(gè)四階的方法,公式為 1、一步法，可以自開始；特點(diǎn)：2、精度較高；3、便于改變步長；4、計(jì)算量較小。53經(jīng)典的R-K方法是一個(gè)四階的方法,公式為 1、一步法，可54例2：用二階R-K方法和四階R-K方法(h=0.1)求解解：二階R-K方法四階R-K方法54例2：用二階R-K方法和四階

17、R-K方法(h=0.1)求解55xn二階R-K誤差四階R-K誤差y=1-e-5x20.10.20.30.40.50.60.70.80.91.0 0.0500 0.0012 0.0488 0.0000 0.0488 0.1830 0.0017 0.1813 0.0000 0.1813 0.3627 0.0004 0.3624 0.0000 0.3624 0.5475 -0.0031 0.5507 -0.0000 0.5507 0.7059 -0.0076 0.7134 -0.0001 0.7135 0.8235 -0.0112 0.8346 -0.0001 0.8347 0.9012 -0.01

18、25 0.9135 -0.0002 0.9137 0.9476 -0.0116 0.9590 -0.0002 0.9592 0.9733 -0.0093 0.9823 -0.0002 0.9826 0.9866 -0.0066 0.9931 -0.0002 0.9933說明：1、高階方法比低階方法精度高，但需要y(x)光滑性好；2、步長h并非越小越好；步長選擇從大到小。55xn二階R-K誤差四階R-K誤差y=1-e-5x20.156數(shù)學(xué)家Butcher于1965年證明了運(yùn)算量與可以達(dá)到的最高精度階數(shù)關(guān)系如下：每步計(jì)算kR的個(gè)數(shù)2345678910可以達(dá)到最高精度階數(shù)23445666R-2可見，

19、通常只選用R 4的公式是合適的。56數(shù)學(xué)家Butcher于1965年證明了運(yùn)算量與可以達(dá)到的57自動(dòng)選步長方法 設(shè)在xn處之前誤差不計(jì)：誤差：57自動(dòng)選步長方法 設(shè)在xn處之前誤差不計(jì)：誤差：58隱式R-K方法類似于顯式R-K公式，稍加改變，就得到隱式R-K方法。顯式公式中對(duì)系數(shù)求和的上限是i-1。而在隱式公式中對(duì)系數(shù)求和的上限是L，需要用迭代法求出Ki。推導(dǎo)隱式公式的思路和方法與顯式R-K法類似。通常，同級(jí)的隱式公式可以獲得比顯式公式更高的階。58隱式R-K方法類似于顯式R-K公式，稍加改變，就得到隱式59通常，同級(jí)的隱式公式獲得比顯式公式更高的階。常用的隱式R-K法有:1級(jí)2階中點(diǎn)公式 :

20、2級(jí)2階梯形公式： 2級(jí)4階R-K公式： 59通常，同級(jí)的隱式公式獲得比顯式公式更高的階。常用的隱式R60隱式龍格-庫塔法其中2階方法 的絕對(duì)穩(wěn)定區(qū)域?yàn)?ReImg而顯式 1 4 階方法的絕對(duì)穩(wěn)定區(qū)域?yàn)閗=1k=2k=3k=4-1-2-3-123ReImg無條件穩(wěn)定龍格-庫塔法穩(wěn)定區(qū)域60隱式龍格-庫塔法其中2階方法 612.4 線性多步法 單步法的優(yōu)缺點(diǎn)：簡單，可以自開始；提高精度時(shí)需要增加中間函數(shù)值；沒有充分利用前幾步得到的信息。多步法： yn-p , yn-p+1 , , yn-1 , yn yn +1612.4 線性多步法 單步法的優(yōu)缺點(diǎn)：簡單，可以自開始62yn-p, yn-p+1,

21、 yn-1, yn yn +1的方法考慮如下形式的求解公式此計(jì)算公式為線性的，所以稱為線性多步法。當(dāng) 時(shí)公式含有 ，這時(shí)公式是隱式的，而當(dāng) 時(shí)公式是顯式的62yn-p, yn-p+1, yn-1, yn 63常微分初值問題與積分公式等價(jià)基本思想：對(duì)f(x,y(x)進(jìn)行多項(xiàng)式插值，利用插值公式計(jì)算右邊積分，可以得到常微分初值問題求解公式。63常微分初值問題與積分公式等價(jià)基本思想：對(duì)f(x,y(x)64三次多項(xiàng)式插值求積分插值余項(xiàng)為令 ，假設(shè)已知 ，取插值節(jié)點(diǎn) ，對(duì)函數(shù) 作三次Lagrange多項(xiàng)式插值64三次多項(xiàng)式插值求積分插值余項(xiàng)為令 65四階Adams內(nèi)插公式內(nèi)插公式局部截?cái)嗾`差： 65四階

22、Adams內(nèi)插公式內(nèi)插公式局部截?cái)嗾`差： 66令 ，假設(shè)已知 ，取插值節(jié)點(diǎn) ，對(duì)函數(shù) 作三次Lagrange多項(xiàng)式插值外插公式局部截?cái)嗾`差： 66令 67內(nèi)插公式外插公式通常計(jì)算公式為（預(yù)估校正系統(tǒng)）： 特點(diǎn)：精度高，但須與同精度的單步法配合使用。67內(nèi)插公式外插公式通常計(jì)算公式為（預(yù)估校正系統(tǒng)）： 特點(diǎn)68例3：用Adams預(yù)測(cè)校正系統(tǒng)求初值問題其中，E=200, L=3, =100, =50解：由于Adams預(yù)測(cè)校正系統(tǒng)為四階公式，選用四階R-K公式計(jì)算開始值。a=t0=0; b=0.05I0=0;h=0.002f(t,I)=(200-50I3-100I)/368例3：用Adams預(yù)測(cè)校正

23、系統(tǒng)求初值問題其中，E=2069再利用Adams預(yù)測(cè)校正系統(tǒng)n45625t0.0080.0100.0120.050I0.4646910.5592060.644933.1.16486169再利用Adams預(yù)測(cè)校正系統(tǒng)n45625t0.008702.5 常微分方程組、高階微分方程及邊值問題的數(shù)值解 考慮如下含二個(gè)未知函數(shù)的方程組：一、常微分方程組的數(shù)值解 則：可依照常微分方程數(shù)值方法來求解。令：702.5 常微分方程組、高階微分方程及邊值問題的數(shù)值解71則其改進(jìn)的Euler格式具有預(yù)估形式校正公式為如對(duì)于方程組71則其改進(jìn)的Euler格式具有預(yù)估形式校正公式為如對(duì)于方程72兩方程微分方程組的四階R-K公式72兩方程微分方程組的四階R-K公式73二、高階微分方程的數(shù)值解 對(duì)于高階微分方程，可以把它化成微分方程