常微分方程模型及其數值解第1頁，課件共68頁，創作于2023年2月0、導言在許多實際問題中，例如物理中的速率問題，人口的增長問題，放射性衰變問題，經濟學中的邊際問題等，常常涉及到兩個變量之間的變化規律。微分方程是研究上述問題的一種機理分析方法，它在科技、工程、生態、環境、人口以及經濟管理等領域中有著十分廣泛的應用。在應用微分方程解決實際問題時，必須經過兩個階段。一是微分方程的建立，建立一個微分方程的實質就是構建函數、自變量以及函數對自變量的導數之間的一種平衡關系。而正確地構建這種平衡關系，需要對實際問題的深入淺出的刻畫，根據物理的和非物理的原理、定律或定理，作出合理的假設和簡化并將它升華成數學問題。第2頁，課件共68頁，創作于2023年2月另一個是方程的求解和結果分析。對一些常系數的或特殊函數形式的微分方程，往往能得到解析解，這對實際問題的分析和應用都是有利的,但是大多數變系數的、非線性函數形式的微分方程都是求不出解析解的,此時就需要應用求解微分方程的另一個重要方法──數值解法。本章簡要介紹有關微分方程模型的概念，微分方程的數值解法和圖解法，主要介紹若干建模實例，通過它們展示微分方程模型的建模步驟及解決實際問題的全過程。第3頁，課件共68頁，創作于2023年2月1、實例及其數學模型

例1海上緝私問題海防某部緝私艇上的雷達發現正東方向c海里處有一艘走私船正以速度a向正北方向行駛，緝私艇立即以最大速度b前往攔截。用雷達進行跟蹤時，可保持緝私艇的速度方向始終指向走私船。建立任意時刻緝私艇的位置和緝私艇航線的數學模型，討論緝私艇能夠追上走私船的條件，求出追上的時間。第4頁，課件共68頁，創作于2023年2月建立直角坐標系如圖，設在t=0時刻緝私艇發現走私船，此時緝私艇的位置在(0,0)，走私船的位置在(c,0)。走私船以速度a平行于y軸正向行駛，緝私艇以速度b按指向走私船的方向行駛。在任意時刻t緝私艇位于P(x,y)點，而走私船到達Q(c,at)點，直線PQ與緝私艇航線（圖中曲線）相切，切線與x軸正向夾角為

。Q(c,at)P(x,y)R(c,y)0yx

c第5頁，課件共68頁，創作于2023年2月緝私艇在x,y方向的速度分別為，由直角三角形PQR寫出sin

和cos

的表達式，得到微分方程：（1）

