导弹追踪问题数学建模matlab

数学实验结课论文——导弹跟踪问题班级：信息与计算科学姓名：孔雪婷学号：2012518083实验一导弹跟踪问题一．实验目的本实验主要涉及常微分方程。通过实验复习微分方程的建模和求解，介绍两种求解微分方程的数值方法：Euler法和改进的Euler法，还介绍了仿真方法。二．实际问题某军的一导弹基地发现正北方向120km处海面上有敌艇一艘以90km/h的速度向正东方向行驶。该基地立即发射导弹跟踪追击敌艇，导弹的速度为450km/h，自动导航系统使导弹在任一时刻都能对准敌艇。试问导弹在何时何处击中敌艇？三．数学模型设坐标系如下，取导弹基地为原点0（0,0）。x轴指向正东方，y轴指向正北方。当t=0时，导弹位于O，敌艇位于点（0，H），（H=120（km））设导弹t时刻的位置为P（)(),(tytx），由题意，（4.1）其中)/(450hkmvw。另外在t时刻，敌艇位置应该为),(HtMve，其中ve=90（km/h）。由于导弹轨迹的切线方向必须指向敌艇，即直线PM的方向就是导弹轨迹上点P的切线方向，故有xtyHdxdyve（4.2）)(xyHdtdxdtdyve（4.3）方程（4.3）初值条件想x（0）=0，y（0）=0(4.4)构成了一个关于时间变量t的一阶微分方程组的初值问题。由（4.2）得xtyHdydxve)(两边对t求导得dtdxdtdydydxyHdtdydxvyde)(22即有vydeyHdtdydx22把（4.1）写为12dydxvedtdy代入上式，就得到轨迹方程。这是一个二阶非线性微分方程，加上初值条件，则初值问题dydxxvvdydxydyyweyHdx0022211上式分别为（4.5），（4.6），（4.7）。就是导弹的轨迹的数学模型。四．解析方法方程（4.5）可以降阶。令vvwdydxpe,记，则式（4.5）化为一介可分离变量方程yHdydpp12易得12pyHpC由式（4.7）得HC，从而yHHpp12于是有HyHyHHp21（4.8）于是积分又可以得到CyHHHyHx1111121利用式（4.6）得211HC,于是导弹轨迹方程为21111121HxyHHHyH（4.9）设导弹击中敌艇于B（L，H），以y=H代入（4.9）得vvvvewHHLwe2221（4.10）而导弹击中敌艇的时刻vvvvewHLTwe22（4.11）将数据代入（4.10），（4.11）式，得L=25（km），T0.2778（h）五．数值方法1.Euler方法Euler方法十分简单，就是用差商代替微商，即将dydx代之以hxxkk1，而将dydp代之以hppkk1。设导弹到达（yxkk,）处的时刻为tk，那么得到计算的迭代格式为0,0100211pxyppppxxkkkkkkHkhh上式分别为（4.15）（4.16）（4.17）于是,xnLvxenT使用MATLAB，编辑文件m4_.m：functionm4_1(n)H=120;h=H/n;lamda=90/450;x(1)=0;p(1)=0;y=0:h:H;fori=0:n-1x(i+2)=x(i+1)+h*p(i+1);p(i+2)=p(i+1)+h*(lamda*sqrt(1+p(i+1)^2)/(H-y(i+1)));end[x；p]’L=x(n+1)T=x(n+1)/90输入m4_1（4）得到ans=0000.05001.50000.11675.00250.217411.52540.4221L=11.5254T=0.1281使用MATLAB，建立m4_2.m：functionm4_2(N)k=1;forn=NH=120;h=H/n;lamda=90/450;x0=0;p0=0;fori=0:n-1x1=x0+h*p0;p1=p0+h*(lamda*sqrt(1+p0^2)/(H-i*h));x0=x1;p0=p1;endL(k)=x1;T(k)=x1/90;k=k+1;end[N;L;T]'键入m4_2([4,8,12,24,48,96,120,240])得到ans=4.000011.52540.12818.000015.95370.177312.000017.97320.199724.000020.55080.228348.000022.24940.247296.000023.32860.2592120.000023.58030.2620240.000024.15100.2683由方程（4.1），（4.3）解出dtdy和dtdx的表达式，取时间步长tt，对应kttk时导弹轨迹上点的坐标为yxkk,，则Euler格式为0,011002121yxyHxtvvyyxtvyHvxxkkkekkekwkkwkk上式分别为（4.21）（4.22）（4.23）当计算到HHyykk1,即停止，于是，vxxekkLTLL,1或使用matlab，编辑m文件：functionm4_3(t)H=120;Ve=90;Vw=450;x(1)=0;y(1)=0;T(1)=0;fori=1:10e6M=(Ve*T(i)-x(i))/(H-y(i));x(i+1)=x(i)+Vw*t/sqrt(1+1/M.^2);y(i+1)=y(i)+Vw*t/sqrt(1+M.^2);T(i+1)=t+T(i);ify(i+1)=Hbreak;endend[T;x;y]'L=x(i+1)T=x(i+1)/Vem4_3(0.1)Inm4_3at6ans=0000.1000045.00000.20005.361589.67950.300022.6750131.2155L=22.6750T=0.2519m4_3(0.05)Warning:Dividebyzero.Inm4_3at6ans=0000.0500022.50000.10001.037444.97610.15003.412167.35040.20007.646289.44840.250014.8679110.75800.300029.1948128.1070L=29.1948T=0.32442.改进的Euler方法改进的Euler迭代式0,011111*1*111211*11*1121002*2*2*12*1yxyHxtvvyyxtvyHvxxyHxtvvyyyxtvyHvxxxkkkekkkekkkkkekkkekkwkwkwkkwkk上式分别为（4.24）（4.25）（4.26）（4.27）（4.28）六．仿真方法设导弹和敌艇在初始时刻分别位于.0,0)0,0(Mp00和。此时，导弹指向M0。而在t时，导弹的位置yxp111,，其中vyxw11,0，敌艇的位置则为HvMe,1。这时导弹沿Mp11方向飞行，Mp11的倾角为vvewHarctan1；在2t时导弹的位置为yxp222,。以此方式，当1kt时，导弹的位置为yxpkkk111,，其中kwkkvxxcos1（4.33）kwkkvyysin1（4.34）yHxvkxvkkekkek22cos（4.35）yHxvkykkeHkk22sin（4.36）而敌艇的位置为HkvMek,11。如之前一样，当HHyykk1,时，仿真停止。这时xkL1，veLT七．实验任务1.应用数学软件MATLAB对问题（4.12）~（4.14）进行数值计算，先运用Euler法，与表4.2以及表4.3的数据比较，并以更小的步长计算结果，再用改进的Euler法计算（步长与Euler法相同）用改进的Euler方法取步长为0.1和0.05时计算，用MATLAB验证结果。m.文件functionm4_5(t)H=120;Ve=90;Vw=450;x(1)=0;y(1)=0;T(1)=0;fori=1:10e6M=(Ve*T(i)-x(i))/(H-y(i));x1(i+1)=x(i)+Vw*t/sqrt(1+1/M.^2);y1(i+1)=y(i)+Vw*t/sqrt(1+M.^2);T(i+1)=i*t;x(i+1)=0.5*(x1(i+1)+x(i)+Vw*t/sqrt(1+((H-y1(i+1))/(Ve*T(i+1)-x1(i+1))).^2));y(i+1)=0.5*(y1(i+1)+y(i)+Vw*t/sqrt(1+((Ve*T(i+1)-x1(i+1))/(H-y1(i+1))).^2));ify(i+1)=Hbreak;endend[T;x;y]'L=x(i+1)T=x(i+1)/Ve运行结果：m4_5(0.1)Warning:Dividebyzero.Inm4_5at6ans=0000.10002.680844.83970.200012.575288.28680.300027.0724130.2557L=27.0724T=0.3008m4_5(0.05)Warning:Dividebyzero.Inm4_5at6ans=0000.05000.518722.48800.10002.106044.92190.15005.098267.20210.200010.161089.06910.250019.6565108.98800.300024.2409130.9903L=24.2409T=0.2693得到的结果与书上表4.6和4.7结果一致。