第四章 微分方程與差分方程方法

1、第四章 微分方程與差分方程方法第一節(jié) 微分方程模型我們在數(shù)學(xué)分析中所研究的函數(shù),是反映客觀現(xiàn)實(shí)世界運(yùn)動(dòng)過程中量與量之間的一種關(guān)系,但我們在構(gòu)造數(shù)學(xué)模型時(shí),遇到的大量實(shí)際問題往往不能直接寫出量與量之間的關(guān)系,卻能比較容易地建立這些變量和它們的導(dǎo)數(shù)(或微分)間的關(guān)系式,這種聯(lián)系著自變量、未知函數(shù)及其導(dǎo)數(shù)(或微分)的關(guān)系式稱為微分方程.§4.1.1 微分方程簡介這一節(jié),我們將介紹關(guān)于微分方程的一些基本概念.一、微分方程的階數(shù)首先我們具體的來看一個(gè)微分方程的例子.例4-1 物體冷卻過程的數(shù)學(xué)模型將某物體放置于空氣中,在時(shí)刻,測量得它的溫度為,10分鐘后測量得溫度為.我們要求決定此物體的溫度和

2、時(shí)間的關(guān)系,并計(jì)算20分鐘后物體的溫度.這里我們假定空氣的溫度保持為.解:根據(jù)物理學(xué)中的牛頓冷卻定律可知,熱量總是從溫度高的物體向溫度低的物體傳導(dǎo);一個(gè)物體的溫度變化速度與這一物體的溫度與其所在介質(zhì)溫度的差值成正比.設(shè)物體在時(shí)刻的溫度為,則溫度的變化速度可以用來表示.我們得到描述物體溫度變化的微分方程 (4.1.1)其中是比例常數(shù).方程(4.1.1)中含有未知函數(shù)及它的一階導(dǎo)數(shù),這樣的方程,我們稱為一階微分方程.微分方程中出現(xiàn)的未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)稱為微分方程的階數(shù).方程 (4.1.2)中未知函數(shù)最高階導(dǎo)數(shù)的階數(shù)是三階,則方程(4.1.2)稱為三階微分方程.二、常微分方程與偏微分方程如果在

3、微分方程中,自變量的個(gè)數(shù)只有一個(gè),我們稱這種微分方程為常微分方程;自變量的個(gè)數(shù)為兩個(gè)或兩個(gè)以上的微分方程稱為偏微分方程.方程 (4.1.3)就是偏微分方程的例子,其中是未知函數(shù),、都是自變量.而方程(4.1.1)(4.1.2)都是常微分方程的例子.三、線性與非線性微分方程如果階常微分方程 (4.1.4)的左端為關(guān)于未知函數(shù)及其各階導(dǎo)數(shù)的線性組合,則稱該方程為線性微分方程,否則稱為非線性方程.一般的階線性微分方程具有形式 (4.1.5)其中是關(guān)于的已知函數(shù).當(dāng)時(shí),稱(4.1.5)為階齊次線性微分方程;當(dāng)時(shí),稱(4.1.5)為階非齊次線性微分方程.例如,方程(4.1.2)是三階線性微分方程.而方程

4、 (4.1.6)是一階齊次非線性方程. 四、微分方程初邊值問題我們把含有個(gè)獨(dú)立的任意常數(shù)的解稱為階微分方程(4.1.4)的通解.為了確定微分方程一個(gè)特定的解,我們通常給出這個(gè)解所必須滿足的條件,這就是定解條件.常見的定解條件是初始條件.階微分方程(4.1.4)的初始條件是指如下的個(gè)條件:求解微分方程滿足初始條件的解,稱為初值問題.求解初值問題的過程,就是通過初始條件確定通解中的常數(shù),從而求得滿足初始條件的特解的過程.若方程所給出的定解條件,既有自變量初始時(shí)刻的值,也有自變量取終值時(shí)的值,則稱該問題為邊值問題.§4.1.2微分方程的求解及Matlab實(shí)現(xiàn)線性微分方程和低階特殊微分方程往

5、往可以通過解析解的方法求解,但一般的非線性微分方程是沒有解析解的,即使可以求得解析解,參數(shù)也很難確定,需要用數(shù)值解的方式求解.具體的求解方法很多,本節(jié)我們主要介紹如何利用matlab來求解微分方程的解析解和數(shù)值解.§4.1.2.1微分方程的解析解一、基本理論有些簡單的常微分方程可用一些技巧,如分離變量法、積分因子法、常數(shù)變異法、降階法等,化為可直接積分的方程從而求得顯示解.例如,一階常系數(shù)線性常微分方程可化為兩邊積分可得的通解為 .一般的常系數(shù)線性微分方程 (4.1.7)的解滿足疊加原理,即方程(4.1.7)的通解是對應(yīng)齊次方程 (4.1.8)的通解與非齊次方程(4.1.7)的一個(gè)特

6、解的和.一階常系數(shù)線性常微分方程總可用這一思路求得顯示解.高階線性常系數(shù)微分方程可用特征根法求對應(yīng)齊次微分方程(4.1.8)的通解.(4.1.8)對應(yīng)的代數(shù)特征方程若方程的特征根均可以求出,且兩兩相異,則方程(4.1.7)的解可以表示為其中,為待定系數(shù),是滿足方程(4.1.8)的一個(gè)特解,這個(gè)特解可以通過常數(shù)變異法,比較系數(shù)法得到.有重根的情況也可以寫出相應(yīng)的解析解形式,這里我們不做過多介紹,因?yàn)槲覀兿旅嬉榻B一種更簡單的計(jì)算機(jī)求解法.二、用matlab求解線性常系數(shù)微分方程解析解dsolve函數(shù)調(diào)用格式其中,既可以描述微分方程,也可以描述初邊值條件.在描述微分方程時(shí),用字母D代表微分,后面加