初始條件為（2）這就是緝私艇位置(x(t),y(t))的數學模型。但是由方程（1）無法得到x(t),y(t)的解析解，需要用數值算法求解。我們將在后面繼續討論這個問題。第6頁，課件共68頁，創作于2023年2月例2弱肉強食問題自然界中在同一環境下的兩個種群之間存在著幾種不同的生存方式，比如相互競爭，即爭奪同樣的食物資源，造成一個種群趨于滅絕，而另一個趨向環境資源容許的最大容量；或者相互依存，即彼此提供部分食物資源，二者和平共處，趨于一種平衡狀態；再有一種關系可稱之為弱肉強食，即某個種群甲靠豐富的自然資源生存，而另一種群乙靠捕食種群甲為生，種群甲稱為食餌（Prey），種群乙為捕食者（Predator），二者組成食餌-捕食者系統。海洋中的食用魚和軟骨魚（鯊魚等）、美洲兔和山貓、落葉松和蚜蟲等都是這種生存方式的典型。這樣兩個種群的數量是如何演變的呢？近百年來許多數學家和生態學家對這一系統進行了深入的研究，建立了一系列數學模型，本節介紹的是最初的、最簡單的一個模型，它是意大利數學家Volterra在上個世紀20年代建立的。第7頁，課件共68頁，創作于2023年2月模型用x(t)表示時刻t食餌（如食用魚）的密度，即一定區域內的數量，y(t)表示捕食者（如鯊魚）的密度。假設食餌獨立生存時的（相對）增長率為常數r>0，即，而捕食者的存在使食餌的增長率減小，設減小量與捕食者密度成正比，比例系數為a>0，則。捕食者離開食餌無法生存，設它獨自存在時死亡率為常數d>0，即，而食餌的存在為捕食者提供了食物，使捕食者的死亡率減小，設減小量與食餌密度成正比，比例系數為b>0，則，實際上，當bx>d時捕食者密度將增長。給定食餌和捕食者密度的初始值x0,y0，由上可知x(t),y(t)滿足以下方程：（3）（3）的解x(t),y(t)描述了食餌和捕食者密度隨時間的演變過程。但是我們同樣得不到x(t),y(t)的解析解，需要用數值算法求解。我們將在§3繼續討論這個問題第8頁，課件共68頁，創作于2023年2月§2歐拉方法和龍格—庫塔方法一階常微分方程初值問題的一般形式為y=(x,y),axb(4)y(a)=其中(x,y)是已知函數,為給定的初值.如果函數(x,y)在區域axb,-<y<上連續且關于y滿足Lipschitz條件其中L>0為Lipschitz常數,則初值問題(4)有唯一解.第9頁，課件共68頁，創作于2023年2月所謂數值解法,就是設法將常微分方程離散化,建立差分方程,給出解在一些離散點上的近似值.a=x0<x1<x2<…<xn<…<xN=b其中剖分節點xn=a+nh,n=0,1,…,N,h稱為剖分步長.數值解法就是求精確解y(x)在剖分節點xn上的近似值yny(xn),n=1,2,…,n.假設初值問題(4)的解y=y(x)唯一存在且足夠光滑.對求解區域[a,b]做剖分我們采用數值積分方法來建立差分公式.2.1構造數值解法的基本思想在區間[xn,xn+1]上對方程(4)做積分,則有第10頁，課件共68頁，創作于2023年2月對右邊的積分應用左矩形公式，則有梯形公式oxyab左矩形公式y=(x)右矩形公式中矩形公式第11頁，課件共68頁，創作于2023年2月對右邊的積分應用左矩形公式，則有因此，建立節點處近似值yn滿足的差分公式稱之為Euler公式.稱為梯形公式.若對(6)式右邊的積分應用梯形求積公式，則可導出差分公式第12頁，課件共68頁，創作于2023年2月利用Euler方法求初值問題

解此時的Euler公式為稱為Euler中點公式或稱雙步Euler公式.若在區間[xn-1,xn+1]上對方程(4)做積分,則有對右邊的積分應用中矩形求積公式，則得差分公式例3的數值解.此問題的精確解是y(x)=x/(1+x2).第13頁，課件共68頁，創作于2023年2月分別取步長h=0.2,0.1,0.05,計算結果如下第14頁，課件共68頁，創作于2023年2月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.491800.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頁，創作于2023年2月Euler中點公式則不然,計算yn+1時需用到前兩步的值yn,yn-1,稱其為兩步方法,兩步以上的方法統稱為多步法.在Euler公式和梯形公式中,為求得yn+1,只需用到前一步的值yn,這種差分方法稱為單步法,這是一種自開始方法.隱式公式中,每次計算yn+1都需解方程,要比顯式公式需要更多的計算量,但其計算穩定性較好.在Euler公式和Euler中點公式中,需要計算的yn+1已被顯式表示出來,稱這類差分公式為顯式公式,而梯形公式中,需要計算的yn+1隱含在等式兩側,稱其為隱式公式.第16頁，課件共68頁，創作于2023年2月從數值積分的角度來看,梯形公式計算數值解的精度要比Euler公式好,但它屬于隱式公式,不便于計算.實際上,常將Euler公式與梯形公式結合使用:2.2改進的Euler方法第17頁，課件共68頁，創作于2023年2月

由迭代法收斂的角度看，當(是給定的精度要求)時,取就可以保證迭代公式收斂,而當h很小時,收斂是很快的.而且,只要通常采用只迭代一次的算法:稱之為改進的Euler方法.這是一種單步顯式方法.改進的Euler方法也可以寫成第18頁，課件共68頁，創作于2023年2月y=y-2x/y,0x1的數值解,取步長h=0.1.[精確解為y(x)=(1+2x)1/2.]例4求初值問題y(0)=1

