常微分方程模型及其數(shù)值解

1、常微分方程模型及其數(shù)值解第1頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三0、導(dǎo)言 在許多實(shí)際問題中，例如物理中的速率問題，人口的增長問題，放射性衰變問題，經(jīng)濟(jì)學(xué)中的邊際問題等，常常涉及到兩個(gè)變量之間的變化規(guī)律。微分方程是研究上述問題的一種機(jī)理分析方法，它在科技、工程、生態(tài)、環(huán)境、人口以及經(jīng)濟(jì)管理等領(lǐng)域中有著十分廣泛的應(yīng)用。在應(yīng)用微分方程解決實(shí)際問題時(shí)，必須經(jīng)過兩個(gè)階段。一是微分方程的建立，建立一個(gè)微分方程的實(shí)質(zhì)就是構(gòu)建函數(shù)、自變量以及函數(shù)對自變量的導(dǎo)數(shù)之間的一種平衡關(guān)系。而正確地構(gòu)建這種平衡關(guān)系，需要對實(shí)際問題的深入淺出的刻畫，根據(jù)物理的和非物理的原理、定律或定理，作出合理的假設(shè)和

2、簡化并將它升華成數(shù)學(xué)問題。第2頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三另一個(gè)是方程的求解和結(jié)果分析。對一些常系數(shù)的或特殊函數(shù)形式的微分方程，往往能得到解析解，這對實(shí)際問題的分析和應(yīng)用都是有利的,但是大多數(shù)變系數(shù)的、非線性函數(shù)形式的微分方程都是求不出解析解的,此時(shí)就需要應(yīng)用求解微分方程的另一個(gè)重要方法數(shù)值解法。本章簡要介紹有關(guān)微分方程模型的概念，微分方程的數(shù)值解法和圖解法，主要介紹若干建模實(shí)例，通過它們展示微分方程模型的建模步驟及解決實(shí)際問題的全過程。第3頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三1、實(shí)例及其數(shù)學(xué)模型 例1 海上緝私問題 海防某部緝私艇上的雷達(dá)發(fā)現(xiàn)正

3、東方向c海里處有一艘走私船正以速度a向正北方向行駛，緝私艇立即以最大速度b前往攔截。用雷達(dá)進(jìn)行跟蹤時(shí)，可保持緝私艇的速度方向始終指向走私船。建立任意時(shí)刻緝私艇的位置和緝私艇航線的數(shù)學(xué)模型，討論緝私艇能夠追上走私船的條件，求出追上的時(shí)間。第4頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三建立直角坐標(biāo)系如圖，設(shè)在t=0時(shí)刻緝私艇發(fā)現(xiàn)走私船，此時(shí)緝私艇的位置在(0, 0)，走私船的位置在(c, 0)。走私船以速度a平行于y軸正向行駛，緝私艇以速度b按指向走私船的方向行駛。在任意時(shí)刻t緝私艇位于P(x, y)點(diǎn)，而走私船到達(dá)Q(c, at)點(diǎn)，直線PQ與緝私艇航線（圖中曲線）相切，切線與x軸

