海上缉私模型 论文 数学建模 MATLAB

海上缉私问题

建模 题目二

组别：第五组 组长：练佳翔 组员：邵 * 组员：***

海上缉私问题

摘要

针对海上缉私问题，要求出缉私船是否能追上走私船，或着是求缉私艇追上走私船的位置和时间，就需要知道走私船和缉私艇的位置坐标、大概的行驶路线、及二者的速度。对于走私船和缉私艇的位置坐标，可以由二者的行驶路线 、速度、行驶时间之间的关系得到。而走私船和缉私艇的位置坐标，可用三角函数、坐标关系、圆的位置关系求解。当缉私船追上走私船时，走私船和缉私艇的位置坐标相同，即二者的横坐标相等，纵坐标相等。在此期间，再加以MATLAB软件进行求解。

关键字： 海上缉私 位置坐标 速度 MATLAB软件

问题重述

分别对以下情况建立缉私船的位置和航线的数学模型，自己设定速度等参数，求数值解： (1) 走私船向正东方向非匀速直线行驶，其速度a(t)按正弦规律变化，如图1．已知缉私船以速度b匀速追击，b1.5d（d为常数），两船初始距离c2d．

图1

(2) 两船速度大小都不变，走私船以速度a沿着与正东方向成角的直线行驶，如图2．已知缉私船的速度b1.6a，两船初始距离ca．取25与65，求数值解，并说明走私船按哪个角度逃跑较快？

2

图2

(3) 两船速度大小都不变，走私船以速度a沿半径为r的圆弧向P点逃跑，现有两种方案，如图3．问两种方案是否都能到达P点？已知圆弧半径ra，缉私船的速度b1.4a，两船初始距离

c0.8a．

方案1 方案2

图3

(4)两船速度都大小不变，走私船以速度a先向正东方向直线行驶，一段时间（设尚未被缉私船追上）后改变方向，沿着与正东方向成角(90180)的直线行驶，如图4．已知缉私船的速度b1.2a，两船初始距离c1.5a．取170，求数值解．

图4

(4)

(5) 开始两船速度大小都不变，走私船以速度a向正东方向沿直线行驶，但当两船距离小于r时，缉私船会发现被人追击，将沿正北方向以速度g加速逃跑，如图5．已知r0.5a，g1.5a，缉私船的速度b1.8a，两船初始距离c3a，求数值解．

3

图5

(6) 实际在追击时，缉私船速度方向的改变并不连续，每隔时间t变换一次角度，在两

次变换之间，缉私船按直线运动．若两船速度大小都不变，走私船以速度a向正东方向沿直线行驶，（海里/小时），缉私船的速度b50（海里/小时），两船初始距离c25a30（海里），t60（秒）．试画出缉私船的航线图，建立此时的追击模型，比较与之前模型有何不同，并求数值解．

问题分析

问题一：要确定缉私船追上走私船的位置及时间，就必须确定缉私船、走私船的坐标。

走私船的速度a(t)按正弦规律变化并向正东行驶，因此走私船的位移即向东行驶的距离可以用a(t)与坐标轴围成的面积表示。而缉私艇的坐标设为p(x,y)，再用缉私艇与走私船之间的位置关系，得出x，y的表达式。再根据两者坐标

即可算出结果。 问题二：走私船以速度a沿着与正东方向成角的直线行驶。根据走私船的起始点c与角的位置关系，及走私船的速度，可以算出走私船坐标。而缉私艇的坐标设为

p(x,y)，再用缉私艇与走私船之间的位置关系，得出x，y的表达式。再根据

两者坐标即可算出结果。

问题三：由于走私船的行驶轨迹是圆弧，所以可以利用圆弧与所对圆心角的关系得出走

私船的坐标的值。而缉私艇的坐标设为p(x,y)，再用缉私艇与走私船之间的位置关系，得出x，y的表达式。再根据两者坐标即可算出结果。

4

模型假设

（1） 在整个追赶拦截的过程中，缉私船和走私船都不会有故障发生导致船不能正常行

驶，在海上没有发生因为天气突然发生变化阻碍船前行的情况。 （2） 走私船的行驶速度一直是呈正弦规律变化，缉私船一直是匀速前行。 （3） 建立直角坐标系，在缉私船发现走私船时计时（t=0），设此时走私船的位置在（0，