解(1)利用Euler方法第19頁，課件共68頁，創作于2023年2月計算結果如下:(2)利用改進Euler方法第20頁，課件共68頁，創作于2023年2月nxnEuler方法yn改進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.73786911.0954451.1832161.2649911.3416411.4142141.4832401.5491931.6124521.6733201.732051第21頁，課件共68頁，創作于2023年2月在節點xn+1的誤差y(xn+1)-yn+1,不僅與yn+1這一步計算有關,而且與前n步計算值yn,yn-1,…,y1都有關.為了簡化誤差的分析,著重研究進行一步計算時產生的誤差.即假設yn=y(xn),求誤差y(xn+1)-yn+1,這時的誤差稱為局部截斷誤差,它可以反映出差分公式的精度.2.3差分公式的誤差分析如果單步差分公式的局部截斷誤差為O(hp+1),則稱該公式為p階方法.這里p為非負整數.顯然,階數越高,方法的精度越高.研究差分公式階的重要手段是Taylor展開式,一元函數和二元函數的Taylor展開式為:第22頁，課件共68頁，創作于2023年2月另外,在yn=y(xn)的條件下,考慮到y(x)=(x,y(x)),則有y(xn)=(xn,y(xn))=(xn,yn)=ny(xn)=(xn,y(xn))=x(xn,yn)+y(xn,yn)(xn,yn)第23頁，課件共68頁，創作于2023年2月yn+1=yn+h(xn,yn)對Euler方法，有=yn+(xn,yn)h+O(h2)從而有:y(xn+1)-yn+1=O(h2)所以Euler方法是一階方法.再看改進Euler方法,因為可得第24頁，課件共68頁，創作于2023年2月所以,改進的Euler方法是二階方法.而從而有:y(xn+1)-yn+1=O(h3)2.4Taylor展開方法設y(x)是初值問題（4)的精確解,利用Taylor展開式可得第25頁，課件共68頁，創作于2023年2月稱之為p階Taylor展開方法.

……

……

……因此，可建立節點處近似值yn滿足的差分公式其中第26頁，課件共68頁，創作于2023年2月所以,此差分公式是p階方法.由于Taylor展開方法涉及很多復合函數(x,y(x))的導數的計算,比較繁瑣,因而很少直接使用,經常用它為多步方法提供初始值.然而,Taylor展開方法給出了一種構造單步顯式高階方法的途徑.Euler方法可寫為可見,公式的局部截斷誤差為:y(xn+1)-yn+1=O(hp+1).2.5Runge-Kutta方法第27頁，課件共68頁，創作于2023年2月構造差分公式

……

……改進的Euler方法可寫為其中i,i,ij為待定參數.若此公式的局部截斷誤差為第28頁，課件共68頁，創作于2023年2月由于yn+1=yn+h1

n+h2(n+hxn+hn

yn)+O(h3)O(h3),稱此公式為p階Runge-kutta方法,簡稱p階R-K方法.對于p=2的情形,應有=yn+h(1+2)n+h2

2(xn+n

yn)+O(h3)所以,只要令

1+2=1,2=1/2,2=1/2(8)第29頁，課件共68頁，創作于2023年2月一般地,參數由(8)確定的一族差分公式(7)統稱為二階R-K方法.稱之為中點公式,或可寫為若取=1,則得1=2=1/2,=1,此時公式(7)就是改進的Euler公式;若取1=0,則得2=1,==1/2,公式(7)為高階R-K公式可類似推導.下面列出常用的三階、四階R-K公式.第30頁，課件共68頁，創作于2023年2月四階標準R-K公式

三階R-K公式第31頁，課件共68頁，創作于2023年2月

解四階標準R-K公式為例3用四階標準R-K方法求初值問題y=y-2x/y,0x1y(0)=1的數值解,取步長h=0.2.計算結果如下:第32頁，課件共68頁，創作于2023年2月nxnyny(xn)nxnyny(xn)0120.00.20.41.001.18321.34171.001.18321.34163450.60.81.01.48331.61251.73211.48321.61251.7321也可以構造隱式R-K方法,其一般形式為稱之為p級隱式R-K方法,同顯式R-K方法一樣確定參數.如第33頁，課件共68頁，創作于2023年2月是二級二階隱式R-K方法,也就是梯形公式.但是p級隱式R-K方法的階可以大于p,例如,一級隱式中點公式為或寫為它是二階方法.2.6變步長Runge-Kutta方法以p階R-K方法為例討論.設從xn以步長h計算y(xn+1)的近似值為,局部截斷誤差為其中,C是與h無關的常數.第34頁，課件共68頁，創作于2023年2月如果將步長減半,取h/2為步長,從xn經兩步計算得到y(xn+1)的近似值記為,其局部截斷誤差為于是有從而,得到事后誤差估計可見,當成立時,可取.否則,應將步長再次減半進行計算.第35頁，課件共68頁，創作于2023年2月求解初值問題的單步顯式方法可一統一寫為如下形式yn+1=yn+h(xn,yn,h)（9)對于Euler方法,有2.7單步方法的收斂性y=(x,y),axby(a)=其中(x,y,h)稱為增量函數.(x,y,h)=(x,y)對于改進的Euler方法,有(x,y,h)=1/2[(x,y)+(x+h,y+h(x,y))]第36頁，課件共68頁，創作于2023年2月設y(x)是初值問題（4)的解,yn是單步法（9)產生的近似解.如果對任意固定的點xn,均有y(xn),則稱單步法（9)是收斂的.可見,若方法（9)是收斂的,則當h0時,整體截斷誤差en=y(xn)-yn將趨于零.