4、正向夾角為。Q(c,at)P(x,y)R(c,y )0yxc第5頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三緝私艇在x, y方向的速度分別為 ，由直角三角形PQR寫出sin 和cos 的表達(dá)式，得到微分方程： （1） 初始條件為 （2） 這就是緝私艇位置(x(t), y(t)的數(shù)學(xué)模型。但是由方程（1）無法得到x(t), y(t)的解析解，需要用數(shù)值算法求解。我們將在后面繼續(xù)討論這個(gè)問題。 第6頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三例2 弱肉強(qiáng)食問題 自然界中在同一環(huán)境下的兩個(gè)種群之間存在著幾種不同的生存方式，比如相互競爭，即爭奪同樣的食物資源，造成一個(gè)種群趨于滅

5、絕，而另一個(gè)趨向環(huán)境資源容許的最大容量；或者相互依存，即彼此提供部分食物資源，二者和平共處，趨于一種平衡狀態(tài)；再有一種關(guān)系可稱之為弱肉強(qiáng)食，即某個(gè)種群甲靠豐富的自然資源生存，而另一種群乙靠捕食種群甲為生，種群甲稱為食餌（Prey），種群乙為捕食者（Predator），二者組成食餌-捕食者系統(tǒng)。海洋中的食用魚和軟骨魚（鯊魚等）、美洲兔和山貓、落葉松和蚜蟲等都是這種生存方式的典型。這樣兩個(gè)種群的數(shù)量是如何演變的呢？近百年來許多數(shù)學(xué)家和生態(tài)學(xué)家對這一系統(tǒng)進(jìn)行了深入的研究，建立了一系列數(shù)學(xué)模型，本節(jié)介紹的是最初的、最簡單的一個(gè)模型，它是意大利數(shù)學(xué)家Volterra在上個(gè)世紀(jì)20年代建立的。第7頁，共6

6、8頁，2022年，5月20日，5點(diǎn)21分，星期三模型 用x(t)表示時(shí)刻t食餌（如食用魚）的密度，即一定區(qū)域內(nèi)的數(shù)量，y(t)表示捕食者（如鯊魚）的密度。假設(shè)食餌獨(dú)立生存時(shí)的（相對）增長率為常數(shù)r0，即 ，而捕食者的存在使食餌的增長率減小，設(shè)減小量與捕食者密度成正比，比例系數(shù)為a0，則 。 捕食者離開食餌無法生存，設(shè)它獨(dú)自存在時(shí)死亡率為常數(shù)d0，即 ，而食餌的存在為捕食者提供了食物，使捕食者的死亡率減小，設(shè)減小量與食餌密度成正比，比例系數(shù)為b0，則 ，實(shí)際上，當(dāng)bxd時(shí)捕食者密度將增長。 給定食餌和捕食者密度的初始值x0, y0，由上可知x(t), y(t)滿足以下方程： （3）（3）的解x(

7、t), y(t)描述了食餌和捕食者密度隨時(shí)間的演變過程。但是我們同樣得不到x(t), y(t)的解析解，需要用數(shù)值算法求解。我們將在3繼續(xù)討論這個(gè)問題 第8頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三2 歐拉方法和龍格庫塔方法一階常微分方程初值問題的一般形式為 y=(x,y) ,axb(4)y(a)=其中(x,y)是已知函數(shù),為給定的初值. 如果函數(shù)(x,y)在區(qū)域axb,-y0為Lipschitz常數(shù),則初值問題(4)有唯一解.第9頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 所謂數(shù)值解法,就是設(shè)法將常微分方程離散化,建立差分方程,給出解在一些離散點(diǎn)上的近似值. a=

8、x0 x1x2xnxN=b其中剖分節(jié)點(diǎn)xn=a+nh,n=0,1,N, h稱為剖分步長.數(shù)值解法就是求精確解y(x)在剖分節(jié)點(diǎn)xn上的近似值yny(xn), n=1,2,n. 假設(shè)初值問題(4)的解y=y(x)唯一存在且足夠光滑.對求解區(qū)域a,b做剖分 我們采用數(shù)值積分方法來建立差分公式. 2.1 構(gòu)造數(shù)值解法的基本思想 在區(qū)間xn,xn+1上對方程(4)做積分,則有第10頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三對右邊的積分應(yīng)用左矩形公式，則有梯形公式oxyab左矩形公式y(tǒng)=(x)右矩形公式中矩形公式第11頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三對右邊的積分應(yīng)用

9、左矩形公式，則有因此，建立節(jié)點(diǎn)處近似值yn滿足的差分公式稱之為Euler公式. 稱為梯形公式. 若對(6)式右邊的積分應(yīng)用梯形求積公式，則可導(dǎo)出差分公式第12頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 利用Euler方法求初值問題 解 此時(shí)的Euler公式為稱為Euler中點(diǎn)公式或稱雙步Euler公式. 若在區(qū)間xn-1,xn+1上對方程(4)做積分,則有對右邊的積分應(yīng)用中矩形求積公式，則得差分公式例3的數(shù)值解.此問題的精確解是y(x)=x/(1+x2).第13頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三分別取步長h=0.2 ,0.1 ,0.05,計(jì)算結(jié)果如下第14頁

10、，共68頁，2022年，5月20日，5點(diǎn)21分，星期三hxnyny(xn)y(xn)-ynh=0.20.000.400.801.201.602.000.000000.376310.542280.527090.466320.406820.000000.344830.487800.491800.449440.400000.00000-0.03148-0.05448-0.03529-0.01689-0.00682h=0.10.000.400.801.201.602.000.000000.360850.513710.509610.458720.404190.000000.344830.487800.4

11、91800.449440.400000.00000-0.01603-0.02590-0.01781-0.00928-0.00419h=0.050.000.400.801.201.602.000.000000.352870.500490.500730.454250.402270.000000.344830.487800.491800.449440.400000.00000-0.00804-0.01268-0.00892-0.00481-0.00227第15頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三Euler中點(diǎn)公式則不然, 計(jì)算yn+1時(shí)需用到前兩步的值yn , yn-1 ,稱其為

12、兩步方法,兩步以上的方法統(tǒng)稱為多步法. 在Euler公式和梯形公式中,為求得yn+1,只需用到前一步的值yn,這種差分方法稱為單步法,這是一種自開始方法. 隱式公式中,每次計(jì)算yn+1都需解方程,要比顯式公式需要更多的計(jì)算量,但其計(jì)算穩(wěn)定性較好. 在Euler公式和Euler中點(diǎn)公式中,需要計(jì)算的yn+1已被顯式表示出來,稱這類差分公式為顯式公式,而梯形公式中,需要計(jì)算的yn+1隱含在等式兩側(cè),稱其為隱式公式.第16頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 從數(shù)值積分的角度來看,梯形公式計(jì)算數(shù)值解的精度要比Euler公式好,但它屬于隱式公式,不便于計(jì)算. 實(shí)際上,常將Euler

13、公式與梯形公式結(jié)合使用: 2.2 改進(jìn)的Euler方法第17頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 由迭代法收斂的角度看，當(dāng) (是給定的精度要求)時(shí), 取 就可以保證迭代公式收斂, 而當(dāng)h很小時(shí), 收斂是很快的. 而且, 只要 通常采用只迭代一次的算法:稱之為改進(jìn)的Euler方法. 這是一種單步顯式方法. 改進(jìn)的Euler方法也可以寫成第18頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 y=y-2x/y , 0 x1的數(shù)值解, 取步長h=0.1 . 精確解為y(x)=(1+2x)1/2.例4 求初值問題 y(0)=1 解 (1) 利用Euler方法第19頁，共68

14、頁，2022年，5月20日，5點(diǎn)21分，星期三計(jì)算結(jié)果如下: (2) 利用改進(jìn)Euler方法第20頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三nxnEuler方法yn改進(jìn)Euler法yn精確解y(xn)01234567891000.10.20.30.40.50.60.70.80.9111.11.1918181.2774381.3582131.4351331.5089661.5803381.6497831.7177791.78477011.0959091.1840961.2662011.3433601.4164021.4859561.5525151.6164761.6781681.7

15、3786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051第21頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 在節(jié)點(diǎn)xn+1的誤差y(xn+1)-yn+1 ,不僅與yn+1這一步計(jì)算有關(guān),而且與前n步計(jì)算值yn,yn-1,y1都有關(guān). 為了簡化誤差的分析,著重研究進(jìn)行一步計(jì)算時(shí)產(chǎn)生的誤差.即假設(shè)yn=y(xn),求誤差y(xn+1)-yn+1,這時(shí)的誤差稱為局部截?cái)嗾`差,它可以反映出差分公式的精度.2.3 差分公式的誤差分析 如果單步差分公式的局部截?cái)嗾`差為O(hp+1

16、),則稱該公式為p階方法.這里p為非負(fù)整數(shù).顯然,階數(shù)越高,方法的精度越高. 研究差分公式階的重要手段是Taylor展開式,一元函數(shù)和二元函數(shù)的Taylor展開式為:第22頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三另外,在yn=y(xn)的條件下,考慮到y(tǒng)(x)=(x,y(x),則有 y(xn)=(xn,y(xn)=(xn,yn)=n y(xn)=(xn,y(xn)=x(xn,yn)+y(xn,yn)(xn,yn)第23頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 yn+1=yn+h(xn,yn) 對Euler方法，有 =yn+(xn,yn)h+O(h2)從而有: y

17、(xn+1)-yn+1=O(h2)所以Euler方法是一階方法.再看改進(jìn)Euler方法, 因?yàn)榭傻玫?4頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三所以, 改進(jìn)的Euler方法是二階方法.而從而有: y(xn+1)-yn+1=O(h3)2.4 Taylor展開方法 設(shè)y(x)是初值問題（4)的精確解, 利用Taylor展開式可得第25頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三稱之為p階Taylor展開方法. 因此，可建立節(jié)點(diǎn)處近似值yn滿足的差分公式其中第26頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三所以,此差分公式是p階方法. 由于Taylor展開方法

18、涉及很多復(fù)合函數(shù)(x,y(x)的導(dǎo)數(shù)的計(jì)算,比較繁瑣,因而很少直接使用,經(jīng)常用它為多步方法提供初始值.然而, Taylor展開方法給出了一種構(gòu)造單步顯式高階方法的途徑. Euler方法可寫為 可見,公式的局部截?cái)嗾`差為: y(xn+1)-yn+1=O(hp+1).2.5 Runge-Kutta方法第27頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 構(gòu)造差分公式 改進(jìn)的Euler方法可寫為其中i,i,ij為待定參數(shù). 若此公式的局部截?cái)嗾`差為第28頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三由于 yn+1=yn+h1n+h2(n+hxn+hn yn)+O(h3)O(h3)

19、,稱此公式為p階Runge-kutta方法,簡稱p階R-K方法. 對于p=2的情形, 應(yīng)有 =yn+h(1+2)n+h22(xn+n yn)+O(h3)所以,只要令 1+2=1, 2=1/2, 2=1/2 (8)第29頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 一般地, 參數(shù)由(8)確定的一族差分公式(7)統(tǒng)稱為二階R-K方法.稱之為中點(diǎn)公式,或可寫為若取=1,則得1=2=1/2,=1,此時(shí)公式(7)就是改進(jìn)的Euler公式; 若取1=0,則得2=1,=1/2,公式(7)為 高階R-K公式可類似推導(dǎo). 下面列出常用的三階、四階R-K公式.第30頁，共68頁，2022年，5月20日

20、，5點(diǎn)21分，星期三 四階標(biāo)準(zhǔn)R-K公式 三階R-K公式第31頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 解 四階標(biāo)準(zhǔn)R-K公式為例3 用四階標(biāo)準(zhǔn)R-K方法求初值問題 y=y-2x/y , 0 x1 y(0)=1的數(shù)值解, 取步長h=0.2 .計(jì)算結(jié)果如下:第32頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三nxnyny(xn)nxnyny(xn)0120.00.20.41.001.18321.34171.001.18321.34163450.60.81.01.48331.61251.73211.48321.61251.7321 也可以構(gòu)造隱式R-K方法,其一般形式為稱

21、之為p級隱式R-K方法,同顯式R-K方法一樣確定參數(shù).如第33頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三是二級二階隱式R-K方法,也就是梯形公式.但是p級隱式R-K方法的階可以大于p,例如,一級隱式中點(diǎn)公式為或?qū)憺樗嵌A方法.2.6 變步長Runge-Kutta方法 以p階R-K方法為例討論.設(shè)從xn以步長h計(jì)算y(xn+1)的近似值為 ,局部截?cái)嗾`差為其中,C是與h無關(guān)的常數(shù).第34頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 如果將步長減半,取h/2為步長, 從xn經(jīng)兩步計(jì)算得到y(tǒng)(xn+1)的近似值記為 ,其局部截?cái)嗾`差為于是有從而,得到事后誤差估計(jì)可見,當(dāng)

