四人追逐游戏问题（matlab实现），产生的思考。

1四人追逐

如下图,在正方形ABCD的四个顶点各有一个人。设在初始时刻t=0时，四人同时出发匀速以v沿顺时针走向下一个人。如果他们始终对准下一个人为目标行进，最终结果会如何。作出各自的运动轨迹。
首先一一分析四个人的运动轨迹：
分析如下：

按照上述分析过程写的代码如下：

clc; clear;
%模拟追逐运动
n=240;
x=zeros(4,n);
y=zeros(4,n);
dt=0.05;%时间间隔
v=10;
x(1,1)=100;y(1,1)=0;%第一个人的坐标
x(2,1)=0;y(2,1)=0;%第2个人坐标
x(3,1)=0;y(3,1)=100;%3
x(4,1)=100;y(4,1)=100;
for j=1:n%第1个人和第2个人距离d1=sqrt((x(2,j)-x(1,j))^2+(y(2,j)-y(1,j))^2);%第i+1个人和第i个人距离cosx1=(x(1,j)-x(2,j))/d1;sinx1=(y(2,j)-y(1,j))/d1;x(1,j+1)=x(1,j)-v*dt*cosx1;y(1,j+1)=y(1,j)+v*dt*sinx1;%第2个人d2=sqrt((x(3,j)-x(2,j))^2+(y(3,j)-y(2,j))^2);cosx2=(x(3,j)-x(2,j))/d2;sinx2=(y(3,j)-y(2,j))/d2;x(2,j+1)=x(2,j)+v*dt*cosx2;y(2,j+1)=y(2,j)+v*dt*sinx2;%第3个人d3=sqrt((x(4,j)-x(3,j))^2+(y(4,j)-y(3,j))^2);cosx3=(x(4,j)-x(3,j))/d2;sinx3=(y(3,j)-y(4,j))/d2;x(3,j+1)=x(3,j)+v*dt*cosx3;y(3,j+1)=y(3,j)-v*dt*sinx3;%第4个人d4=sqrt((x(1,j)-x(4,j))^2+(y(1,j)-y(4,j))^2);cosx4=(x(4,j)-x(1,j))/d4;sinx4=(y(4,j)-y(1,j))/d4;x(4,j+1)=x(4,j)-v*dt*cosx4;y(4,j+1)=y(4,j)-v*dt*sinx4;end
%画出动态效果图
% for j=1:n
%  plot(x(1,j),y(1,j),'ro',x(2,j),y(2,j),'bo',x(3,j),y(3,j),'go',x(4,j),y(4,j),'yo')
%  hold on
%  pause(0.1)
% end%非动态效果图
plot(x(1,:),y(1,:),'r*',x(2,:),y(2,:),'b*',x(3,:),y(3,:),'g*',x(4,:),y(4,:),'y*')

结果图：

下面是书上的分析过程和代码实现；
代码如下：

clc; clear;
%模拟追逐运动
n=240;
x=zeros(4,n);
y=zeros(4,n);
dt=0.05;%时间间隔
v=10;
x(1,1)=100;y(1,1)=0;%第一个人的坐标
x(2,1)=0;y(2,1)=0;%第2个人坐标
x(3,1)=0;y(3,1)=100;%3
x(4,1)=100;y(4,1)=100;
for j=1:nfor i=1:3d=sqrt((x(i+1,j)-x(i,j))^2+(y(i+1,j)-y(i,j))^2);%第i+1个人和第i个人距离cosx=(x(i+1,j)-x(i,j))/d;sinx=(y(i+1,j)-y(i,j))/d;x(i,j+1)=x(i,j)+v*dt*cosx;y(i,j+1)=y(i,j)+v*dt*sinx;endd=sqrt((x(1,j)-x(4,j))^2+(y(1,j)-y(4,j))^2);cosx=(x(1,j)-x(4,j))/d;sinx=(y(1,j)-y(4,j))/d;x(4,j+1)=x(4,j)+v*dt*cosx;y(4,j+1)=y(4,j)+v*dt*sinx;plot(x(1,j),y(1,j),'ro',x(2,j),y(2,j),'bo',x(3,j),y(3,j),'go',x(4,j),y(4,j),'yo')hold onpause(0.1)
end