定理設單步方法（9)是p1階方法,增量函數(x,y,h)在區域axb,-<y<+,0hh0上連續,且關于y滿足Lipschitz條件,初始近似y0=y(a)=,則方法（9)是收斂的,且存在與h無關的常數C,使|y(xn)-yn|Chp

證明因為單步方法（9)是p階方法,則y(x)滿足定義第37頁，課件共68頁，創作于2023年2月其中,局部截斷誤差|Rn(h)|C1hp+1,記en=y(xn)-yn,則有利用Lipschitz條件得y(xn+1)=y(xn)+h(xn,y(xn),h)+Rn(h)遞推得到注意到en+1=en+h[(xn,y(xn),h)-(xn,yn,h)]+Rn(h)|en+1|(1+hL)|en|+C1hP+1

1+hLehL,(1+hL)nenhLeL(b-a)第38頁，課件共68頁，創作于2023年2月由于e0=y(a)-y0=0,所以有則有設(x,y)連續且關于y滿足Lipschitz條件,對于Euler方法,由于(x,y,h)=(x,y),故Euler方法是收斂的.對于改進的Euler方法,由(x,y)的Lipschitz條件有|en||e0|eL(b-a)+C1hp/L(eL(b-a)-1)|en|C1hp/L(eL(b-a)-1)=Chp|(x,y,h)-(x,y*,h)|1/2|(x,y)-(x,y*)|+1/2|(x+h,y+h(x,y))-(x+h,y*+h(x,y*))|1/2L(1+hL)|y-y*|則當hh0時,關于y滿足常數為1/2L(1+h0L)的Lipschitz第39頁，課件共68頁，創作于2023年2月條件,因此改進的Euler方法是收斂的.可類似驗證各階R-K方法是收斂的.2.8單步方法的穩定性

定義對于初值問題（4),取定步長h,用某個差分方法進行計算時,假設只在一個節點值yn上產生計算誤差,即計算值yn=yn+,如果這個誤差引起以后各節點值ym(m>n)的變化均不超過,則稱此差分方法是絕對穩定的.討論數值方法的穩定性,通常僅限于典型的試驗方程y=y其中是復數且Re()<0.在復平面上,當方法穩定時要求變量h的取值范圍稱為方法的絕對穩定域,它與實軸的交集稱為絕對穩定區間.第40頁，課件共68頁，創作于2023年2月將Euler方法應用于方程y=y,得到設在計算yn時產生誤差n,計算值yn=yn+n,則n將對以后各節點值計算產生影響.記ym=ym+m,mn,由上式可知誤差m滿足方程

m=(1+h)m-1=…=(1+h)m-n

n,mn對隱式單步方法也可類似討論.如將梯形公式用于方程y=y,則有yn+1=yn+h/2(yn+yn+1)

yn+1=(1+h)yn

可見,若要|m|<|n|,必須且只須|1+h|<1,因此Euler法的絕對穩定域為|1+h|<1,絕對穩定區間是-2<Re()h<0.解出yn+1得

第41頁，課件共68頁，創作于2023年2月類似前面分析,可知絕對穩定區域為由于Re()<0,所以此不等式對任意步長h恒成立,這是隱式公式的優點.一些常用方法的絕對穩定區間為方法方法的階數穩定區間Euler方法梯形方法改進Euler方法二階R-K方法三階R-K方法四階R-K方法122234(-2,0)(-,0)(-2,0)(-2,0)(-2.51,0)(-2.78,0)第42頁，課件共68頁，創作于2023年2月