c），缉私船位置在（0，0）。

（4） 设在任意时刻缉私船的坐标（x，y），走私船到达Q（at，c）点，直线PQ与缉

私船航线相切，切线PQ与y轴正方向夹角为

（5） 设常数d=20，缉私船b=1.5d=30海里/小时，初始距离c=2d=40海里

符号定义

1、任意时刻缉私船的坐标（x，y） 2、走私船到达Q（at，c）点

3、切线PQ与y轴正方向夹角为 4、初始位置c 5、缉私船速度b 6、走私船速度a

7、走私船行驶t时间的路程k

模型建立与求解

模型一：

模型建立：

（1）根据题目所给走私船的速度a(t)的图，可以求解出速度a(t)的表达式：

a(t)0.05dsin(t/5)d

走私船的位移即向东行驶的距离可以用a(t)与坐标轴围成的面积表示如下：

ka(t)d(t)(0.05dsin(t/5)d)d(t)

00ttt0.25dcos()dt0.25d

5又因为走私船从c(0,c)处向正东行驶，因此走私船在t时刻的位置在：

B(a(t)d(t),c)

0t 5

（2）缉私船的速度b1.5d，两船的初始位置c=2d。缉私船在x，y方向是速度分别为

dxbsin1.5dsin dtdydt 得到微分方程：

bcos1.5dcos

sin初始条件为：

kx(kx)(cy)22,cosky(kx)(cy)22

x(0)0，y(0)0

模型求解：

用MATLAB求数值解，可得结果如图所示。

图1.1 x（t），y（t）曲线

6

图1.2 y（x）曲线

模型二

模型建立：

( 1 )走私船以速度a沿着与正东方向成角的直线行驶，时间t时走私船的行驶距离为

ka*t

根据角的三角函数关系可以得出，时间t时走私船的位置：

ACatsint CQatcos

ADAQQXcatsinya

BQCQCBatcosx

ABAD2BQ2(atsinya)(acosx)2

7

（2）缉私船的速度b1.5d，两船的初始位置c=2d。缉私船在x，y方向是速度分别为

dxbsin1.5dsin dtdydt 得到微分方程：

bcos1.5dcos

sin初始条件为：

kx(kx)(cy)22,cosky(kx)(cy)22

x(0)0，y(0)0

（3）把25与65代入（1）中分别得出走私船的坐标

模型求解

当 25时

8

当65时

9

10

模型三

模型建立：

根据我们的观测，方案一的走私船不能够到达p点 又因为辑私船在x，y方向的速度分别为

dydxbsink，bcosk dtdt

既需求sink和cosk

(1) 已知圆弧半径ra，缉私船的速度b1.4a，两船初始距离c0.8a (2)假设走私船和辑私船在t时在在坐标轴上的点为z（m，n） q（x，y） ，通过船经过的圆心角度p来表示点z

P=a*t*360/(2*pi*r) r =a p= 180*t/pi m=a-a*cos(p) n=0.8*a+a*sin(p)

(3)zp连线与水平的夹角设为k，构成一个三角形，需求三边，根据两座标，可得，两

直角边为a-a*cos(p)-x和0.8*a+a*sin(p)-y，勾股定理得第三边为

22（a-a*cos(p)-x）（0.8*aa*sin(p)-y）

（4）编程时记x（1）=x, x(2)=y, x=(x(1),x(2))T,a=20 方案一如下：

缉私船在x,y方向的速度分别为

dxdyb*scos(),b*sin()。 dtdta*tpia*t *2*pi2r4

11

在方案二如下：

缉私船在x,y方向的速度分别为

dxdyb*scos(),b*sin()。 dtdt 12

p点的坐标为(r,rc)即为(30,54)