7、數(shù)字表示微分的階數(shù),例如D4y表示,也可以用D2y(0)=3 這類記號來表示初值條件 .一般而言任何D后面所跟的字母為因變量,自變量默認(rèn)為,也可以在函數(shù)調(diào)用時(shí)指明.例4-2 求方程的通解.解:輸入命令>>syms t,x; %定義符號變量x=dsolve(D2x-2*Dx-3*x=exp(-t) %求解方程 latex(x) %用 LATEX語句顯示結(jié)果運(yùn)行結(jié)果x =exp(-t)*C2+exp(3*t)*C1-1/4*t*exp(-t)ans =e3tit C2+e-tit C1-1/4,te-t原方程的通解為其中,為任意常數(shù).若給出初始或邊界條件,則可以通過這些條件建立方程,求

8、出的值.仍考慮上面的微分方程,假設(shè)已知 ,則可以通過下面的命令求滿足該方程的特解.>> x=dsolve(D2x-2*Dx-3*x=exp(-t),x(0)=3,Dx(0)=2)運(yùn)行結(jié)果x =27/16*exp(-t)+21/16*exp(3*t)-1/4*t*exp(-t)所以，滿足這一初值問題的特解是對于線性方程組我們可以通過定義多維的因變量同樣求解.例4-3 試求解下面的線性微分方程解:輸入命令>>x,y,z=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z&#

9、39;)；x=collect(simple(x) %化簡x的表達(dá)式y(tǒng)=collect(simple(y)z=collect(simple(z)結(jié)果為部分特殊的非線性微分方程也可以用dsolve()函數(shù)來求解析解.這樣的方程描述方式和前面介紹的線性微分方程是一致的.下面我們通過例子來演示非線性方程的解析解求解問題,同時(shí)還將演示不能求解的例子.例4-4 試求出一階非線性微分方程的解析解.解: 輸入命令>>x=dsolve('Dx=x*(1-x2)')運(yùn)行結(jié)果x = 1/(1+exp(-2*t)*C1)(1/2) -1/(1+exp(-2*t)*C1)(1/2)即該方程的

10、解析解為.但是如果稍微改變原方程,例如將等號右側(cè)加上1,則可以用下面語句試解方程.我們會(huì)發(fā)現(xiàn)該方程是沒有解析解的.>>x=dsolve('Dx=x*(1-x2)+1')運(yùn)行結(jié)果Warning: Explicit solution could not be found; implicit solution returned.> In dsolve at 312 x = t+Int(-1/(_a-_a3+1),_a = . x)+C1 = 0從上例我們可以看出,可以求得解析解的方程是很有限的,對于一般的非線性方程只能用數(shù)值解法去求解.下面我們介紹微分方程數(shù)值解的求

11、解方法.§4.1.2.2微分方程的數(shù)值解一、基本理論(一) 數(shù)值解的定義考慮一階常微分方程組初值問題 (4.1.9)其中 ,這里表示轉(zhuǎn)置.微分方程的數(shù)值解,是指不求出解的解析表達(dá)式,而是尋求在一系列離散結(jié)點(diǎn)上的近似值.稱為步長.通常,取為常數(shù),此時(shí)結(jié)點(diǎn)變?yōu)榈染嘟Y(jié)點(diǎn),有.(二) 求解數(shù)值解的方法1.歐拉法基本思想:在結(jié)點(diǎn)處用差商近似替代導(dǎo)數(shù)這樣導(dǎo)出歐拉格式的計(jì)算公式 (4.1.10)歐拉法可以求解各種形式的微分方程,但是它只有一階精度,所以實(shí)際應(yīng)用效率比較差.2.改進(jìn)的歐拉法基本思想: 使用數(shù)值積分離散化.對方程,兩邊由到積分,得到對右端的定積分用數(shù)值積分方法做離散化,可得計(jì)算公式,如

12、用矩形公式可得歐拉公式,若用梯形公式可得改進(jìn)的歐拉公式,下面我們用梯度公式來做數(shù)值積分 故有如下迭代計(jì)算公式: (4.1.11)實(shí)際應(yīng)用時(shí),與歐拉公式結(jié)合使用: (4.1.12)對于已給的精度,當(dāng)滿足時(shí),取,然后繼續(xù)下一步的計(jì)算.3.Runge-Kutta 方法基本思想:利用Taylor展式進(jìn)行離散化.假設(shè)(4.1.9)中的充分光滑,將在點(diǎn)作Taylor展開:其中對照標(biāo)準(zhǔn)形式.若取并用代替,則得到一個(gè)p階近似公式 (4.1.13)顯然p=1時(shí),式(4.1.13)就是(4.1.10),即為歐拉方法.但當(dāng)時(shí),如果直接利用公式(4.1.13)求解數(shù)值解問題,就需要計(jì)算的高階微商.這個(gè)計(jì)算量是很大的.

13、因此,直接利用式(4.1.13)構(gòu)造高階公式是不實(shí)用的.通常,在求解高階精度的數(shù)值解時(shí),我們并不是直接使用Taylor級數(shù),而是利用它的思想,即計(jì)算在不同結(jié)點(diǎn)的函數(shù)值,然后作這些函數(shù)值的線性組合,構(gòu)造近似公式,式中有一些可供選擇的參數(shù).將近似公式與Taylor展開式相比較,使前面的若干項(xiàng)密合,從而使近似公式達(dá)到一定的精度.若線性組合選取的函數(shù)值個(gè)數(shù)為,計(jì)算達(dá)到的精度階數(shù)為,這樣的Runge-Kutta 方法就是我們通常所說的p階n級Runge-Kutta方法.后來,德國學(xué)者Felhberg對傳統(tǒng)的Runge-Kutta 方法進(jìn)行了改進(jìn),使原來的定步長算法變?yōu)樽詣?dòng)變換步長的Runge-Kutta