22、成立時(shí),可取 .否則,應(yīng)將步長再次減半進(jìn)行計(jì)算.第35頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 求解初值問題的單步顯式方法可一統(tǒng)一寫為如下形式 yn+1=yn+h(xn,yn,h) （9) 對于Euler方法,有2.7 單步方法的收斂性 y=(x,y) ,axb y(a)= 其中(x,y,h)稱為增量函數(shù). (x,y,h)=(x,y)對于改進(jìn)的Euler方法,有 (x,y,h)=1/2(x,y)+(x+h,y+h(x,y)第36頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 設(shè)y(x)是初值問題（4)的解 ,yn是單步法 （9)產(chǎn)生的近似解.如果對任意固定的點(diǎn)xn,均

23、有y(xn),則稱單步法（9)是收斂的. 可見,若方法（9)是收斂的,則當(dāng)h0時(shí),整體截?cái)嗾`差en=y(xn)-yn將趨于零. 定理 設(shè)單步方法（9)是p1階方法, 增量函數(shù)(x,y,h)在區(qū)域axb,-yn)的變化均不超過 ,則稱此差分方法是絕對穩(wěn)定的. 討論數(shù)值方法的穩(wěn)定性,通常僅限于典型的試驗(yàn)方程 y=y 其中是復(fù)數(shù)且Re()0. 在復(fù)平面上,當(dāng)方法穩(wěn)定時(shí)要求變量h的取值范圍稱為方法的絕對穩(wěn)定域,它與實(shí)軸的交集稱為絕對穩(wěn)定區(qū)間. 第40頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 將Euler方法應(yīng)用于方程y=y, 得到 設(shè)在計(jì)算yn時(shí)產(chǎn)生誤差n,計(jì)算值yn=yn+n,則n將

24、對以后各節(jié)點(diǎn)值計(jì)算產(chǎn)生影響.記ym=ym+m ,mn,由上式可知誤差m滿足方程 m=(1+h)m-1=(1+h)m-nn , mn 對隱式單步方法也可類似討論.如將梯形公式用于方程y=y,則有 yn+1=yn+h/2 (yn+yn+1) yn+1=(1+h)yn 可見,若要|m|n|,必須且只須|1+h|1 ,因此Euler法的絕對穩(wěn)定域?yàn)閨1+h|1,絕對穩(wěn)定區(qū)間是-2Re()h0.解出yn+1得 第41頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三類似前面分析,可知絕對穩(wěn)定區(qū)域?yàn)橛捎赗e()0,所以此不等式對任意步長h恒成立,這是隱式公式的優(yōu)點(diǎn). 一些常用方法的絕對穩(wěn)定區(qū)間為方

25、法方法的階數(shù)穩(wěn) 定 區(qū) 間Euler方法梯形方法改進(jìn)Euler方法二階R-K方法三階R-K方法四階R-K方法122234(-2 , 0)(- , 0)(-2 , 0)(-2 , 0)(-2.51 , 0)(-2.78 , 0)第42頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 解 因y0=1,計(jì)算得y10=1024,而y(1)=9.35762310-14.例4 考慮初值問題 y=-30y , 0 x1 y(0)=1取步長h=0.1 ,利用Euler方法計(jì)算y10y(1). y(x)=e-30 x 這是因?yàn)閔=-3不屬于Euler方法的絕對穩(wěn)定區(qū)間. 若取h=0.01,計(jì)算得y100

26、=3.23447710-16. 若取h=0.001,計(jì)算得y1000=5.91199810-14. 若取h=0.0001,計(jì)算得y10000=8.94505710-14. 若取h=0.00001,計(jì)算得y100000=9.315610-14.第43頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 單步顯式方法的穩(wěn)定性與步長密切相關(guān), 在一種步長下是穩(wěn)定的差分公式,取大一點(diǎn)步長就可能是不穩(wěn)定的. 收斂性是反映差分公式本身的截?cái)嗾`差對數(shù)值解的影響;穩(wěn)定性是反映計(jì)算過程中舍入誤差對數(shù)值解的影響.只有即收斂又穩(wěn)定的差分公式才有實(shí)用價(jià)值.2.9 線性多步方法 由于在計(jì)算yn+1時(shí) ,已經(jīng)知道yn

