笔者偏微分的一次小实验，利用微分对策解决目标拦截问题。自给参数求解，不足之处望多多指正。

问题描述

导弹从坐标原点以的速度
发射，发射角
，在高度H=100m处有一架飞机从水平居原点l=100m处以
的速度水平向x轴负方向飞来，求导弹能打到高度H时的最小射角以及导弹来拦截飞机的最小射程角。如下图所示。（取重力加速度g=9.8m/s^2,取
）

求解思路

求解导弹速度关于时间的数值解序列点

由物理知识可以建立如下的微分等式：

其中：u1(t)为导弹x方向的关于t的函数；u3(t)为导弹y方向上关于时间t的函数。代入数据得：
对上述微分方程使用4阶的龙格—库塔算法(算法原理如下)：$$\{\begin{array}{l}{y}_{i+1}=y_i+\frac{1}{6}\left(K_1+2K_2+2K_3+K_4\right)\\ K_1=hf\left(x_i,y_i\right)\\ K_2=hf\left(x_i+\frac{h}{2},y_i+\frac{1}{2}{K}_1\right)\\ K_3=hf\left(x_i+\frac{h}{2},y_i+\frac{1}{2}{K}_2\right)\\ K_4=hf\left(x_i+h,y_i+K_3\right)\end{array}$$
求解出u2(t)、u4(t)关于时间t的序列点。

角度求解

利用二分法思想，取初始角度为0，步长取
，利用4阶的龙格—库塔算法得到导弹对应的相应时间点的速度的数值解，通过插值型积分得到不同角度下的导弹的位移情况，并求出导弹在y轴方向刚好最大射高为100m时的角度。具体设计的求解算法程序运行步骤如下：

求解结果与程序

导弹能打到高度H时的最小射角

首先将龙格-库塔算法通过程序实现：

function [x,y]=runge_kutta1(ufunc,y0,h,a,b)
%参数表顺序依次是微分方程组的函数名称，初始值向量，步长，时间起点，时间终点
n=floor((b-a)/h);%求步数
x(1)=a;%时间起点
y(:,1)=y0;%赋初值，可以是向量，但是要注意维数
for ii=1:n

x(ii+1)=x(ii)+h;

k1=ufunc(x(ii),y(:,ii));

k2=ufunc(x(ii)+h/2,y(:,ii)+h*k1/2);

k3=ufunc(x(ii)+h/2,y(:,ii)+h*k2/2);

k4=ufunc(x(ii)+h,y(:,ii)+h*k3);

y(:,ii+1)=y(:,ii)+h*(k1+2*k2+2*k3+k4)/6;
%按照龙格库塔方法进行数值求解
end

代入函数运行尝试：

微分程序

%导弹函数x y方向的速度的方程y(1)是速度v0 y（2）为射角
function dy=Missile(x,y)
dy = zeros(2,1);%始化列向量
%定义y方向
dy(1) = -(y(1)*0.025+9.8);
% dy(2) = y(2);
% %定义x方向
% dy(3) =-y(3)*0.01;
dy(2)=-(y(2)*0.025);

[T1,F1]=runge_kutta1(@Missile,[300*sin(pi/3) 300*cos(pi/3)],0.01,0,2);

%subplot(122)
plot(T1,F1)%自编的龙格库塔函数效果
legend('y方向的速度','x方向的速度')
title('自编的   龙格库塔函数')

运行结果(假设射角为pi/3时的尝试)：

再利用龙格-库塔方法与角度迭代思想程序求取导弹能打到高度H时的最小射角程序为：

%定义速度为300m/s,飞机从高100m的地方水平距离为100米以50m/s的速度向y轴方向飞行
%最大飞行时间是2s，取射角为c，按区间值pi/180长度进行迭代

H=100;
v=300;
%飞机的速度
v1=50;
%H=0;
t=0;
%x轴方向l方向
l=100;
%导弹的水平路程
x=0;
%找到一个初始角
for c=0:pi/180:pi/2
    [T,F1]=runge_kutta1(@Missile,[300*sin(c) 300*cos(c) ],0.01,0,2);
    h=0;
    j=1;
    F11=F1(1,:);
     for i=1:0.01:2
         d=0;
         h=h+0.01*F11(j);
         j=j+1;
         if h<H+0.1&&h>H-0.1
             d=1;
             fprintf('高度正好到h时角度：c=%d',c);
             break;
         else
             fprintf('角度：c=%d,达不到高度\n',c);
         end
     end
     if d==1
         break;
     end
end

通过步长为pi/180迭代尝试求解，当h<H+0.1&&h>H-0.1是视作符合高度条件，迭代停止得到结果，最终的运行结果为：

达高的最小角度为c=5.061455e-01。

拦截飞机的最小射程角

利用二分法思想求取最小射程角，同时画出两者的运动轨迹的实现（）

二分程序

程序运行的两者运动轨迹图结果如下图：

可得最小拦截射程角为1.00，此时导弹拦截到飞机

小结

本次实验主要通过龙格-库塔法与插值型积分的思想对导弹拦截问题进行求解，求解得当射角取得弧度为0.5061时，导弹的最大射高可取得100m;最小拦截射程角为1.00，导弹可以拦截到飞机。本次实验在求解最小拦截射程角时，利用了二分法的思想求解，减少的求解的计算量，但在误差分析与灵敏度分析方面，仍是将来的改进提高处，有待进一步的研究使实验得到更精确的结果。