解因y0=1,計算得y10=1024,而y(1)=9.35762310-14.例4考慮初值問題y=-30y,0x1y(0)=1取步長h=0.1,利用Euler方法計算y10y(1).[y(x)=e-30x]這是因為h=-3不屬于Euler方法的絕對穩定區間.若取h=0.01,計算得y100=3.23447710-16.若取h=0.001,計算得y1000=5.91199810-14.若取h=0.0001,計算得y10000=8.94505710-14.若取h=0.00001,計算得y100000=9.315610-14.第43頁，課件共68頁，創作于2023年2月單步顯式方法的穩定性與步長密切相關,在一種步長下是穩定的差分公式,取大一點步長就可能是不穩定的.收斂性是反映差分公式本身的截斷誤差對數值解的影響;穩定性是反映計算過程中舍入誤差對數值解的影響.只有即收斂又穩定的差分公式才有實用價值.2.9線性多步方法由于在計算yn+1時,已經知道yn,yn-1,…,及(xn,yn),(xn-1,yn-1),…,利用這些值構造出精度高、計算量小的差分公式就是線性多步法.2.9.1利用待定參數法構造線性多步方法r+1步線性多步方法的一般形式為第44頁，課件共68頁，創作于2023年2月當-10時,公式為隱式公式,反之為顯式公式.參數i,i的選擇原則是使方法的局部截斷誤差為y(xn+1)-yn+1=O(h)r+2

選取參數,0,1,2,使三步方法yn+1=yn+h(0

n+1

n-1+2

n-2)

這里,局部截斷誤差是指,在yn-i=y(xn-i),i=0,1,…,r的前提下,誤差y(xn+1)-yn+1.為三階方法.

例5

解設yn=y(xn),yn-1=y(xn-1),yn-2=y(xn-2),則有

第45頁，課件共68頁，創作于2023年2月

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+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頁，創作于2023年2月=1,0+1+2=1,1+22=-1/2,1+42=1/3于是有三步三階顯式差分公式設pr(x)是函數(x,y(x))的某個r次插值多項式,則有解之得:yn+1=yn+h/12(23n-16n-1+5n-2)因為2.9.2利用數值積分構造線性多步方法其中第47頁，課件共68頁，創作于2023年2月選取不同的插值多項式pr(x),就可導出不同的差分公式.下面介紹常用的Adams公式.設已求得精確解y(x)在步長為h的等距節點xn-r,…,xn上的近似值yn-r,…,yn,記k=(xk,yk),利用r+1個數據(xn-r,n-r),…,(xn,n)構造r次Lagrange插值多項式由此,可建立差分公式1.Adams顯式公式其中第48頁，課件共68頁，創作于2023年2月由此,可建立差分公式由于

hrj

則有稱之為r+1步Adams顯式公式.第49頁，課件共68頁，創作于2023年2月下面列出幾個帶有局部截斷誤差主項的Adams顯式公式r=0yn+1=yn+hn+(1/2)h2y(xn)2.Adams隱式公式r=1yn+1=yn+(h/2)(3n-n-1)+(5/12)h3y(xn)r=2yn+1=yn+(h/12)(23n-16n-1+5n-2)+(3/8)h4y(4)(xn)r=3yn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3)+(251/720)h5y(5)(xn)如果利用r+1個數據(xn-r+1,n-r+1),…,(xn+1,n+1)構造r次Lagrange插值多項式pr(x),則可導出數值穩定性好的隱式公式,稱為Adams隱式公式,其一般形式為第50頁，課件共68頁，創作于2023年2月其中系數為下面列出幾個帶有局部截斷誤差主項的Adams隱式公式r=0yn+1=yn+hn+1-(1/2)h2y

(xn)r=1yn+1=yn+(h/2)(n+n+1)-(1/12)h3y

(xn)r=2yn+1=yn+(h/12)(5n+1+8n-n-1)-(1/24)h4y(4)(xn)r=3yn+1=yn+(h/24)(9n+1+19n-5n-1+n-2)-(19/720)h5y(5)(xn)第51頁，課件共68頁，創作于2023年2月3.Adams預估-校正公式由顯式公式提供一個預估值，再用隱式公式校正一次，求得數值解，稱為預估校正方法。校正yn+1=yn+(h/24)(9n+1+19n-5n-1+n-2)一般預估公式和校正公式都采用同階公式。例如：預估yn+1=yn+(h/24)(55n-59n-1+37n-2-9n-3)