14、-Felhberg方法,保證了更高的精度和數(shù)值穩(wěn)定性,它已經(jīng)成為求解數(shù)值解問題中最常見的方法.具體的算法我們在這里不做介紹,下面我們主要來討論如何通過matlab軟件來運(yùn)用這一方法求解微分方程的數(shù)值解問題.二、用matlab求解微分方程數(shù)值解(一) 一階微分方程組初值問題Matlab下求解一階微分方程組初值問題(4.1.9)的數(shù)值解的最常見的方法ode45()函數(shù),該函數(shù)實(shí)現(xiàn)了前面介紹的變步長四階五級Runge-Kutta-Felhberg算法,可以采用變步長的算法求解微分方程.ode45()函數(shù)的調(diào)用格式參數(shù)說明t, y 中t 表示自變量,是標(biāo)量.y表示函數(shù),可以是標(biāo)量也可以是向量.表示初值

15、.當(dāng) y為 n維向量時(shí),表明求解的問題是有n個(gè)未知函數(shù)的方程組,則也應(yīng)為n維向量. 表示微分方程的求解區(qū)間,即自變量的取值范圍.如果只給出一個(gè)值 則表示初始時(shí)刻為.另外,該函數(shù)還允許,即可以認(rèn)為為終值時(shí)刻,為初始時(shí)刻,則表示狀態(tài)變量的終值,而該函數(shù)可以直接求解這樣一個(gè)終值問題.odefun 表示待解的微分方程(4.1.9)中所寫成的m-文件的名稱,m-文件的編寫格式為 function ydot=odefun(t,y) ydot=f(t,y);也可以不把函數(shù)單獨(dú)用m-文件保存,而是直接用inline()函數(shù)輸入,從而獲得函數(shù)名,編寫格式為 odefun=inline(f(t,y),t,y);其

16、中,t是自變量,即使所求方程是自治系統(tǒng),也需要寫出t,否則matlab在變量傳遞過程中將出現(xiàn)問題.y為未知函數(shù),ydot表示未知函數(shù)的導(dǎo)數(shù).如果求解的問題是微分方程組問題,即我們之前討論過的y為向量的情況,則在書寫上述m-文件或inline()函數(shù)時(shí),待解方程應(yīng)以ydot的分量形式寫成,或者仍用單個(gè)方程的編寫格式,多個(gè)方程用 括起來,并以;隔開,在后面的例題中我們將具體的看到.options 表示控制選項(xiàng).在微分方程求解中有時(shí)需要對求解算法及控制條件進(jìn)行進(jìn)一步設(shè)置,這可以通過求解過程中的options變量來進(jìn)行修改.初始o(jì)ptions變量可以通過odeset()函數(shù)來獲取,該變量是一個(gè)結(jié)構(gòu)體變

17、量,其中有眾多成員變量,表4-1列出了常用的一些成員變量.表4-1 微分方程求解函數(shù)的控制參數(shù)表參數(shù)名參數(shù)說明RelTol為相對誤差容許上限,默認(rèn)值為,在一些特殊的微分方程求解時(shí),為了保證較高的精度,可以適當(dāng)減小該值.AbsTol為絕對誤差容許上限,默認(rèn)值為,它是一個(gè)向量,每個(gè)分量表示各未知函數(shù)分量允許的絕對誤差.MaxStep為求解方程允許的最大步長Mass微分代數(shù)方程中的質(zhì)量矩陣,用于描述微分代數(shù)方程Jacobian為描述Jacobi矩陣函數(shù)的函數(shù)名,如果已知Jacobi矩陣,可以加速仿真過程修改option變量有兩種方式:用odeset()函數(shù)設(shè)置,其格式為 options=odeset

18、(RelTol,1e-7);達(dá)到的目的是,將相對誤差上限設(shè)置為較小的直接修改options的成員變量,也可以達(dá)到這樣的目的,其格式為 options=odeset; options.RelTol=1e-7;表示附加參數(shù).即在函數(shù)表達(dá)式中,除了還有其他的可變的參數(shù),引入這樣的附加參數(shù),使得微分方程的某些參數(shù)可以選擇不同的值,對不同值求解時(shí),考慮附加參數(shù)可以避免每次修改模型文件.但是值得注意的是,這種情況下的函數(shù)定義格式有特殊要求,格式為 function ydot=odefun(t,y,flag,p1,p2,pm) ydot=f(t,y,p1,p2,pm);或 odefun=inline(f(t

19、,y,p1,p2, ,pm),t,y,flag,p1,p2,pm);注意,這里的flag變量不能省略,用于占位.例4-5 解微分方程組解: 該方程是非線性微分方程,不能用我們之前討論的求解解析解的方法來處理,只能用數(shù)值解法求解.首先,我們將待解方程寫成m-文件,用文件名eg1fun.m保存.我們將2個(gè)未知變量x,y,以分量的形式寫成一個(gè)向量變量x.%函數(shù)eg1fun.mfunction f=eg1fun(t,x)f(1)=-x(1)3-x(2);f(2)=x(1)-x(2)3;f=f(:); %保證f為列向量這個(gè)m-文件還有另一種等價(jià)的書寫方式.%函數(shù)eg1fun.m的另一種書寫方式funct