结果图：

拓展问题1：1.1三人追逐问题

如下图,在边长为100的等边三角形ABC的三个顶点各有一个人。设在初始时刻t=0时，三人同时出发匀速以v沿顺时针走向下一个人。如果他们始终对准下一个人为目标行进，最终结果会如何。作出各自的运动轨迹。
同样先一个个分析每个人的运动轨迹：
如下分析：
按照上述分析过程写的代码如下：

clc; clear;
%模拟追逐运动
n=240;
x=zeros(3,n);
y=zeros(3,n);
dt=0.05;%时间间隔
v=10;
x(1,1)=100;y(1,1)=0;%第一个人初始的坐标
x(2,1)=0;y(2,1)=0;%第2个人坐标
x(3,1)=50;y(3,1)=50*sqrt(3);%3
for j=1:n%第1个人和第2个人距离d1=sqrt((x(2,j)-x(1,j))^2+(y(2,j)-y(1,j))^2);%第i+1个人和第i个人距离cosx1=(x(1,j)-x(2,j))/d1;sinx1=(y(2,j)-y(1,j))/d1;x(1,j+1)=x(1,j)-v*dt*cosx1;y(1,j+1)=y(1,j)+v*dt*sinx1;%第2个人d2=sqrt((x(3,j)-x(2,j))^2+(y(3,j)-y(2,j))^2);cosx2=(x(3,j)-x(2,j))/d2;sinx2=(y(3,j)-y(2,j))/d2;x(2,j+1)=x(2,j)+v*dt*cosx2;y(2,j+1)=y(2,j)+v*dt*sinx2;%第3个人d3=sqrt((x(3,j)-x(1,j))^2+(y(3,j)-y(1,j))^2);cosx3=(x(1,j)-x(3,j))/d2;sinx3=(y(3,j)-y(1,j))/d2;x(3,j+1)=x(3,j)+v*dt*cosx3;y(3,j+1)=y(3,j)-v*dt*sinx3;end
%画出动态效果图
for j=1:nplot(x(1,j),y(1,j),'ro',x(2,j),y(2,j),'bo',x(3,j),y(3,j),'go')hold onpause(0.1)
end%非动态效果图
% plot(x(1,:),y(1,:),'r*',x(2,:),y(2,:),'b*',x(3,:),y(3,:),'g*')

结果如下：

第二种代码写法如下：

clc; clear;
%模拟追逐运动
n=240;
x=zeros(3,n);
y=zeros(3,n);
dt=0.05;%时间间隔
v=10;
x(1,1)=100;y(1,1)=0;%第一个人的坐标
x(2,1)=0;y(2,1)=0;%第2个人坐标
x(3,1)=50;y(3,1)=50*sqrt(3);%3
for j=1:nfor i=1:2d=sqrt((x(i+1,j)-x(i,j))^2+(y(i+1,j)-y(i,j))^2);%第i+1个人和第i个人距离cosx=(x(i+1,j)-x(i,j))/d;sinx=(y(i+1,j)-y(i,j))/d;x(i,j+1)=x(i,j)+v*dt*cosx;y(i,j+1)=y(i,j)+v*dt*sinx;endd=sqrt((x(1,j)-x(3,j))^2+(y(1,j)-y(3,j))^2);cosx=(x(1,j)-x(3,j))/d;sinx=(y(1,j)-y(3,j))/d;x(3,j+1)=x(3,j)+v*dt*cosx;y(3,j+1)=y(3,j)+v*dt*sinx;plot(x(1,j),y(1,j),'ro',x(2,j),y(2,j),'bo',x(3,j),y(3,j),'go')hold onpause(0.1)
end