n+1=(xn+1,yn+1),n=3,4,…稱為四階Adams預估校正公式.實際計算時通常用四階單步方法(如四階R-K公式)為它提供起始值y1,y2,y3.例6用四階Adams預估校正公式求解初值問題第52頁，課件共68頁，創作于2023年2月y

=y-2x/y,0x1y(0)=1取步長h=0.1.解用四階R-K公式提供起始值,計算結果如下xnR-k法yn預估值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頁，創作于2023年2月§3R—K方法的MATLAB實現對于微分方程（組）的初值問題

龍格—庫塔方法可用如下MATLAB命令實現其計算：[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頁，創作于2023年2月命令的輸入f是待解方程寫成的函數m文件：functiondx=f(t,x)dx=[f1;f2;

;fn];若ts=[t0,t1,t2,…,tf]，則輸出在指定時刻t0,t1,t2,…,tf的函數值；等分點時用ts=t0:k:tf，輸出在[t0,tf]內等分點處的函數值。x0為函數初值（n維向量）。options可用于設定誤差限（options缺省時設定相對誤差10-3，絕對誤差10-6），命令為：options=odeset(‘reltol’,rt,’abstol’,at)其中rt，at分別為設定的相對誤差和絕對誤差。命令的輸出t為指定的ts，x為相應的函數值（n維向量）。注意，計算步長是根據誤差限自動調整的，并不是輸入中指定的ts的分點。第55頁，課件共68頁，創作于2023年2月下面用MATLAB軟件解決§1提出的兩個問題例1海上緝私（續）模型的數值解1.

設a=20(海里/小時)，b=40(海里/小時)，c=15(海里)，由模型（1），（2）求任意時刻緝私艇的位置及緝私艇航線。對于給出的a,b,c用MATLAB求數值解時，記x(1)=x,x(2)=y,x=(x(1),x(2))T。編寫如下m文件：functiondx=jisi(t,x)%建立名為jisi的函數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頁，創作于2023年2月然后運行以下程序：ts=0:0.05:0.5;%設定t的起終點及中間的等分點，終點可先作試探，再按照x(t)

c=15調整到0.5x0=[0,0];%輸入x，y的初始值（2）[t,x]=ode45(@jisi,ts,x0);%調用ode45計算[t,x]%輸出t,x(t),y(t)plot(t,x),grid,%按照數值輸出作x(t),y(t)的圖形gtext(‘x(t)’),gtext(‘y(t)’),pauseplot(x(:,1),x(:,2)),grid,%作y(x)的圖形gtext(‘x’),gtext(‘y’)第57頁，課件共68頁，創作于2023年2月得到的數值結果x(t),y(t)為緝私艇的位置，列入表1。走私船的位置記作x1(t),y1(t)，顯然x1(t)=c=15，y1(t)=at=20t，將y1(t)列入表1最后一列?？芍攖=0.5（小時），x,y與x1,y1幾乎一致，認為緝私艇追上走私船。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.4013.98066.17738.00.4514.74518.02739.00.5015.00469.997910.0a=20，b=40，c=15的數值解x(t),y(t)和y1(t)第58頁，課件共68頁，創作于2023年2月a=20，b=40，c=15時x(t),y(t)和y(x)的圖形

第59頁，課件共68頁，創作于2023年2月2.設b，c不變，而a變大為30，35，

接近40(海里/小時)，觀察解的變化修改a的輸入，并相應地延長t的終點。設a=35，t的終點經試探，調整為1.6合適。表2是計算結果，其中x(t),y(t)有兩列數字，左邊的是用“缺省”精度（即相對誤差10-3，絕對誤差10-6）計算的，中間的y1(t)=at=35t是走私船到達的位置?？芍猼=1.3,1.4,1.5時緝私艇的位置x

15,但y與y1(t)相差甚遠，t=1.6時x,y與x1,y1也有差距，這是累積誤差造成的?？衫胦de45的控制參數options提高精度（上面的“調用ode45計算”用以下程序代替），如設第60頁，課件共68頁，創作于2023年2月opt=odeset('reltol',1e-6,'abstol',1e-9);[t,x]=ode45(@jisi,ts,x0,opt);得到的x(t),y(t)，與走私船到達的位置x1(t)，y1(t)相對照，可知t=1.6時x,y與x1,y1幾乎一致，認為緝私艇追上走私船。x(t),y(t)及y