27、 ,yn-1 ,及(xn,yn), (xn-1,yn-1),利用這些值構(gòu)造出精度高、計(jì)算量小的差分公式就是線性多步法.2.9.1 利用待定參數(shù)法構(gòu)造線性多步方法 r+1步線性多步方法的一般形式為第44頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三當(dāng)-10時(shí),公式為隱式公式,反之為顯式公式.參數(shù)i,i的選擇原則是使方法的局部截?cái)嗾`差為 y(xn+1)-yn+1=O(h)r+2 選取參數(shù),0,1,2,使三步方法 yn+1=yn+h(0n+1n-1+2n-2) 這里,局部截?cái)嗾`差是指 ,在yn-i=y(xn-i),i=0,1,r的前提下,誤差y(xn+1)-yn+1.為三階方法. 例5 解

28、 設(shè)yn=y(xn),yn-1=y(xn-1),yn-2=y(xn-2),則有 第45頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 n=(xn,y(xn)=y(xn) y(xn+1)=y(xn)+hy(xn)+1/2h2y(xn)+1/6h3y(xn) 于是有若使: y(xn+1)-yn+1=O(h4) ,只要,0,1,2滿足: n-1=(xn-1,y(xn-1)=y(xn-1)=y(xn-h) =y(xn)-hy(xn)+1/2h2y(xn)-1/6h3y(4)(xn)+O(h4) n-2=y(xn)-2hy(xn)+2h2y(xn)-4/3h3y(4)(xn)+O(h4) yn

29、+1=y(xn)+h(0+1+2)y(xn)-h2(1+22)y(xn) +h3(1/21+22)y(xn)-h4/6(1+82)y(4)(xn)+O(h5) +1/24h4y(4)(xn)+O(h5) 第46頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 =1, 0+1+2=1, 1+22=-1/2 , 1+42=1/3于是有三步三階顯式差分公式設(shè)pr(x)是函數(shù)(x,y(x)的某個(gè)r次插值多項(xiàng)式,則有解之得: yn+1=yn+h/12(23n-16n-1+5n-2) 因?yàn)?2.9.2 利用數(shù)值積分構(gòu)造線性多步方法 其中 第47頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期

30、三 選取不同的插值多項(xiàng)式pr(x),就可導(dǎo)出不同的差分公式.下面介紹常用的Adams公式. 設(shè)已求得精確解y(x)在步長為h的等距節(jié)點(diǎn)xn-r,xn上的近似值yn-r ,yn , 記k=(xk,yk) ,利用r+1個(gè)數(shù)據(jù)(xn-r,n-r),(xn,n)構(gòu)造r次Lagrange插值多項(xiàng)式由此,可建立差分公式 1.Adams顯式公式 其中 第48頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三由此,可建立差分公式 由于 hrj 則有 稱之為r+1步Adams顯式公式. 第49頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三下面列出幾個(gè)帶有局部截?cái)嗾`差主項(xiàng)的Adams顯式公式 r=

31、0 yn+1=yn+hn+(1/2)h2y(xn) 2.Adams隱式公式 r=1 yn+1=yn+(h/2)(3n-n-1)+(5/12)h3y(xn) r=2 yn+1=yn+(h/12)(23n-16n-1+5n-2) +(3/8)h4y(4)(xn) r=3 yn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3) +(251/720)h5y(5)(xn) 如果利用r+1個(gè)數(shù)據(jù)(xn-r+1,n-r+1),(xn+1,n+1)構(gòu)造r次Lagrange插值多項(xiàng)式pr(x),則可導(dǎo)出數(shù)值穩(wěn)定性好的隱式公式,稱為Adams隱式公式,其一般形式為第50頁，共68頁，2022年，

32、5月20日，5點(diǎn)21分，星期三其中系數(shù)為 下面列出幾個(gè)帶有局部截?cái)嗾`差主項(xiàng)的Adams隱式公式 r=0 yn+1=yn+hn+1-(1/2)h2y(xn) r=1 yn+1=yn+(h/2)(n+n+1)-(1/12)h3y(xn) r=2 yn+1=yn+(h/12)(5n+1+8n-n-1) -(1/24)h4y(4)(xn) r=3 yn+1=yn+(h/24)(9n+1+19n-5n-1+n-2) -(19/720)h5y(5)(xn)第51頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 3.Adams預(yù)估-校正公式 由顯式公式提供一個(gè)預(yù)估值，再用隱式公式校正一次，求得數(shù)值解

33、，稱為預(yù)估校正方法。 校正 yn+1=yn+(h/24)(9n+1+19n-5n-1+n-2) 一般預(yù)估公式和校正公式都采用同階公式。例如： 預(yù)估 yn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3) n+1=(xn+1,yn+1) , n=3,4,稱為四階Adams預(yù)估校正公式.實(shí)際計(jì)算時(shí)通常用四階單步方法(如四階R-K公式)為它提供起始值y1,y2,y3 . 例6 用四階Adams預(yù)估校正公式求解初值問題 第52頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三 y=y-2x/y , 0 x1 y(0)=1取步長h=0.1. 解 用四階R-K公式提供起始值,計(jì)算結(jié)