结果如下：
思考与总结：
对于5人，6人，还是更多的人，解决思路与上面没有太大的差别，基本上是把比如6人追逐的初始位置坐标确定就基本解决了一半问题了。有时间的话，大家可以尝试写写看。下面再写一个八人追逐的问题。

拓展问题1：1.2 八人追逐问题**

如下图,在边长为100m正方形边上有8个人，A,C,E,G四个人位于正方形的顶点，B,D,F,H四个人位于正方形边长的中点。设在初始时刻t=0时，八人同时出发匀速以v沿顺时针走向下一个人。如果他们始终对准下一个人为目标行进，最终结果会如何。作出各自的运动轨迹。
分析过程类似，这里就直接给出代码了：

clc; clear;
%模拟追逐运动
n=350;
x=zeros(8,n);
y=zeros(8,n);
dt=0.05;%时间间隔
v=10;
x(1,1)=100;y(1,1)=0;%第一个人的坐标
x(2,1)=50;y(2,1)=0;%第2个人坐标
x(3,1)=0;y(3,1)=0;%3
x(4,1)=0;y(4,1)=50;x(5,1)=0;y(5,1)=100;
x(6,1)=50;y(6,1)=100;
x(7,1)=100;y(7,1)=100;
x(8,1)=100;y(8,1)=50;
for j=1:nfor i=1:7d=sqrt((x(i+1,j)-x(i,j))^2+(y(i+1,j)-y(i,j))^2);%第i+1个人和第i个人距离cosx=(x(i+1,j)-x(i,j))/d;sinx=(y(i+1,j)-y(i,j))/d;x(i,j+1)=x(i,j)+v*dt*cosx;y(i,j+1)=y(i,j)+v*dt*sinx;end
%     d5=sqrt((x(5,j)-x(6,j))^2+(y(5,j)-y(6,j))^2);
%     cosx5=(x(6,j)-x(5,j))/d5;
%     sinx5=(y(5,j)-y(6,j))/d5;
%     x(5,j+1)=x(5,j)+v*dt*cosx5;
%     y(5,j+1)=y(5,j)-v*dt*sinx5;d8=sqrt((x(1,j)-x(8,j))^2+(y(1,j)-y(8,j))^2);cosx8=(x(8,j)-x(1,j))/d8;sinx8=(y(8,j)-y(1,j))/d8;x(8,j+1)=x(8,j)-v*dt*cosx8;y(8,j+1)=y(8,j)-v*dt*sinx8;
%     plot(x(1,j),y(1,j),'ro',x(2,j),y(2,j),'bo',x(3,j),y(3,j),'go',x(4,j),y(4,j),'yo', x(5,j),y(5,j),'ko',x(6,j),y(6,j),'co',x(7,j),y(7,j),'mo',x(8,j),y(8,j),'y+')
%     hold on
%     pause(0.1)
end%动态效果图
for j=1:nplot(x(1,j),y(1,j),'ro',x(2,j),y(2,j),'bo',x(3,j),y(3,j),'go',x(4,j),y(4,j),'yo', x(5,j),y(5,j),'ko',x(6,j),y(6,j),'co',x(7,j),y(7,j),'mo',x(8,j),y(8,j),'y+')hold onpause(0.1)
end
%非动态效果图
% plot(x(1,:),y(1,:),'ro',x(2,:),y(2,:),'bo',x(3,:),y(3,:),'go',x(4,:),y(4,:),'yo', x(5,:),y(5,:),'ko',x(6,:),y(6,:),'co',x(7,:),y(7,:),'mo',x(8,:),y(8,:),'y+')

结果如下：

拓展问题2：舰艇追击问题

先写我自己的理解过程吧，后面给出标准解法和代码。
我的分析过程如下：

按照我的分析过程编写的代码如下：

