緝私艇追走私船模型實驗課件

1、2、了解微分方程數(shù)值解思想，掌握兩種簡單的微分方程數(shù)值解方法。3、學會根據(jù)實際問題建立簡單微分方程數(shù)學模型。實驗目的1、學會用MATLAB軟件求解微分方程的初值問題.4、了解計算機數(shù)據(jù)仿真、數(shù)據(jù)模擬的基本方法。 實驗十二 緝私艇追趕走私船模型實驗用MATLAB軟件求微分方程解析解、數(shù)值解理論. 編程求解微分方程數(shù)值解. 微分方程模型實驗:緝私艇追趕走私船. 實驗內(nèi)容1. 微分方程的解一階常微分方程的初值問題: 用matlab軟件求微分方程的解析解dsolve(equation,condition) %求方程equation在初始條件condition下的解，自變量默認為tdsolve(equa

2、tion) %求方程equation的通解%一階導數(shù)用Dy表示，二階導數(shù)用D2y表示， A=dsolve(Dy=5) B=dsolve(Dy=x,x) %求方程的通解，指定自變 量為x C=dsolve(D2y=1+Dy) D=dsolve(D2y=1+Dy,y(0)=1,Dy(0)=0) x,y=dsolve(Dx=y+x, Dy=2*x,x(0)=0, y(0)=1) %解微分方程組例1 (1)dsolve(D2y=exp(x),x)dsolve(D2y=exp(x)dsolve(D2y=exp(x),y(0)=1,Dy(0)=2,x) 微分方程數(shù)值解的計算 歐拉方法、歐拉兩步法、改進的歐

3、拉方法歐拉方法：是一種求解給定初值的常微分方程的一階數(shù)值解方法。對于方程設y=y(x)是它的解，則在解y=y(x)的曲線上每一點處的切線的斜率等于f(x,y)在該點處的值。所以，過初始點A0(x0,y0)，做切線與直線x=x1交點的縱坐標y1近似代替y=y(x)在x=x1處的函數(shù)值y(x1)，即然后再過點A1(x1,y1)，以f(x1,y1)為斜率做切線，與直線x=x2交點的縱坐標y2近似代替y=y(x)在x=x2處的函數(shù)值y(x2)，即A0 x0 x1x2A1A2這樣依次類推得到上述微分方程的解在x=xn+1處函數(shù)值的近似計算公式，待求的曲線為藍色，它的折線近似為紅色A0 x0 x1x2x3

4、A1A2A3A4x4特別地，如果節(jié)點之間是等分的，則這就是著名的歐拉(Euler)公式。事實上，上述公式是用函數(shù)y(x)在xn處的向前差商代替函數(shù)在這點處的導數(shù)，從而將微分方程離散化。同樣，如果用函數(shù)y(x)在xn+1處的向后差商代替導數(shù)值，就得后退的歐拉公式 注意：這兩個公式雖然具有相同的精度，但也是有區(qū)別的，前者是計算yn+1的顯示公式，后者是隱式公式。A0 x0 x1x2x3A1A2A3A4x4將上兩式算術平均即得梯形公式：梯形公式消去了歐拉公式和后退的歐拉公式誤差的主要部分，因此比它們具有更高的精度。這是一個隱式公式，可采用迭代法求解。通常先用歐拉公式提供一個yn附近的初始迭代值，然后

5、再用梯形公式做迭代。2）歐拉兩步法用函數(shù)y(x)在xn處的中心差商代替微分方程中的導數(shù)，這就是歐拉兩步法公式。它比歐拉公式、后退的歐拉公式具有更高的精度。3）改進的歐拉公式 歐拉公式雖然簡單，但是誤差較大，并且隨著迭代次數(shù)的增加，誤差越來越大。梯形公式雖然提高了精度，但因其是隱式的，所以計算麻煩，工作量較大。通常將這兩種方法結合，得到改進的歐拉公式，具體如下：步驟1 先由歐拉公式得到y(tǒng)(xn+1)的一個初始近似值，稱為預測值步驟2 對預測值再用梯形公式校正一次得yn+1，稱為校正值整個過程（稱預測-校正系統(tǒng)）可以表示成這就是改進的歐拉公式。當然也可以將歐拉兩步法與梯形方法結合。再解釋：歐拉公式

6、用函數(shù)y(x)在區(qū)間xn,xn+1上的平均斜率代替了曲線在點xn的斜率。得到啟發(fā)：只要對曲線在區(qū)間xn,xn+1上的平均斜率提供一種算法，就可以得到一種計算yn+1的公式。如果設法在區(qū)間內(nèi)多預測幾個點的斜率值，然后將這些點處斜率值的加權平均值作為曲線在區(qū)間上的平均值，由此構造出由yn計算yn+1的精度更高的計算公式，這就是龍格-庫塔方法的基本思想。歐拉兩步法用區(qū)間xn,xn+1和區(qū)間xn-1,xn上的平均斜率的算術平均值代替了曲線在xn的斜率。編程計算微分方程數(shù)值解例2 已知初值問題（1）用MATLAB軟件求解析解，算出解析解在結點x=pi:0.1:2*pi處的縱坐標；（2）取步長0.1，用歐

7、拉公式求數(shù)值解；（3）取步長0.1，用改進的歐拉公式求數(shù)值解；（4）比較兩種數(shù)值解和解析解，并畫出數(shù)值解和解析解曲線。解：（1）直接在命令窗執(zhí)行命令dsolve(Dy=3/x*y+x3*(exp(x)+cos(x)-2*x, y(pi)=(exp(pi)+2/pi)*pi3,x)ans= x3*exp(x)+x3*sin(x)+2*x2clc;clf;szy_eu=;y=(exp(pi)+2/pi)*pi3;for x=pi:0.1:2*piy=y+0.1*(3/x*y+x3*(exp(x)+cos(x)-2*x);szy_eu=szy_eu,y;endszy_eu t=pi:0.1:2*pi

8、;f=t.3.*exp(t)+t.3.*sin(t)+2*t.2;plot(t,f,b*-) hold onplot(t,szy_eu,r-,linewidth,2)(2) 取步長0.1，用歐拉公式求數(shù)值解x=pi:0.1:2*pi;h=0.1; y=(exp(pi)+2/pi)*pi3;szy_imeu=y; for i=1:length(x)-1 jzj=y+h*(3/x(i)*y+x(i)3*(exp(x(i)+cos(x(i)-2*x(i); yy=y+h/2*(3/x(i)*y+x(i)3*(exp(x(i)+cos(x(i)- 2*x(i)+3/x(i+1)*jzj+x(i+1)3

9、*(exp(x(i+1)+cos(x(i+1)-2*x(i+1); szy_imeu=szy_imeu,yy; y=yy; end szy_imeu plot(x,szy_imeu,r-,linewidth,2) f=x.3.*exp(x)+x.3.*sin(x)+2*x.2; hold on plot(x,f,b*-)(3) 取步長0.1，用改進的歐拉公式求數(shù)值解 2. 用MATLAB軟件求微分方程數(shù)值解 MATLAB軟件求常微分方程數(shù)值解的指令有 ode23，ode23t和ode45，分別表示二三階龍格庫塔法，二三階龍格庫塔梯形法和四五階龍格庫塔法。還有其他命令ode113，ode15i，

10、ode15s，ode23s等，讀者可通過help命令查詢。 xb,yb=ode23(fname,xs,xe,sv) 返回結點處的橫坐標xb和縱坐標yb。ode23(fname,xs,xe,sv) 其中用單引號引起的是存儲微分方程的函數(shù)文件名，xs表示自變量的初始值，xe表示自變量的終止值，sv表示迭代初始向量值。該命令執(zhí)行后自動畫出微分方程數(shù)值解曲線。最簡單的使用格式為注意：（1）命令ode求解的是形如y=f(t,y)的微分方程，我們稱它為一階導數(shù)可解出的微分方程。而對于一階導數(shù)解不出的形如f(t,y,y)=0的微分方程，可以用命令ode15i求解，有興趣的讀者可以用help命令獲得。（2）求

11、解微分問題時必須為其指定初值，即上面的sv。 （3）ode只能直接求解一階微分方程，高階微分方程必須等價地轉(zhuǎn)化成一階微分方程組，才能用ode命令求解。比如 y=f(t,y,y,y)， 設y1=y, y2= y, y3= y ， 那么如上的三階微分方程就可以表示為如下的三個一階方程組了。 y1=y2 y2=y3 y3=f(t,y1,y2,y3)例3使用格式: fc.mfunction f=fc(x,y)f=2*y+x+2;在MATLAB命令窗中執(zhí)行命令ode23(fc,0,1,1)例4 求解微分方程建立m文件 erjie.mfunction f=erjie(t,y)f=y(2); 1000*(1

12、-y(1)2)*y(2)-y(1);在MATLAB命令窗中執(zhí)行命令t,y=ode15s(erjie,0,3000,2,0)plot(t,y(:,1),-)4. 模型實驗: 緝私艇追趕走私船 海上邊防緝私艇發(fā)現(xiàn)距c公里處有一走私船正以勻速a沿直線行駛,緝私艇立即以最大速度b追趕,在雷達的引導下，緝私艇的方向始終指向走私船。問緝私艇何時追趕上走私船?并求出緝私艇追趕的路線方程。xyco解 (1) 建立模型走私船初始位在點(0,0)，方向為y軸正方向，緝私艇的初始位在點(c,0)，在時刻t，走私船的位置到達點R(0,at)，緝私艇的位置到達D(x,y)可得緝私艇追擊走私船過程的數(shù)學模型：(2) 模型

13、求解求解析解令：表示追趕時間。表示走私船跑過的距離。a=0.4;c=3;for b=0.6:0.3:1.2r=a/b;x=0:0.05:c;y=c/2*(1/(r+1)*(x/c).(1+r)-1/(1-r)*(x/c).(1-r)+c*r/(1-r2);plot(x,y,r-,linewidth,2)pausehold onendaxis(0 4 0 3.5)gtext(b=0.6)gtext(b=0.9)gtext(b=1.2)取c=3千米，a=0.4千米/分鐘，分別取b=0.6,0.9,1.2千米/分鐘時，編程繪制緝私艇追趕路線圖形.定義m文件函數(shù)：function y=zx(t,y)y

14、=0.5*(t/3)0.5-(3/t)0.5); ode23(zx,3,0.0005,0) (3)用MATLAB軟件求數(shù)值解設c=3千米，a=0.4千米/分鐘，b=0.8千米/分鐘，r=a/b=0.5,用二三階龍格-庫塔算法計算數(shù)值解:(4) 動態(tài)仿真 當建立動態(tài)系統(tǒng)的微分方程模型很困難時，可以用計算機仿真法對系統(tǒng)進行分析研究。所謂計算機仿真就是利用計算機對實際動態(tài)系統(tǒng)的結構和行為進行編程、模擬和計算，以此來預測系統(tǒng)的行為效果。走私船初始位在點(0,0)，方向為y軸正方向，緝私艇的初始位在點(c,0)，走私船的位置到達點緝私艇的位置到達點 追趕方向可用方向余弦表示為： 時間步長緝私艇的位置：

15、則仿真算法： 第二步：計算動點緝私艇D在時刻 時的坐標 計算走私船R在時刻 時的坐標 第一步：設置時間步長, 速度a, b及初始位置仿真算法步驟：第三步：計算緝私艇與走私船這兩個動點之間的距離： 根據(jù)事先給定的距離，判斷緝私艇是否已經(jīng)追上了走私船，從而判斷退出循環(huán)還是讓時間產(chǎn)生一個步長，返回到第二步繼續(xù)進入下一次循環(huán)； 第四步：當從上述循環(huán)退出后，由點列 和可分別繪制成兩條曲線即為緝私艇和走私船走過的軌跡曲線。 取c=3千米，a=0.4千米/分鐘，b=0.8千米/分鐘，r=a/b=0.5顯示船與艇行進路線程序clc;clear;clf;c=3; a=0.4/60; b=0.8/60;d=0.0

16、1;dt=2;t=0;jstx=c;jsty=0;zscx=0;zscy=0;while (sqrt(jstx-zscx)2+(jsty-zscy)2)d) plot(jstx,jsty,ro,markersize,15); hold on plot(zscx,zscy,b*,markersize,15); pause(0.01) t=t+dt; jstx=jstx-b*dt*jstx/sqrt(jstx2+(a*t-jsty)2); jsty=jsty+b*dt*(a*t-jsty)/sqrt(jstx2+(a*t-jsty)2); zscy=a*t;endfprintf(追上所需時間t=%.0fn,t);fprintf(緝私艇的位置=(%.5f,%.5f)n,jstx,jsty);fprintf(走私船的位置=(%.5f,%.5f)n,zscx,zscy);結果分析計算機仿真法計算的結果依賴于時間迭代步長的選取和程序終止條件的設定，修改終止條件的設定和減小時間迭代步長可以提高計算精度，減小誤差。