緝私艇追走私船模型實驗解析

MathematicsLaboratory阮小娥博士ExperimentsinMathematics數(shù)學(xué)實驗2、了解微分方程數(shù)值解思想，掌握兩種簡單的微分方程數(shù)值解方法。3、學(xué)會根據(jù)實際問題建立簡單微分方程數(shù)學(xué)模型。實驗?zāi)康?、學(xué)會用MATLAB軟件求解微分方程的初值問題.4、了解計算機(jī)數(shù)據(jù)仿真、數(shù)據(jù)模擬的基本方法。實驗十二緝私艇追趕走私船模型實驗微分方程數(shù)值解(理論).用MATLAB軟件求微分方程解析解、數(shù)值解.微分方程模型實驗:緝私艇追趕走私船.問題的進(jìn)一步擴(kuò)展.實驗內(nèi)容1.微分方程數(shù)值解一階常微分方程的初值問題:數(shù)值解方法:歐拉方法、歐拉兩步法、改進(jìn)的歐拉方法歐拉方法：是一種求解給定初值的常微分方程的一階數(shù)值解方法。對于方程設(shè)y=y(x)是它的解，則在解y=y(x)的曲線上每一點處的切線的斜率等于f(x,y)在該點處的值。所以，過初始點A0(x0,y0)，做切線與直線x=x1交點的縱坐標(biāo)y1近似代替y=y(x)在x=x1處的函數(shù)值y(x1)，即然后再過點A1(x1,y1)，以f(x1,y1)為斜率做切線，與直線x=x2交點的縱坐標(biāo)y2近似代替y=y(x)在x=x2處的函數(shù)值y(x2)，即A0x0x1x2A1A2這樣依次類推得到上述微分方程的解在x=xn+1處函數(shù)值的近似計算公式，待求的曲線為藍(lán)色，它的折線近似為紅色A0x0x1x2x3A1A2A3A4x4特別地，如果節(jié)點之間是等分的，則有這就是著名的歐拉(Euler)公式。事實上，上述公式是用函數(shù)y(x)在xn處的向前差商代替函數(shù)在這點處的導(dǎo)數(shù)，從而將微分方程離散化。同樣，如果用函數(shù)y(x)在xn+1處的向后差商代替導(dǎo)數(shù)值，就得后退的歐拉公式注意：這兩個公式雖然具有相同的精度，但也是有區(qū)別的，前者是計算yn+1的顯示公式，后者是隱式公式。將上兩式算術(shù)平均即得梯形公式：梯形公式消去了歐拉公式和后退的歐拉公式誤差的主要部分，因此比它們具有更高的精度。這是一個隱式公式，可采用迭代法求解。通常先用歐拉公式提供一個yn附近的初始迭代值，然后再用梯形公式做迭代。2）歐拉兩步法用函數(shù)y(x)在xn處的中心差商代替微分方程中的導(dǎo)數(shù)，這就是歐拉兩步法公式。它比歐拉公式、后退的歐拉公式具有更高的精度。3）改進(jìn)的歐拉公式歐拉公式雖然簡單，但是誤差較大，梯形公式雖然提高了精度，但因其是隱式的，所以計算麻煩，工作量較大。通常將這兩種方法結(jié)合，得到改進(jìn)的歐拉公式，具體如下：步驟1先由歐拉公式得到y(tǒng)(xn+1)的一個初始近似值，稱為預(yù)測值步驟2對預(yù)測值再用梯形公式校正一次得yn+1，稱為校正值整個過程（稱預(yù)測-校正系統(tǒng)）可以表示成這就是改進(jìn)的歐拉公式。當(dāng)然也可以將歐拉兩步法與梯形方法結(jié)合。再解釋：則用函數(shù)y(x)在區(qū)間[xn,xn+1]上的平均斜率代替了曲線在點xn的斜率。得到啟發(fā)：只要對曲線在區(qū)間[xn,xn+1]上的平均斜率提供一種算法，就可以得到一種計算yn+1的公式。如果設(shè)法在區(qū)間內(nèi)多預(yù)測幾個點的斜率值，然后將這些點處斜率值的加權(quán)平均值作為曲線在區(qū)間上的平均值，由此構(gòu)造出由yn計算yn+1的精度更高的計算公式，這就是龍格-庫塔方法的基本思想。2.用MATLAB軟件求微分方程的解析解例1dsolve('D2y=exp(x)','x')ans=exp(x)+C1*x+C2dsolve('D2y=exp(x)')ans=1/2*exp(x)*t^2+C1*t+C2dsolve('D2y=exp(x)','y(0)=1,Dy(0)=2','x')

ans=exp(x)+x3.編程計算微分方程數(shù)值解例2（1）直接在命令窗執(zhí)行命令dsolve('Dy=3/x*y+x^3*(exp(x)+cos(x))-2*x','y(pi)=(exp(pi)+2/pi)*pi^3','x')ans=x^3*exp(x)+x^3*sin(x)+2*x^2szy_eu=[];y=(exp(pi)+2/pi)*pi^3;forx=pi:0.1:2*piy=y+0.1*(3/x*y+x^3*(exp(x)+cos(x))-2*x);szy_eu=[szy_eu,y]；endszy_eut=pi:0.1:2*pi;f=t.^3.*exp(t)+t.^3.*sin(t)+2*t.^2;plot(t,szy_eu,'r-','linewidth',2)holdonplot(t,f,'b*--')(2)取步長0.1，用歐拉公式求數(shù)值解