20、ion f=eg1fun(t,x)f=-x(1)3-x(2);x(1)-x(2)3;這樣我們就完成了對待解方程的描述,然后輸入命令>>t,x=ode45('eg1fun',0 30,1;0.5); %求解方程 subplot(1,2,1); %分2部分繪圖,編輯第一部分 plot(t,x(:,1),t,x(:,2),':'); %在第一個(gè)圖中繪制x,y的函數(shù)圖象title('函數(shù)曲線圖'); %圖象名稱 xlabel('t'); %橫軸說明 legend('x(t)','y(t)');

21、%說明圖例 axis square %顯示正方形的坐標(biāo)系subplot(1,2,2); %編輯第二部分圖象plot(x(:,1),x(:,2); %在第二個(gè)圖中繪制相平面圖title('相平面圖'); %圖象名稱xlabel('x'); %橫軸說明ylabel('y'); %縱軸說明axis square %顯示正方形的坐標(biāo)系運(yùn)行結(jié)果:圖4-1例4-5函數(shù)曲線和相平面圖運(yùn)行程序,得到函數(shù)圖象和相平面圖.根據(jù)相平面上的相軌線,可以看出解從初值開始的變化趨勢,這對于研究微分方程平衡點(diǎn)的穩(wěn)定性是很有意義的,在后面的章節(jié)中我們將具體討論.例4-6 解微分

22、方程組解:此問題仍為非線性微分方程組,我們求它的數(shù)值解,不妨將求解區(qū)間取為.下面我們用inline()函數(shù)來描述方程組,這樣就不用單獨(dú)編寫m-文件來描述方程.輸入命令:>>eg2fun=inline('y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)','t','y');tf=20;y0=0;1;1;t,y=ode45(eg2fun,0,tf,y0);plot(t,y(:,1),'-',t,y(:,2),'*',t,y(:,3),'+'); %繪出解的圖象lege

23、nd('y1','y2','y3'); %說明圖例運(yùn)行得到函數(shù)曲線圖圖4-2 例4-6函數(shù)曲線圖例4-7 編寫帶有附加參數(shù)的matlab函數(shù),來描述食餌-捕食者模型并對和分別求解,畫出解的曲線圖和相軌線圖.解:選定附加參數(shù)為,編寫如下m-文件來描述題述方程%函數(shù)eg3fun.mfunction xdot=eg3fun(t,x,flag,a1,b1,a2,b2)xdot=x(1)*(a1-b1*x(2);-x(2)*(a2-b2*x(1);注意,這里的flag變量不能省略.這樣我們就完成了對方程的描述.我們不妨取自變量區(qū)間為0,10,對方程求數(shù)值解

24、.輸入命令>>tf=10;x0=1;1;a1=3;b1=2;a2=2.5;b2=1;t,x=ode45(eg3fun,0,tf,x0,a1,b1,a2,b2);subplot(1,2,1);plot(t,x(:,1),t,x(:,2),':');title('函數(shù)曲線圖');xlabel('t');legend('x1','x2');axis square;subplot(1,2,2);plot(x(:,1),x(:,2);title('相平面圖');xlabel('x1'

25、;);ylabel('x2');axis square; 運(yùn)行結(jié)果圖4-3 當(dāng)時(shí),只需要修改上述命令中的給a1,b1,a2,b2的賦值語句,就可以得到需要的結(jié)果.圖4-4 (二) 一階微分方程組邊值問題對于一階常微分方程組的邊值問題 (4.1.14)這里y,f,g都可以是向量.用matlab求(4.1.14)的數(shù)值解的基本格式 sinit=bvpinit(tint,yinit);由在粗略結(jié)點(diǎn)tinit的預(yù)估解yinit生成粗略解網(wǎng)絡(luò)sinit. sol=bvp4c(odefun,bcfun,sinit)其中odefun是微分方程組函數(shù),bcfun為邊值條件函數(shù),sinit是由b

26、vpinit得到的粗略解網(wǎng)絡(luò).求得的邊值問題的解sol是一個(gè)結(jié)構(gòu)體變量,sol.x為求解結(jié)點(diǎn),sol.y為數(shù)值解. sx=deval(sol,ti)計(jì)算由bvp4c得到的解在ti的值.例4-8 求解邊值問題解:輸入命令>>sinit=bvpinit(0:4,1;0);odefun=inline('y(2);-abs(y(1)','t','y');bcfun=inline('ya(1);yb(1)+2','ya','yb');sol=bvp4c(odefun,bcfun,sinit);t=

27、linspace(0,4,101);y=deval(sol,t);plot(t,y(1,:),sol.x,sol.y(1,:),'o',sinit.x,sinit.y(1,:),'x');legend('解曲線','求解點(diǎn)','粗略解');運(yùn)行結(jié)果圖4-5 例4-8函數(shù)曲線圖(三) 高階微分方程由于matlab只能處理一階微分方程組問題,對于高階微分方程的初邊值問題,我們必須通過變量代換,使它轉(zhuǎn)換為一階方程組,才能用matlab求解數(shù)值解.我們考慮高階常微分方程的初值問題 (4.1.15)且初始值為已知.令,通過這

28、樣的變量代換,高階微分方程(4.1.15)就轉(zhuǎn)化為一階方程組的形式 (4.1.16)且初值條件為.這樣,我們就可以通過之前討論過的一階微分方程組的數(shù)值解法來研究這一問題. 對于高階微分方程組,只要依次將每個(gè)未知變量的各階導(dǎo)數(shù)均單獨(dú)定義變量,就可以同樣的轉(zhuǎn)化為一階微分方程組.例4-9 用數(shù)值解法求解微分方程解:首先我們需要將這個(gè)2階微分方程轉(zhuǎn)化為一階微分方程組,我們令,則原方程轉(zhuǎn)換為對0,20區(qū)間求數(shù)值解,輸入命令>>tf=20;x0=-0.2;-0.7;eg5fun=inline('x(2);-(x(1)2-1)*x(2)-x(1)','t',

29、9;x');t,y=ode45(eg5fun,0,tf,x0);plot(t,y);legend('y','y''');運(yùn)行得到函數(shù)圖象圖4-6 例4-9函數(shù)曲線圖例4-10 求解微分方程組(豎直加熱板的自然對流問題)解:這是一個(gè)高階方程組問題,我們需要通過變量代換轉(zhuǎn)化為一階微分方程組,從而尋求數(shù)值解.我們令原方程轉(zhuǎn)化為如下一階方程組輸入命令>>y0=0;0;0.68;1;-0.5;tf=5;eg6fun=inline('y(2);y(3);-3*y(1)*y(3)+2*y(2)2-y(4);y(5);-2.1*y(1

30、)*y(5)','t','y');t,y=ode45(eg6fun,0,tf,y0);plot(t,y(:,1),t,y(:,4),':');legend('f','T');運(yùn)行得到函數(shù)圖象圖4-7 例4-10函數(shù)曲線圖§4.1.3導(dǎo)彈追蹤問題一、導(dǎo)彈追擊問題一艘敵艦在某海域內(nèi)沿正北方向航行時(shí),我方戰(zhàn)艦恰位于敵艦的正西方向1km處.我艦向敵艦發(fā)射導(dǎo)彈,導(dǎo)彈頭始終對準(zhǔn)敵艦,敵艦速度為=1km/min,導(dǎo)彈速度為敵艦速度的5倍.問需要多長時(shí)間,在何處導(dǎo)彈擊中敵艦?以我艦位置為坐標(biāo)原點(diǎn),正北方向?yàn)閥軸方

31、向建立坐標(biāo)系,設(shè)t時(shí)刻導(dǎo)彈所處的位置為,敵艦所處的位置為.圖4-8 導(dǎo)彈追擊問題由于導(dǎo)彈頭始終對準(zhǔn)敵艦,因此直線PQ是導(dǎo)彈運(yùn)行軌跡OP在P點(diǎn)的切線,即從而 (4.1.17)又因?yàn)閷?dǎo)彈速度為敵艦速度的5倍,所以O(shè)P弧的長度為AQ長度的5倍.即 (4.1.18)聯(lián)立(4.1.17)和(4.1.18)得到兩邊對x求導(dǎo),得到微分方程 (4.1.19)該問題的初始條件為首先,我們將二階方程化為一階方程組,令,(4.1.19)轉(zhuǎn)化為 (4.1.20)對于這一問題,我們可以求得方程的解析解,從而確定導(dǎo)彈擊中敵艦位置和時(shí)間,也可以通過求數(shù)值解的方法得到同樣的結(jié)論.1) 解析解.關(guān)于的方程分離變量,有兩邊積分,

32、并代入初始條件,得即 (4.1.21)又因?yàn)?(4.1.22)(4.1.21)和(4.1.22)兩式相加,得到下面,我們求解關(guān)于的方程直接積分,并代入初始條件得到當(dāng)時(shí),即當(dāng)敵艦航行到點(diǎn)(1,)處時(shí),被導(dǎo)彈擊中.被擊中的時(shí)間為.2) 數(shù)值解.首先,編寫出描述方程(4.1.20)的m-文件function dy=equ20(x,y)dy=y(2);1/5*sqrt(1+y(2)2)/(1-x);由于方程在處沒有定義,我們?nèi)∏蠼鈪^(qū)間為0,1-1e-8,輸入命令xf=1-1e-8;x,y=ode45('equ20',0,xf,0;0)plot(x,y(:,1);hold ony=0:0

33、.01:2;plot(1,y,'*');運(yùn)行結(jié)果y = 0 00.0000 0.0001 0.2083 16.1097 0.2083 20.0688圖4-9 導(dǎo)彈擊中曲線說明當(dāng)時(shí),導(dǎo)彈擊中敵艦的位置大致在(1,0.2083)處,擊中的時(shí)間=0.2083min=12.498s.這與求解析解得到的結(jié)果是一致的.比較上面2種求解方法,我們發(fā)現(xiàn),解析解可以得到解的真實(shí)值,但是方法比較復(fù)雜,技巧性強(qiáng),只能處理比較簡單的方程問題;而數(shù)值解雖然得到的解與真實(shí)值之間會(huì)有一些允許范圍內(nèi)的誤差,但是方法比較容易掌握,更具一般性.二、導(dǎo)彈系統(tǒng)的改進(jìn)現(xiàn)根據(jù)情報(bào),這種敵艦?zāi)茉谖遗灠l(fā)射導(dǎo)彈后T分鐘做出反應(yīng)

34、并摧毀導(dǎo)彈.因此,我們要求改進(jìn)電子導(dǎo)彈系統(tǒng),使其根據(jù)敵艦與我艦的距離,行使方向和速度,能自動(dòng)判斷出敵艦是否在有效打擊范圍之內(nèi).對于更一般的情況,我們做出如下假設(shè),設(shè)敵艦在我艦正東方向d km處,行駛速度為 km/min,行駛方向與正東方向的夾角為,導(dǎo)彈的飛行速度為v km/min.問題的關(guān)鍵是計(jì)算出導(dǎo)彈擊中敵艦所需要的時(shí)間,并將與T比較,若<T,則敵艦在打擊范圍內(nèi).我們?nèi)砸晕遗炍恢脼樽鴺?biāo)原點(diǎn),以正北方向?yàn)閥軸建立坐標(biāo)系,設(shè)t時(shí)刻導(dǎo)彈所處的位置為,敵艦所處位置為.圖4-10 導(dǎo)彈擊中問題由于導(dǎo)彈頭始終對準(zhǔn)敵艦,因此直線PQ是導(dǎo)彈運(yùn)行軌跡OP在P點(diǎn)的切線,即 (4.1.23)又因?yàn)閷?dǎo)彈速度為

35、敵艦速度的倍,所以O(shè)P弧的長度為AQ長度的倍.即 (4.1.24)聯(lián)立方程(4.1.23)和(4.1.24),消去,再對方程兩邊對x求導(dǎo),得到微分方程 (4.1.25)同樣,我們令,(4.1.25)轉(zhuǎn)化為一階微分方程組 (4.1.26)從圖象上來判斷,我們知道,當(dāng)和兩點(diǎn)的運(yùn)動(dòng)曲線相遇時(shí),導(dǎo)彈擊中敵艦.因此我們可以認(rèn)為,若滿足方程(4.1.25)且 (4.1.27)則點(diǎn)()為導(dǎo)彈擊中敵艦的擊中點(diǎn).再根據(jù)Q點(diǎn)的表達(dá)式,可以計(jì)算出擊中時(shí)間.若<T,則敵艦在打擊范圍內(nèi),可以發(fā)射.從上述分析,我們知道,當(dāng)、v、d、T已知時(shí),根據(jù)(4.1.25)和(4.1.27)可以計(jì)算出導(dǎo)彈擊中敵艦所需要的時(shí)間,

36、且當(dāng)<T時(shí),敵艦在打擊范圍內(nèi).這種對于每組、v、d、T分別計(jì)算從而作出判斷的方法,稱為在線算法.如果在敵艦和我艦的情況已知,即、v、T已知,計(jì)算出所有在打擊范圍內(nèi)的d和,之后要判斷敵艦是否在打擊范圍,只需要在計(jì)算結(jié)果中查詢,這種算法稱為離線算法.從原理上來看,我們發(fā)現(xiàn),在線算法靈活,容易調(diào)整參數(shù)和模型,但是速度比較慢;離線算法事先計(jì)算好,實(shí)時(shí)使用查詢方式,不需要計(jì)算,速度快.(1)在線算法現(xiàn)已知=90km/h,v=450km/h,d=50km,T=0.1h.首先,我們編寫帶參數(shù)的m-文件描述方程(4.1.26).為了使方程始終有意義,我們通常在分母上加上一個(gè)無窮小量1e-8.%方程(4.

37、1.26)function dy=equ26(x,y,v0,v,d,theta)dy(1)=y(2);dy(2)=(v0/v)*(1+(y(2)2)(1/2)*(sin(theta)-y(2)*cos(theta)2)/(d-x)*(sin(theta)-y(2)*cos(theta)+cos(theta)*(y(2)*(d-x)+y(1)+1e-8);dy=dy(:);輸入命令clear;close;v0=90;v=450;d=50;theta=pi/4;T=0.1;%初始化參數(shù)x,y=ode45(equ26,0,60,0;0,v0,v,d,theta);%計(jì)算導(dǎo)彈運(yùn)行軌跡方程的數(shù)值解z=(

38、x-d)*tan(theta);%計(jì)算敵艦運(yùn)行的軌跡n=length(x);for i=1:n if abs(z(i)-y(i,1)<1e-6 xk=x(i); yk=y(i,1); break; endendxkyk%求出擊中點(diǎn)xk,yk坐標(biāo)for i=1:n if z(i)>0 z1(i)=z(i); else z1(i)=0; endendplot(x,y(:,1),x,z1(:),xk,yk,'*');legend('導(dǎo)彈運(yùn)行軌跡','敵艦運(yùn)行軌跡','擊中點(diǎn)');tk=yk/(v0*sin(theta)+1e

39、-8)%計(jì)算擊中時(shí)間if tk<Tdisp('敵艦在打擊范圍內(nèi)，擊中地點(diǎn)在', num2str(xk),num2str(yk),'擊中時(shí)間為', num2str(tk);else disp('敵艦不在打擊范圍內(nèi)');end得到結(jié)果擊中點(diǎn) xk=58.4056,yk=8.4056,擊中時(shí)間 tk =0.1321,所以敵艦不在打擊范圍內(nèi),應(yīng)等接近一些再發(fā)射.圖4-11 在線算法(2)離線算法在介紹離線算法之前,我們先來看方程(4.1.26)的另一種表達(dá).因?yàn)閷?dǎo)彈的線速度 (4.1.28)方程(4.1.28)與(4.1.23)聯(lián)立,可以建立以t為

40、參數(shù)的關(guān)于x,y的參數(shù)方程 (4.1.29)當(dāng)滿足,則導(dǎo)彈擊中了敵艦.由于我們所要求的打擊范圍是在中討論的,所以考慮以t為參數(shù)的方程(4.1.29)能使求解過程大大簡化.現(xiàn)已知=90km/h,v=450km/h,T=0.1h,要計(jì)算出所有在打擊范圍內(nèi)的d和.依題意,.據(jù)此,我們可以確定d的取值范圍.當(dāng)時(shí),敵艦正好背向行駛,導(dǎo)彈直線運(yùn)行,擊中時(shí)間求得.當(dāng)時(shí),敵艦迎面駛來,導(dǎo)彈直線運(yùn)行,擊中時(shí)間則.所以.這樣,我們對于所有可能的和的取值,計(jì)算擊中所需時(shí)間,從而對不同的,得到的臨界值.具體應(yīng)用時(shí)直接查詢判斷.編寫如下m-文件描述方程(4.1.29)%方程(4.1.29)function dy=equ

41、29(t,y,v0,v,d,theta)dydx=(v0*t*sin(theta)-y(2)/(d+v0*t*cos(theta)-y(1)+1e-8);dy(1)=v/(1+dydx2)(1/2);dy(2)=v/(1+dydx(-2)(1/2);dy=dy(:);輸入命令clear;close;v0=90;v=450;d=50;theta=pi/4;T=0.1;i=1;for d=54:-1:36 for theta=0:0.1:pi t,y=ode45(equ29,0,T,0;0,v0,v,d,theta); if max(y(:,1)-d-v0*t*cos(theta)>0 ra

42、nge(i,:)=d,theta; i=i+1; break; end endendfigure;plot(range(:,1),range(:,2);xlabel('d');ylabel('theta');運(yùn)行得到臨界曲線,即在該曲線上方的和值,所對應(yīng)的敵艦位置在打擊范圍內(nèi),曲線下方不在打擊范圍內(nèi).圖4-12 離線算法因此d=50km,不在打擊范圍內(nèi).§4.1.4微分方程穩(wěn)定性理論簡介在處理實(shí)際問題時(shí),對于有些微分方程模型我們不僅要得到問題的解,有時(shí)還需要研究解的穩(wěn)定性,即解對初始值的連續(xù)依賴性,如果解在一定范圍內(nèi)是穩(wěn)定的,那么初始條件發(fā)生一些小的擾

43、動(dòng)(如實(shí)驗(yàn)測量誤差等),對問題的解不會(huì)造成影響.另外,還有一些問題我們并不需要求解,而通過解的變化趨勢的研究,并分析一些特殊解的穩(wěn)定性就可以解決問題.本節(jié)僅介紹幾類特殊,但常用的常微分方程穩(wěn)定性分析的基本方法.一、單個(gè)常微分方程的平衡點(diǎn)及穩(wěn)定性若微分方程 (4.1.30)方程右端不顯含自變量,稱為自治方程.代數(shù)方程的實(shí)根稱為方程(4.1.30)的平衡點(diǎn)(或奇點(diǎn)).注意到,平衡點(diǎn)也是方程的解(奇解).如果從一定范圍內(nèi)的初始條件出發(fā),方程(4.1.30)的解都滿足則稱平衡點(diǎn)是穩(wěn)定的.對于一些不易求解的問題,我們可以不求方程(4.1.30)的解,不用定義來判斷平衡點(diǎn)的穩(wěn)定性,下面我們來介紹這種方法.

44、將在點(diǎn)作泰勒(Taylor)展開,只取一次項(xiàng),得到方程(4.1.30)的近似線性方程 (4.1.31)也是方程(4.1.31)的平衡點(diǎn),方程(4.1.31)的通解為.關(guān)于點(diǎn)穩(wěn)定性有如下結(jié)論: 若,則對于方程(4.1.31)和(4.1.30)都是穩(wěn)定的; 若,則對于方程(4.1.31)和(4.1.30)都是不穩(wěn)定的.二、二階常微分方程組的平衡點(diǎn)及穩(wěn)定性現(xiàn)在討論二階微分方程組 (4.1.32)它的解在以為坐標(biāo)的歐氏空間中決定了一條曲線,如果把時(shí)間看作參數(shù),僅考慮為坐標(biāo)的歐氏空間,此空間稱為方程(4.1.32)的相平面(若方程組是高階,則稱為相空間).對于右端函數(shù)不顯含時(shí)間的自治系統(tǒng) (4.1.33

45、)代數(shù)方程組的實(shí)根稱為方程(4.1.33)的平衡點(diǎn),記作.它也是方程(4.1.33)的解.如果從一定的范圍內(nèi)的初始條件出發(fā),方程(4.1.33)的解都滿足則稱平衡點(diǎn)是穩(wěn)定的,否則稱是不穩(wěn)定的.與單個(gè)方程的討論類似,對于不易求解的微分方程組,我們可以通過研究與其對應(yīng)的近似線性方程組,從而得到平衡點(diǎn)和穩(wěn)定性.在這里,我們省略證明過程,只給出判別平衡點(diǎn)是否穩(wěn)定的判別準(zhǔn)則.令關(guān)于點(diǎn)穩(wěn)定性有如下結(jié)論: 當(dāng)且時(shí),平衡點(diǎn)是穩(wěn)定的; 當(dāng)或時(shí),平衡點(diǎn)是不穩(wěn)定的.§4.1.5 動(dòng)物種群模型分析一、單一種群模型(一) Malthus增長模型英國經(jīng)濟(jì)學(xué)家和人口統(tǒng)計(jì)學(xué)家Malthus(1766-1834)根據(jù)

46、百余年的統(tǒng)計(jì)資料, 在1798年提出了聞名于世的人口指數(shù)增長模型(Malthus人口模型).此模型也適用于單種族問題.此模型有兩個(gè)基本假設(shè)前提： 假設(shè)人口數(shù)量x (t )是時(shí)間t的連續(xù)可微函數(shù), 且x (0) = x0. 雖然人口數(shù)量總是離散變化的, 但實(shí)際上增加一個(gè)人所引起的變化與龐大的人口數(shù)量相比是微不足道的. 這里實(shí)際上是將離散變量作連續(xù)化處理, 在差分方程中則恰好相反. 人口數(shù)量的增長速度與現(xiàn)有人口數(shù)量成正比, 比例系數(shù)為r. 這個(gè)假設(shè)在當(dāng)時(shí)是比較合理的.在此基礎(chǔ)上, Malthus提出了如下的人口模型(即指數(shù)增長模型)： = rx, x (0) = x0 .不難求得其解為x (t )

47、 = x0ert, 這個(gè)模型比較準(zhǔn)確地反映了在1700-1961年期間的人口總數(shù).(二)Logistic增長模型因?yàn)閞0, x00, x (t ) = x0ert表明人口數(shù)量x (t ) 將隨時(shí)間t按指數(shù)規(guī)律無限地增長, 對于有限的地球資源來說, 這當(dāng)然是不合理的. 即假設(shè)在自然環(huán)境條件較好, 人口密度較小時(shí)是合理的, 而對于人口密度較大時(shí)它是不合理的, 這是因?yàn)槿祟惿婵臻g及可利用資源(食物、水、空氣)等環(huán)境影響對人口的增長起著阻滯作用. 從這里我們可以看出, 在進(jìn)行數(shù)學(xué)建模時(shí), 我們必須作合理的假設(shè).現(xiàn)在我們對指數(shù)增長模型加以改進(jìn).設(shè)人類生存空間及可利用資源(食物、水、空氣)等環(huán)境因素所能

48、容納的最大人口容量為K (稱為飽和系數(shù)). 人口數(shù)量x的增長速率不僅與現(xiàn)有人口數(shù)量成正比, 而且還與人口尚未實(shí)現(xiàn)的部分(相對最大容量K而言)所占比例(K - x)/K成正比, 比例系數(shù)為固有增長率r. 修改后的模型為不難求得其解為.這就是非常著名的Logistic增長模型(即阻滯增長模型), 它是荷蘭數(shù)學(xué)家、生物學(xué)家Verhulst在1839年首次提出的.因子(K - x)/K的生物學(xué)含義是“剩余空間”或尚未利用的增長機(jī)會(huì).圖4-13描繪了與x之間的關(guān)系, 它表明人口變化率隨著人口數(shù)量x的增加而先增后減, 在x = K / 2處到達(dá)最大值, 在x = K處, = 0. 圖4-14描繪了Logi

49、stic曲線x (t ), 且有, 即x = K是穩(wěn)定的平衡點(diǎn), 從模型本身的意義看這是明顯的結(jié)果. 圖4 -13 圖4 -14 Logistic模型的曲線Logistic模型用途十分廣泛, 除了用于預(yù)測人口增長外, 也可完全類似地用于蟲口增長、疾病的傳播、謠言的傳播、技術(shù)革新的推廣、銷售預(yù)測等.但這個(gè)模型的最大缺點(diǎn)是飽和系數(shù)K不易確定.二、兩種物種相互競爭與相互依存模型當(dāng)某個(gè)自然環(huán)境中只有一種生物的群體(生態(tài)學(xué)上稱為種群)生存時(shí), 人們常用Logistic模型來描述這個(gè)種群數(shù)量的演變過程.即x(t ) 是種群在時(shí)刻t的數(shù)量, r是固有增長率, K是環(huán)境資源容許的種群最大數(shù)量, 且x = K是

50、穩(wěn)定的平衡點(diǎn).如果一個(gè)自然環(huán)境中有兩個(gè)或兩個(gè)以上種群生存, 那么它們之間就要存在著或是相互競爭, 或是相互依存, 或是弱肉強(qiáng)食(捕食-被捕食模型)的關(guān)系. 種群的相互競爭當(dāng)兩個(gè)種群為了爭奪有限的同一種食物來源和生活空間而進(jìn)行生存競爭時(shí), 最常見的結(jié)局是競爭力較弱的種群滅絕, 競爭力較強(qiáng)的種群達(dá)到環(huán)境容許的最大數(shù)量.設(shè)有甲乙兩個(gè)種群, 當(dāng)他們獨(dú)自在一個(gè)自然環(huán)境中生存時(shí), 數(shù)量的演變均遵從Logistic規(guī)律, x1 (t ), x2 (t )是兩個(gè)種群的數(shù)量, r1, r2是它們固有增長率, N1, N2是它們的最大容量. 于是對種群甲有, 其中因子反映由于甲對有限資源的消耗導(dǎo)致對它本身增長的阻

51、滯作用, 可解釋為相對于N1而言單位數(shù)量的甲消耗的供養(yǎng)甲的食物量(設(shè)食物總量為1).當(dāng)兩個(gè)種群在同一自然環(huán)境中生存時(shí), 考察由于乙消耗同一種有限資源對甲的增長產(chǎn)生的影響, 可以合理地在因子中再減去一項(xiàng), 該項(xiàng)與種群乙的數(shù)量x2 (相對于N2而言)成正比, 得到種群甲的增長方程, 這里的s1意義是, 單位數(shù)量乙(相對N2而言)消耗的供給甲的食物量為單位數(shù)量甲(相對N1)消耗的供給甲的食物量的s1倍.類似地, 得到種群乙的增長方程, 對s2可作相應(yīng)地解釋.關(guān)于種群的相互競爭模型的穩(wěn)定性分析見表4-1.表4-2種群相互競爭模型的平衡點(diǎn)及其穩(wěn)定性平衡點(diǎn)pq穩(wěn)定條件P1 (N1, 0)r1 - r2 (1- s2)- r1 r2 (1- s2)s11, s 21N1 (1- s1) N2 (1- s2) 1- s1s2 1- s1s2P2 (0, N2 )r2 - r1 (1- s1)- r1 r2 (1- s1)s11, s 21P3 (, )r1 (1- s1) + r2 (1- s2)1- s1 s2r1 r2 (1- s1) (1- s2)1- s1