clear;
n=151;
dt=0.01;
u=8;
v=12;
x=zeros(2,n);
y=zeros(2,n);x(1,1)=10;y(1,1)=0;%缉私舰的初始坐标
x(2,1)=0;y(2,1)=0;%走私舰的初始坐标for j=1:n-1d=sqrt( (x(2,j)-x(1,j))^2+(y(2,j)-y(1,j)^2) );
cosx=(x(1,j)-x(2,j))/d;
sinx=(y(2,j)-y(1,j))/d;
x(1,j+1)=x(1,j)-v*dt*cosx;
y(1,j+1)=y(1,j)+v*dt*sinx;x(2,j+1)=0;
y(2,j+1)=y(2,j)+u*dt;% plot(x(1,j),y(1,j),'ro',x(2,j),y(2,j),'bo')
% hold on
% pause(0.1)
end% disp(['所需时间为：' num2str(T)])
plot(x(1,:),y(1,:),'r-',x(2,:),y(2,:),'b-')
T=y(2,151)/u;
disp(['所需时间为：' num2str(T)])

结果如下图：
所需时间为：1.5

下面采用理论求解和计算机模拟两种方法求解：

dt=0.01;
n=151;
d=10;
u=8;
v=12;
T=d*v/(v*v-u*u);%理论时间
x1=zeros(n,1);y1=zeros(n,1);
x2=zeros(n,1);y2=zeros(n,1);
x1(1)=0;y1(1)=0;%走私船开始位置
x2(1)=d;y2(1)=0;%缉私船开始位置%仿真曲线
for j=1:n-1x1(j+1)=0;%走私船横坐标y1(j+1)=(j+1)*dt*u;%走私船纵坐标d1=sqrt((x1(j)-x2(j))^2+(y1(j)-y2(j)^2));cosx=(x1(j)-x2(j))/d1;sinx=(y1(j)-y2(j))/d1;x2(j+1)=x2(j)+v*dt*cosx %缉私船横坐标y2(j+1)=y2(j)+v*dt*sinx%缉私船纵坐标
end
%理论曲线
x=d:-0.01:0;
k=u/v;
y=d/2*((x/d).^(1+k)/(1+k)-(x/d).^(1-k)/(1-k))+d*k/(1-k^2);
plot(x1,y1,'r',x2,y2,'--r',x,y,'b')
disp(['理论时间：',num2str(T)])
disp(['仿真时间：',num2str(y1(151)/u)])

结果如下图
理论时间：1.5
仿真时间：1.51

*

*拓展问题3：导弹追踪问题**

位于坐标原点的A船向位于其正东方20个单位的B船发射导弹，导弹始终对准B船，B船以时速V单位（常数）沿东北方向逃逸。若导弹的速度为3V，导弹的射程是50个单位，画出导弹运行的曲线，导弹是否能在射程内击中B船？
我的分析如下：

按照上述的分析过程写的代码如下：

clear;clc
v=200; % 任意给定B船的速度
dt=0.0001; % 定义时间间隔
n=500;
x=zeros(2,n);
y=zeros(2,n);
x(1,1)=0;y(1,1)=0;%A船的初始坐标
x(2,1)=20;y(2,1)=0;%B船的初始坐标t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;
dd=sqrt((x(2,1)-x(1,1))^2+(y(2,1)-y(1,1))^2); % 导弹与B船的初始距离for j=1:nt=t+dt; % 更新导弹击落B船的时间d=d+3*v*dt; % 更新导弹飞行的距离%B船的新坐标x(2,j+1)=x(2,j)+v*m*dt;y(2,j+1)=y(2,j)+v*m*dt;%更新导弹和B船的距离dd=sqrt( (x(2,j)-x(1,j))^2+(y(2,j)-y(1,j))^2 );cosx=(x(2,j)-x(1,j))/dd;sinx=(y(2,j)-y(1,j))/dd;%A导弹新的坐标x(1,j+1)=x(1,j)+3*v*dt*cosx;y(1,j+1)=y(1,j)+3*v*dt*sinx;if dd<=0.05 && d<50disp(['导弹飞行',num2str(d),'单位后击中B船'])disp(['导弹飞行的时间为',num2str(t*60),'分钟'])breakendendplot(x(1,:),y(1,:),'r-',x(2,:),y(2,:),'b-')