.3849)处就被抓到了，从两张初值表中可看出，走私船按照方案一逃跑在点(29.7814,50

13

而按照方案二逃跑则是在点(41.1054,51所以按照方案一不可能到达p点，.8688)处被抓。按方案二可以到达p点 模型四

走私船速度与时间的关系a(t)30海里/小时 缉私船以速度b匀速追击 b1.5*d36海里/小时

模型建立

建立判断模型：

此模型用于确定走私船在何时改变方向

b*(a*t-x)dxdt(a*t-x)2(c-y)2 b*(c-y)dy22dt(a*t-x)(c-y)

利用数学软件matlab求得数值解知道（见附件四（1）），若是走私船一直沿着正东

方向逃跑，缉私船需要4.1小时才可追上走私船。因此我们在建立航线模型时假设走私船在t2小时时改变航向 建立航线模型

14

15

缉私船在x,y方向的速度分别为

dxdyb*scos(),b*sin() dtdtdxdtdydtdxdtdydtb*(a*t-x)(a*t-x)2(c-y)2b*(c-y)(a*t-x)2(c-y)2b*(a*t1a*(t-t1)*cos()-x)(a*t1a*(t-t1)*cos()-x)2(a*(t-t1)*sin(w)c-y)2b*(a*(t-t1)*sin()c-y)(a*t1a*(t-t1)*cos()-x)2(a*(t-t1)*sin()c-y)2t20t2

从数值解中，可以得知在该种情况下，缉私船在发现走私船开始追击走私船需要3.9

小时才可以追上走私船。

参考文献

[1] 陈笑缘、张国勇.数学建模.杭州：中国财政经济出版社，2011.1

附件

第一题

function dx=js(t,x) d=20;b=1.5*d;c=2*d; k=-5*cos(t/5)+20*t+5;

s=sqrt((k-x(1))^2+(c-x(2))^2); dx=[b*(k-x(1))/s;b*(c-x(2))/s];

>> ts=0:0.05:2.8; >> x0=[0,0];

16

>> options=odeset('reltol',1e-6,'abstol',1e-9); >> [t,x]=ode45(@js,ts,x0,options); >> d=20;b=1.5*d;c=2*d; >> plot(t,x(:,1),'--',t,x(:,2),'-'),grid

>> ts=0:0.05:2.8; >> x0=[0,0];

>> options=odeset('reltol',1e-6,'abstol',1e-9); >> [t,x]=ode45(@js,ts,x0,options); >> d=20;b=1.5*d;c=2*d; >> zx=-5*cos(t/5)+20*t+5; >> zy=c*ones(length(t),1);

plot(t,x(:,1),'--',t,x(:,2),'-'),grid pause,plot(zx,zy,x(:,1),x(:,2)),grid

0 0 0 0.0500 0.7500 0.0299 0.1000 1.5000 0.1194 0.1500 2.2500 0.2678 0.2000 3.0000 0.4745 0.2500 3.7500 0.7385 0.3000 4.5000 1.0591 0.3500 5.2500 1.4351 0.4000 6.0000 1.8655 0.4500 6.7500 2.3488 0.5000 7.5000 2.8839 0.5500 8.2500 3.4693 0.6000 9.0000 4.1034 0.6500 9.7500 4.7847 0.7000 10.5000 5.5115 0.7500 11.2500 6.2821 0.8000 12.0000 7.0948 0.8500 12.7500 7.9476 0.9000 13.5000 8.8388 0.9500 14.2500 9.7664 1.0000 15.0000 10.7285 1.0500 15.7500 11.7230 1.1000 16.5000 12.7479 1.1500 17.2500 13.8011

17

40.0000 41.2990 42.5981 43.8971 45.1962 46.4952 47.7942 49.0933 50.3923 51.6913 52.9904 54.2894 55.5885 56.8875 58.1865 59.4856 60.7846 62.0836 63.3827 64.6817 65.9808 67.2798 68.5788 69.8779 0 2.3998 4.7980 7.1934 9.5844 11.9698 14.3482 16.7185 19.0796 21.4303 23.7699 26.0974 28.4120 30.7132 33.0005 35.2734 37.5316 39.7749 42.0033 44.2167 46.4155 48.5997 50.7698 52.9264

1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500 1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500 2.0000 18.0000 18.7500 19.5000 20.2500 21.0000 21.7500 22.5000 23.2500 24.0000 24.7500 25.5000 26.2500 27.0000 27.7500 28.5000 29.2500 30.0000 14.8805 15.9839 17.1090 18.2534 19.4149 20.5906 21.7778 22.9733 24.0000 24.7500 25.5000 26.2500 27.0000 27.7500 28.5000 29.2500 30.0001 71.1769 72.4760 73.7750 75.0740 76.3731 77.6721 78.9711 80.2702 81.5692 82.8683 84.1673 85.4663 86.7654 88.0644 89.3634 90.6625 91.9615 55.0699 57.2012 59.3212 61.4307 63.5309 65.6232 77.7090 79.7901 81.5692 82.8682 84.1673 85.4663 86.7654 88.0644 89.3634 90.6624 91.9615

第二题

function dx=js(t,x)

a=30;b=1.6*a;c=a;d=5/36*pi;

s=sqrt((a*t*cos(d)-x(1))^2+(a*t*sin(d)+c-x(2))^2); dx=[b*(a*t*cos(d)-x(1))/s;b*(a*t*sin(d)+c-x(2))/s];

>> ts=0:0.05:2; >> x0=[0,0];

>> options=odeset('reltol',1e-6,'abstol',1e-9); >> [t,x]=ode45(@js,ts,x0,options); >> a=30; >> b=1.6*a; >> c==a;

>> d=5/36*pi; >> zx=a*t*cos(d);

>> zy=(c+a*t*sin(d)).*ones(length(t),1); >> [t,zx,x(:,1),zy,x(:,2)]

18

t

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 1.2 1.25 1.3 1.35 1.4 1.45 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85

0

zt

0 0.0551 0.2228 0.5067 0.9099 1.4348 2.0829 2.855 3.7509 4.769 5.9071 7.1614 8.5273 9.9991

1.3595 2.7189 4.0784 5.4378 6.7973 8.1568 9.5162

30 30.6339 31.2679 31.9018 32.5357 33.1696 33.8036 34.4375 35.0714 35.7053 36.3393 36.9732 37.6071 38.2411 38.875 39.5089 40.1428 40.7768 41.4107 42.0446 42.6785 43.3125 43.9464 44.5803 45.2143 45.8482 46.4821 47.116 47.75 48.3839 49.0178 49.6517 50.2857 50.9196 51.5535 52.1875 52.8214 53.4553

0 2.3992 4.7931 7.176 9.5416 11.8832 14.1937 16.4658 18.692 20.865 22.9777 25.0235 26.9965 28.8918 30.7056 32.4352 34.0793 35.638 37.1128 38.5064 39.8231 41.0683 42.2488 43.3726 44.4492 45.4897 46.4819 47.116 47.7499 48.3838 49.0178 49.6516 50.2856 50.9195 51.5535 52.1873 52.8214 53.4553

10.8757 12.2352 13.5946 14.9541 16.3135 17.673

19.0325 11.5704 20.3919 13.2339 21.7514 14.9819 23.1108 16.8066 24.4703 18.6998 25.8298 20.6535 27.1892 22.6599 28.5487 24.7115 29.9082

26.801

31.2676 28.9216 32.6271 31.0665 33.9865 33.2292 35.346

35.346

36.7055 36.7054 38.0649 38.0649 39.4244 39.4244 40.7839 40.7838 42.1433 42.1433 43.5028 43.5027 44.8622 44.8622 46.2217 46.2216 47.5812 47.5811 48.9406 48.9406 50.3001

50.3

19

1.9 1.95 2

>> plot(t,x(:,1),'--',t,x(:,2),'-'),grid

51.6595 51.6595 53.019

53.019

54.3785 54.3783

54.0892 54.7232 55.3571

54.0892 54.7231 55.3572

pause,plot(zx,zy,x(:,1),x(:,2)),grid

20

第三题

function dx=js(t,x)

a=20;b=1.4*a;c=0.8*a;p=180*t/pi;

s=sqrt((a-a*cos(p)-x(1))^2+(0.8*a+a*sin(p)-x(2))^2) dx=[b*(0.8*a+a*sin(p)-x(2))/s;b*(a-a*cos(p)-x(1))/s];

>> ts=0:0.02:3; >> x0=[0,0];

>> options=odeset('reltol',1e-6,'abstol',1e-9); >> [t,x]=ode45(@js,ts,x0,options);

21