1常微分方程數(shù)值解1常微分方程數(shù)值解2實(shí)例1海上緝私

海防某部緝私艇上的雷達(dá)發(fā)現(xiàn)正東方向c海里處有一艘走私船正以速度a向正北方向行駛，緝私艇立即以最大速度b(>a)前往攔截。如果用雷達(dá)進(jìn)行跟蹤時(shí)，可保持緝私艇的速度方向始終指向走私船。

建立任意時(shí)刻緝私艇位置及航線的數(shù)學(xué)模型,并求解;

求出緝私艇追上走私船的時(shí)間。

a北bc艇船2實(shí)例1海上緝私海防某部緝私艇上的雷達(dá)發(fā)現(xiàn)正東方向c海3實(shí)例1海上緝私

建立坐標(biāo)系如圖:t=0艇在(0,0),船在(c,0);船速a,艇速b

時(shí)刻t艇位于P(x,y),船到達(dá)Q(c,at)模型:0yxcR(c,y)Q(c,at)P(x,y)b由方程無法得到x(t),y(t)的解析解需要用數(shù)值解法求解

3實(shí)例1海上緝私建立坐標(biāo)系如圖:t=0艇在(0,4龍格—庫塔方法的MATLAB實(shí)現(xiàn)[t,x]=ode23(@f,ts,x0,opt)3級(jí)2階龍格-庫塔公式[t,x]=ode45(@f,ts,x0,opt)5級(jí)4階龍格-庫塔公式f是待解方程寫成的函數(shù)m文件：functiondx=f(t,x)dx=[f1;f2;;fn];ts=[t0,t1,…,tf]輸出指定時(shí)刻t0,t1,…,tf的函數(shù)值ts=t0:k:tf輸出[t0,tf]內(nèi)等分點(diǎn)處的函數(shù)值x0為函數(shù)初值(n維)輸出t=ts,x為相應(yīng)函數(shù)值(n維)opt為選項(xiàng)，缺省時(shí)精度為：相對(duì)誤差10-3，絕對(duì)誤差10-6,計(jì)算步長按精度要求自動(dòng)調(diào)整4龍格—庫塔方法的MATLAB實(shí)現(xiàn)[t,x]=ode25實(shí)例1海上緝私(續(xù))模型的數(shù)值解0yxc(x(t),y(t))ab設(shè):船速a=20(海里/小時(shí))

艇速b=40(海里/小時(shí))

距離c=15(海里)求:緝私艇的位置x(t),y(t)緝私艇的航線y(x)jisi.m,seajisi.m5實(shí)例1海上緝私(續(xù))模型的數(shù)值解0yxc(x(t),6%Creatthefunctionforjisi%Letx(1)=x,x(2)=yfunctiondx=jisi(t,x,a,b,c)s=sqrt((c-x(1))^2+(a*t-x(2))^2);%dx=[b*(c-x(1))/s;b*(a*t-x(2))/s];dx=[1;1];dx(1)=b*(c-x(1))/s;dx(2)=b*(a*t-x(2))/s;6%Creatthefunctionforjisi7clear,clf,shg%Setthedefiniedtime%ts=0:0.05:0.5;%ts=0:0.1:1.6;n=length(ts);x0=[00];a=35;b=40;c=15;opt1=odeset('RelTol',1e-6,'AbsTol',1e-9);[t,x]=ode45(@jisi,ts,x0,opt1,a,b,c);%a=35;b=40;c=15;%opt1=odeset('RelTol',1e-6,'AbsTol',1e-9);%[t,x]=ode45(@jisi,ts,x0,opt1,a,b,c);%exactsolutionx1=cy1=a*t;%outputt,x(t),y(t)anddrawx(t),y(t)[t,x,y1]plot(t,x),grid,gtext('x(t)','FontSize',16),gtext('y(t)','FontSize',16),pause%drawy(x):thepositionoftatchjisiplot(x(:,1),x(:,2),'r*'),gridxlabel('x','FontSize',16),ylabel('y','FontSize',16)7clear,clf,shg8

實(shí)例1海上緝私(續(xù))模型的數(shù)值解

a=20,b=40,c=15走私船的位置x1(t)=c=15y1(t)=at=20t

t=0.5時(shí)緝私艇追上走私船緝私艇的航線y(x)tx(t)y(t)0000.051.99840.06980.103.98540.29240.155.94450.69060.207.85151.28990.259.67052.11780.3011.34963.20050.3512.81704.55520.4013.98066.17730.4514.74518.02730.5015.00469.9979y1(t)01.02.03.04.05.06.07.08.09.010.08實(shí)例1海上緝私(續(xù))模型的數(shù)值解a=20,b=9實(shí)例1海上緝私(續(xù))模型的數(shù)值解

設(shè)b，c不變，a變大為30,35,接近40,觀察解的變化:a=35,b=40,c=15t=?緝私艇追上走私船tx(t)y(t)y1(t)00000.13.95610.50583.50.27.59282.13087.00.310.52404.828310.50.412.53848.275514.00.513.755112.083017.51.214.998640.016442.01.314.999644.016545.51.415.011748.018349.01.515.002352.014652.51.614.986655.948656.0累積誤差較大提高精度!9實(shí)例1模型的數(shù)值解設(shè)b，c不變，a變大為30,35,10實(shí)例1

海上緝私(續(xù))模型的數(shù)值解

a=35,b=40,c=15opt=odeset('RelTol',1e-6,'AbsTol',1e-9);[t,x]=ode45(@jisi,ts,x0,opt);t=1.6時(shí)緝私艇追上走私船

緝私艇的航線y(x)判斷“追上”的有效方法?tx(t)y(t)y1(t)00000.13.9561040.5058133.50.27.5928222.1306787.00.310.5219214.82930810.50.412.5394548.26984014.00.513.75397412.07534417.51.214.99961640.00000542.01.314.99996344.00000545.51.414.99999348.00000549.01.514.99999852.00000552.51.615.00002055.99993156.010實(shí)例1模型的數(shù)值解a=35,b=40,c=111實(shí)例1海上緝私(續(xù))模型的解析解

11實(shí)例1海上緝私(續(xù))模型的解析解12實(shí)例1海上緝私(續(xù))模型的解析解

緝私艇的航線y(x)的解析解x=c時(shí)

緝私艇追上走私船的y坐標(biāo)

緝私艇追上走私船的時(shí)間:

a=20,b=40,c=15

t1=0.5

a=35,b=40,c=15

t1=1.6

12實(shí)例1海上緝私(續(xù))模型的解析解緝私艇的航線y(13實(shí)例2

弱肉強(qiáng)食問題自然界中同一環(huán)境下兩個(gè)種群之間的生存方式相互競爭相互依存弱肉強(qiáng)食弱肉強(qiáng)食種群甲靠豐富的自然資源生存食餌(Prey)種群乙靠捕食種群甲為生捕食者(Predator)兩個(gè)種群的數(shù)量如何演變?

歷史背景——一次世界大戰(zhàn)期間地中海漁業(yè)的捕撈量下降(食用魚和鯊魚同時(shí)捕撈)，但是其中鯊魚的比例卻增加，為什么？13實(shí)例2弱肉強(qiáng)食問題自然界中同一環(huán)境14實(shí)例2

弱肉強(qiáng)食模型食餌(甲)的密度x(t),捕食者(乙)的密度y(t)甲獨(dú)立生存的增長率r乙使甲的增長率減小,減小量與

y

成正比乙獨(dú)立生存的死亡率d甲使乙的死亡率減小,減小量與

x成正比Volterra模型x(t),y(t)無解析解14實(shí)例2弱肉強(qiáng)食模型食餌(甲)的密15實(shí)例2

弱肉強(qiáng)食模型的數(shù)值解

猜測(cè)x(t),y(t)是周期函數(shù);

y(x)是封閉曲線

數(shù)值積分計(jì)算一個(gè)周期的平均值：shier.m,shier1.m15實(shí)例2弱肉強(qiáng)食模型的數(shù)值解猜測(cè)x(t),y(t16functiony=shier(t,x)r=1;d=0.5;a=0.1;b=0.02;%y=diag([r-a*x(2),-d+b*x(1)])*x;y=[(r-a*x(2))*x(1),(-d+b*x(1))*x(2)]';16functiony=shier(t,x)17ts=0:0.1:15;x0=[25,2];[t,x]=ode45('shier',ts,x0);[t,x], plot(t,x),grid,gtext('x(t)','FontSize',16),gtext('y(t)','FontSize',16),pause,plot(x(:,1),x(:,2)),grid,xlabel('x','FontSize',16),ylabel('y','FontSize',16)%avergevalueofx:Tisabout10.7%trapz(t(1:108),x(1:108,:))/10.717ts=0:0.1:15;18實(shí)例2

弱肉強(qiáng)食模型的解析解

c由初始條件確定

相軌線是封閉曲線(c在一定范圍內(nèi))求x(t),y(t)一周期的平均值:可以證明x(t),y(t)是周期函數(shù)

(周期記作T)18實(shí)例2弱肉強(qiáng)食模型的解析解c由初始條件確定相19實(shí)例2

弱肉強(qiáng)食模型的解析解

x(t),y(t)一周期的平均值:r~食餌增長率a~捕食者對(duì)食餌的捕獲能力d~捕食者死亡率b~食餌對(duì)捕食者的喂養(yǎng)能力結(jié)果解釋與計(jì)算結(jié)果同既相互制約又相互依存19實(shí)例2弱肉強(qiáng)食模型的解析解x(t),y(t)一周20?T2T3T4T1PT1

T2

T3

T4x(t)的“相位”領(lǐng)先y(t)進(jìn)一步分析初值相軌線的方向20?T2T3T4T1PT1T2T321一次大戰(zhàn)期間地中海漁業(yè)的捕撈量下降，但是其中鯊魚的比例卻在增加，為什么？rr-1,dd+1捕撈戰(zhàn)時(shí)捕撈rr-2,dd+2,2<1???xy食餌(魚)減少，捕食者(鯊魚)增加自然環(huán)境

還表明：對(duì)害蟲(食餌)—益蟲(捕食者)系統(tǒng)，使用滅兩種蟲的殺蟲劑,會(huì)使害蟲增加，益蟲減少。21一次大戰(zhàn)期間地中海漁業(yè)的捕撈量下降，但是其中鯊魚的比例卻22食餌-捕食者模型(Volterra)的缺點(diǎn)與改進(jìn)Volterra模型改寫多數(shù)食餌—捕食者系統(tǒng)觀察不到周期震蕩,而是趨向某個(gè)平衡狀態(tài),即存在穩(wěn)定平衡點(diǎn)加