结果如下：
导弹飞行27.84单位后击中B船
导弹飞行的时间为2.784分钟
红色线为导弹运动轨迹，左边多出的线不知道咋回事
哎哎，经过大概一个小时的盘查代码终于让我找到了问题所在，只需要改一个变量n,
n是一个经验值，我一开始设为500，后来看看导弹和B船在运行中走了465次（即有465个坐标），就相遇了，所以从465次后导弹和B船就一直被我初始化为0的原点坐标。如下图：
这是x:第一行是A导弹的横坐标变化，第二行是B船横坐标的变化情况
这是y:第一行是A导弹的纵坐标变化，第二行是B船纵坐标
其实，正好从他们坐标来看，他们相遇（击中）坐标为：
（26.5776 , 6.5620），也可印证下面通过微分方程结果的正确性。
故上面代码只改成n=465，其他不变。
代码另写如下：

clear;clc
v=200; % 任意给定B船的速度
dt=0.0001; % 定义时间间隔
n=465;
x=zeros(2,n);
y=zeros(2,n);
x(1,1)=0;y(1,1)=0;%A船的初始坐标
x(2,1)=20;y(2,1)=0;%B船的初始坐标t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;
dd=sqrt((x(2,1)-x(1,1))^2+(y(2,1)-y(1,1))^2); % 导弹与B船的初始距离for j=1:nt=t+dt; % 更新导弹击落B船的时间d=d+3*v*dt; % 更新导弹飞行的距离%B船的新坐标x(2,j+1)=x(2,j)+v*m*dt;y(2,j+1)=y(2,j)+v*m*dt;%更新导弹和B船的距离dd=sqrt( (x(2,j)-x(1,j))^2+(y(2,j)-y(1,j))^2 );cosx=(x(2,j)-x(1,j))/dd;sinx=(y(2,j)-y(1,j))/dd;%A导弹新的坐标x(1,j+1)=x(1,j)+3*v*dt*cosx;y(1,j+1)=y(1,j)+3*v*dt*sinx;if dd<=0.05 && d<50disp(['导弹飞行',num2str(d),'单位后击中B船'])disp(['导弹飞行的时间为',num2str(t*60),'分钟'])breakendend
plot(x(1,:),y(1,:),'r-',x(2,:),y(2,:),'b-')

结果如下：
导弹飞行27.84单位后击中B船
导弹飞行的时间为2.784分钟

其他人的分析过程如下：

其代码如下：

clear;clc
v=200; % 任意给定B船的速度（后期我们可以再改的）
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;   % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
for i=1:2plot(x(i),y(i),'.k','MarkerSize',1);  % 画出导弹和B船所在的坐标，点的大小为1，颜色为黑色(k)，用小点表示grid on;  % 打开网格线hold on;  % 不关闭图形，继续画图
end
axis([0 30 0 10])  % 固定x轴的范围为0-30  固定y轴的范围为0-10
k = 0;  % 引入一个变量  为了控制画图的速度（因为Matlab中画图的速度超级慢）
while(dd>=0.001)  % 只要两者的距离足够大，就一直循环下去。（两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）t=t+dt; % 更新导弹击落B船的时间d=d+3*v*dt; % 更新导弹飞行的距离x(2)=20+t*v*m;  y(2)=t*v*m;   % 计算新的B船的位置 （注：m=sqrt(2)/2）dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);  % 更新导弹与B船的距离tan_alpha=(y(2)-y(1))/(x(2)-x(1));   % 计算斜率，即tan(α)cos_alpha=sqrt(1/(1+tan_alpha^2));   % 利用公式：sec(α)^2 = (1+tan(α)^2)  计算出cos(α)sin_alpha=sqrt(1-cos_alpha^2);  % 利用公式： sin(α)^2 +cos(α)^2 = 1  计算出sin(α)x(1)=x(1)+3*v*dt*cos_alpha;   y(1)=y(1)+3*v*dt*sin_alpha;   % 计算新的导弹的位置k = k +1 ;  if mod(k,500) == 0   % 每刷新500次时间就画出下一个导弹和B船所在的坐标  mod(m,n)表示求m/n的余数for i=1:2plot(x(i),y(i),'.k','MarkerSize',1);hold on; % 不关闭图形，继续画图endpause(0.001);  % 暂停0.001s后再继续下面的操作endif d>50  % 导弹的有效射程为50个单位disp('导弹没有击中B船');break;  % 退出循环endif d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）disp(['导弹飞行',num2str(d),'个单位后击中B船'])disp(['导弹飞行的时间为',num2str(t*60),'分钟'])end
end