zhao1201.mx=pi:0.1:2*pi;h=0.1;y=(exp(pi)+2/pi)*pi^3;szy_imeu=[y];fori=1:length(x)-1jzj=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*(exp(x(i+1))+cos(x(i+1)))-2*x(i+1));szy_imeu=[szy_imeu,yy];y=yy;endszy_imeuplot(x,szy_imeu,'r-','linewidth',2)f=x.^3.*exp(x)+x.^3.*sin(x)+2*x.^2;holdonplot(t,f,'b*--')(3)取步長0.1，用改進(jìn)的歐拉公式求數(shù)值解

zhao12024.用MATLAB軟件求微分方程數(shù)值解MATLAB軟件求常微分方程數(shù)值解的指令有ode23，de23t和ode45，分別表示二三階龍格—庫塔法，二三階龍格—庫塔梯形法和四五階龍格—庫塔法。最簡單的使用格式為ode23('fname',[xs,xe],sv)其中用單引號引起的是存儲微分方程的函數(shù)文件名，xs表示自變量的初始值，xe表示自變量的終止值，sv表示迭代初始向量值。該命令執(zhí)行后自動畫出微分方程數(shù)值解曲線。例3使用格式:fc.mfunctionf=fc(x,y)f=2*y+x+2在MATLAB命令窗中執(zhí)行命令ode23(‘fc’,[0,1],1)5.模型實驗:緝私艇追趕走私船海上邊防緝私艇發(fā)現(xiàn)距c公里處有一走私船正以勻速a沿直線行駛,緝私艇立即以最大速度b追趕,在雷達(dá)的引導(dǎo)下，緝私艇的方向始終指向走私船。問緝私艇何時追趕上走私船?并求出緝私艇追趕的路線方程。xyco(1)建立模型走私船初始位在點(0,0)，方向為y軸正方向，緝私艇的初始位在點(c,0)，在時刻t，走私船的位置到達(dá)點R(0,at)，緝私艇的位置到達(dá)D(x,y)可得緝私艇追擊走私船過程的數(shù)學(xué)模型：(2)模型求解求解析解令：--zhao1203.m表示追趕時間。表示走私船跑過的距離。取c=3千米，a=0.4千米/分鐘，分別取b=0.6,0.9,1.2千米/分鐘時，緝私艇追趕路線的圖形為定義m文件函數(shù)：functiony=zx(t,y)y=0.5*((t/3)^0.5-(3/t)^0.5)ode23('zx',3,0.0005,0)

用MATLAB軟件求數(shù)值解設(shè)c=3千米，a=0.4千米/分鐘，b=0.8千米/分鐘，r=a/b=0.5,用二三階龍格-庫塔算法計算數(shù)值解:6動態(tài)仿真當(dāng)建立動態(tài)系統(tǒng)的微分方程模型很困難時，可以用計算機(jī)仿真法對系統(tǒng)進(jìn)行分析研究。所謂計算機(jī)仿真就是利用計算機(jī)對實際動態(tài)系統(tǒng)的結(jié)構(gòu)和行為進(jìn)行編程、模擬和計算，以此來預(yù)測系統(tǒng)的行為效果。走私船初始位在點(0,0)，方向為y軸正方向，緝私艇的初始位在點(c,0)，走私船的位置到達(dá)點緝私艇的位置到達(dá)點

追趕方向可用方向余弦表示為：時間步長緝私艇的位置：則仿真算法：

第二步：計算動點緝私艇D在時刻

時的坐標(biāo)

計算走私船R在時刻

時的坐標(biāo)第一步：設(shè)置時間步長,速度a,b及初始位置仿真算法步驟：第三步：計算緝私艇與走私船這兩個動點之間的距離：

根據(jù)事先給定的距離，判斷緝私艇是否已經(jīng)追上了走私船，從而判斷退出循環(huán)還是讓時間產(chǎn)生一個步長，返回到第二步繼續(xù)進(jìn)入下一次循環(huán)；

第四步：當(dāng)從上述循環(huán)退出后，由點列

和可分別繪制成兩條曲線即為緝私艇和走私船走過的軌跡曲線。

取c=3千米，a=0.4千米/分鐘，b=0.8千米/分鐘，r=a/b=0.5顯示船與艇行進(jìn)路線程序clc;clear;c=3;a=0.4/60;b=0.8/60;d=0.01;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);holdonplot(zscx,zscy,'b*','markersize',15);pause(0.01)t=t+dt;jstx=jstx-b*dt*jstx/sqrt(jstx^2+(a*t-jsty)^2);jsty=jsty+b*dt*(a*t-jsty)/sqrt(jstx^2+(a*t-jsty)^2);zscy=a*t;endfprintf('追上所需時間t=%.0f\n',t);fprintf('緝私艇的位置=(%.5f,%.5f)\n',jstx,jsty);fprintf('走私船的位置=(%.5f,%.5f)\n',zscx,zscy);zhao1204.m結(jié)果分析計算機(jī)仿真法計算的結(jié)果依賴于時間迭代步長的選取和程序終止條件的設(shè)定，修改終止條件的設(shè)定和減小時間迭代步長可以提高計算精度，減小誤差。

上機(jī)任務(wù)李繼成（書）