34、果如下xnR-k法yn預(yù)估值yn校正值yn精確值y(xn)00.10.20.30.40.50.60.70.80.9111.0954461.1832171.2649121.3415511.4140451.4830171.5489171.6121141.6729141.7315661.3416411.4142131.4832391.5491921.6124501.6733181.73204811.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051第53頁，共68頁，2022年，5月20日，5點(diǎn)

35、21分，星期三3 RK方法的MATLAB實(shí)現(xiàn)對于微分方程（組）的初值問題 龍格庫塔方法可用如下MATLAB命令實(shí)現(xiàn)其計(jì)算：t，x=ode23(f，ts，x0，options)t，x=ode45(f，ts，x0，options)其中ode23用的是3級2階龍格庫塔公式，ode45用的是以Runge-Kutta-Fehberg命名的 5級4階公式。第54頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三命令的輸入f是待解方程寫成的函數(shù)m文件：function dx=f(t,x)dx=f1;f2;fn;若ts=t0,t1,t2, ,tf，則輸出在指定時(shí)刻t0,t1,t2, ,tf的函數(shù)值；等分

36、點(diǎn)時(shí)用ts=t0:k:tf，輸出在t0,tf內(nèi)等分點(diǎn)處的函數(shù)值。x0為函數(shù)初值（n維向量）。options可用于設(shè)定誤差限（options缺省時(shí)設(shè)定相對誤差10-3，絕對誤差10-6），命令為： options=odeset(reltol,rt,abstol,at)其中rt，at分別為設(shè)定的相對誤差和絕對誤差。命令的輸出t為指定的ts，x為相應(yīng)的函數(shù)值（n維向量）。注意，計(jì)算步長是根據(jù)誤差限自動(dòng)調(diào)整的，并不是輸入中指定的ts的分點(diǎn)。第55頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三下面用MATLAB軟件解決1提出的兩個(gè)問題例1 海上緝私（續(xù)）模型的數(shù)值解 1. 設(shè)a=20 (海里/

37、小時(shí))，b=40 (海里/小時(shí))，c=15 (海里)，由模型（1），（2）求任意時(shí)刻緝私艇的位置及緝私艇航線。 對于給出的a, b, c用MATLAB求數(shù)值解時(shí)，記x(1)=x, x(2)=y, x=(x(1), x(2)T。編寫如下m文件：function dx=jisi(t,x) % 建立名為jisi的函數(shù)m文件a=20; b=40; c=15;s=sqrt(c-x(1)2+(a*t-x(2)2); dx=b*(c-x(1)/s;b*(a*t-x(2)/s; % 以向量形式表示方程（1）第56頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三然后運(yùn)行以下程序：ts=0:0.05:0.

38、5; % 設(shè)定t的起終點(diǎn)及中間的等分點(diǎn)，終點(diǎn)可先作試探，再按照x(t)c=15調(diào)整到0.5x0=0,0; % 輸入x，y的初始值（2）t,x=ode45(jisi,ts,x0); % 調(diào)用ode45計(jì)算t,x % 輸出t, x(t), y(t)plot(t,x), grid, % 按照數(shù)值輸出作x(t), y(t)的圖形gtext(x(t), gtext(y(t), pauseplot(x(:,1),x(:,2), grid, % 作y( x) 的圖形gtext(x), gtext(y)第57頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三得到的數(shù)值結(jié)果x(t), y(t)為緝私艇的位

39、置，列入表1。走私船的位置記作x1(t), y1(t)，顯然x1(t)= c=15，y1(t)=at=20t，將y1(t)列入表1最后一列?？芍?dāng)t=0.5（小時(shí)），x, y與x1, y1幾乎一致，認(rèn)為緝私艇追上走私船。x(t), y(t)及y(x)的圖形見圖2，y(x)為緝私艇的航線。tx(t)y(t)y1(t)00000.051.99840.06981.00.103.98540.29242.00.155.94450.69063.00.207.85151.28994.00.259.67052.11785.00.3011.34963.20056.00.3512.81704.55527.00.4

40、013.98066.17738.00.4514.74518.02739.00.5015.00469.997910.0a=20， b=40，c=15的數(shù)值解x(t), y(t)和y1(t) 第58頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三a=20，b=40，c=15 時(shí)x(t), y(t) 和y( x)的圖形 第59頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三2. 設(shè)b，c不變，而a變大為30，35，接近40 (海里/小時(shí))，觀察解的變化 修改a的輸入，并相應(yīng)地延長t的終點(diǎn)。設(shè)a=35，t的終點(diǎn)經(jīng)試探，調(diào)整為1.6合適。表2是計(jì)算結(jié)果，其中x(t), y(t)有兩列數(shù)字，左邊的是用“缺省”精度（即相對誤差10-3，絕對誤差10-6）計(jì)算的，中間的y1(t) =at=35t是走私船到達(dá)的位置。可知t=1.3, 1.4, 1.5時(shí)緝私艇的位置x15, 但y與y1(t)相差甚遠(yuǎn)，t=1.6時(shí)x, y與x1, y1也有差距，這是累積誤差造成的?？衫胦de45的控制參數(shù)options提高精度（上面的“調(diào)用ode45計(jì)算”用以下程序代替），如設(shè) 第60頁，共68頁，2022年，5月20日，5點(diǎn)21分，星期三opt=odeset(reltol,1e-6,abstol,1e-9);t,x=ode45(jisi,ts,x0,opt)