结果如下：
导弹飞行27.8019个单位后击中B船
导弹飞行的时间为2.7802分钟
其实清风老师在计算导弹位置的代码有点啰嗦，不够精炼，主要是这三句

tan_alpha=(y(2)-y(1))/(x(2)-x(1)); cos_alpha=sqrt(1/(1+tan_alpha^2));
sin_alpha=sqrt(1-cos_alpha^2);

其实没必要求正切值，直接用导弹(x(1),y(1))和B船(x(2),y(2))的当时时刻的坐标就可以获得余弦值和正弦值。
为了方便直观理解，我还是画个示意图，来解释吧。
经过上面的分析，代码可以精简成两句：
cos_alpha=(x(2)-x(1))/dd;
sin_alpha=(y(2)-y(1))/dd;
另外把击中的位置坐标打印出来
修改后的代码如下：

clear;clc
v=200; % 任意给定B船的速度
dt=0.0000001; % 定义时间间隔
x=[0,20]; % 定义导弹和B船的横坐标分别为x(1)和x(2)
y=[0,0]; % 定义导弹和B船的纵坐标分别为y(1)和y(2)
t=0; % 初始化导弹击落B船的时间
d=0; % 初始化导弹飞行的距离
m=sqrt(2)/2;   % 将sqrt(2)/2定义为一个常量，使后面看起来很简洁
dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2); % 导弹与B船的距离
for i=1:2plot(x(i),y(i),'.r','MarkerSize',1);  % 画出导弹和B船所在的坐标，点的大小为1，颜色为红色，用小点表示grid on;  % 打开网格线hold on;  % 不关闭图形，继续画图
end
axis([0 30 0 10])  % 固定x轴的范围为0-30  固定y轴的范围为0-10
k = 0;  % 引入一个变量  为了控制画图的速度（因为Matlab中画图的速度超级慢）
while(dd>=0.001)  % （两者距离足够小时表示导弹击中，这里的临界值要结合dt来取，否则可能导致错过交界处的情况）t=t+dt; % 更新导弹击落B船的时间d=d+3*v*dt; % 更新导弹飞行的距离x(2)=20+t*v*m;  y(2)=t*v*m;   % 计算新的B船的位置 （注：m=sqrt(2)/2）dd=sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);  % 更新导弹与B船的距离cos_alpha=(x(2)-x(1))/dd;sin_alpha=(y(2)-y(1))/dd;x(1)=x(1)+3*v*dt*cos_alpha; % 计算新的导弹的位置y(1)=y(1)+3*v*dt*sin_alpha;   k = k +1 ;  if mod(k,500) == 0   % 每刷新500次时间就画出下一个导弹和B船所在的坐标  mod(m,n)表示求m/n的余数for i=1:2plot(x(i),y(i),'.r','MarkerSize',1);hold on; % 不关闭图形，继续画图endpause(0.001);  % 暂停0.001s后再继续下面的操作endif d>50  % 导弹的有效射程为50个单位disp('导弹没有击中B船');break;  % 退出循环endif d<=50 & dd<0.001   % 导弹飞行的距离小于50个单位且导弹和B船的距离小于0.001（表示击中）disp(['导弹飞行',num2str(d),'个单位后击中B船'])disp(['导弹飞行的时间为',num2str(t*60),'分钟'])disp('击中点的坐标：');  disp([x(1),y(1)])  %end
end

结果如下：
导弹飞行的时间为2.7802分钟
击中的坐标为 (26.5523,6.5523)
击中点的坐标：
26.5523 6.5523

关于导弹追踪问题，网上还有采用微分方程求解的，下面是分析过程。
下面是网上给出的代码：

clear; clc
options = odeset('reltol',1e-4,'abstol',1e-8);  %  设定相对误差和绝对误差，这样可以提高微分方程数值解的精度
% 下面的[0,0.1]表示时间t的范围，因为我们导弹速度假设的是200，所以这里的范围给小点
[t,y]=ode45('df5' ,[0,0.1],[0 0], options); % y是因变量，第一列为导弹运行的横坐标，第二列为导弹运行的纵坐标
% [t,y]=ode45('df5' ,[0:0.0001:0.1],[0 0], options); % y是因变量，第一列为导弹运行的横坐标，第二列为导弹运行的纵坐标
x = y(:,1);  y =y(:,2);
plot(x, y,'*','MarkerSize',1)  % 画出导弹的运行轨迹
hold on
plot([20,30],[0,10]) % (20,0) (30,10)画出B船的运行轨迹: x-y-20=0
hold on
% 下面我们找到导弹与B船相撞的点（由于Matlab浮点数计算的原因，距离足够近时即可认为相撞）
n =length(t);  % 找到Matlab计算微分方程的数值解时一共有多少个时间点d = 0;  % 初始化导弹飞行的距离
d = 0;  % 初始化导弹飞行的距离
for i = 1:n % 开始循环dd = abs(x(i) - y(i) - 20) / sqrt(2);  % 利用点到直线 x-y-20=0的距离公式计算导弹和船的距离if dd < 0.001 % 如果这个距离足够小了，我们就认为相撞了，但再此之前别忘了判断导弹是否达到了有效射程for k = 2:id = sqrt((x(k)-x(k-1))^2+(y(k)-y(k-1))^2) + d;  % 以直代曲的思想求曲线的长度，即导弹飞行的距离endif d <= 50 % 导弹的有效射程为50个单位,如果没有达到50单位disp(['导弹飞行',num2str(d),'单位后击中B船'])disp(['导弹飞行的时间为',num2str(t(i)*60 ),'分钟']) % 输出导弹击中B船的时间(转换为分钟)disp('击中点的坐标：');  disp([x(i),y(i)])  % 输出导弹击中B船的坐标plot(x(i),y(i),'r*');text(x(i)+0.5,y(i)+0.1,'击中点')endbreak; % 跳出循环end
end
if d >50 || dd >= 0.001 % 如果射程大于50或导弹与B船的距离四种都没有小于0.001（这个数需要根据实际情况调整）disp('导弹没有击中B船');
endt(i) * 200 * 3   % 更快计算导弹飞行距离的公式：速度*时间
% 得到的结果和我们上面以直代曲的结果很接近

function dy=df5(t,y)v = 200; dy=zeros(2,1);dy(1)=3*v*(20+sqrt(2)/2*v*t-y(1))/sqrt((20+sqrt(2)/2*v*t-y(1))^2+(sqrt(2)/2*v*t-y(2))^2);dy(2)=3*v*(sqrt(2)/2*v*t-y(2))/sqrt((20+sqrt(2)/2*v*t-y(1))^2+(sqrt(2)/2*v*t-y(2))^2);
end

结果如下图所示：
导弹飞行27.7509单位后击中B船
导弹飞行的时间为2.7751分钟
击中点的坐标：
26.5163 6.5166

参考资料：

本文参考了mooc上的西北工业大学的肖老师课件和B站上的清风老师以及网上的部分资料。
注：如有理解错误之处，请指出。