实验二-导弹跟踪问题

实验二导弹跟踪问题成员：葛宇202305020035徐帅202311090001一、实验目的本实验主要涉及常微分方程。通过实验复习微分方程的建模和求解；介绍两种微分方程的数值方法：Euler法和改良的Euler法；还介绍了仿真方法。二、实验方法1、数学模型某军的一导弹基地发现正北方向Hkm处海面上有敌舰一艘以km/h的速度向正东方向行驶。该基地立即发射导弹跟踪追击敌艇，导弹速度为km/h，自动导航系统使导弹在任一时刻都能对准敌艇。试问导弹在何时何处击中敌舰？设坐标系如图3.1所示，取导弹基地为原点O(0,0)，x轴指向正东方,y轴指向正北方。可列出一个关于时间变量t的一阶微分方程组的初值问题2、解析方法解上述初值问题，得导弹轨迹方程为设导弹击中敌艇于B(L,H),以Y=H代入(3.9)式,得而导弹击中敌艇的时刻3、Euler方法(数值方法)Euler方法十分简单,就是用差商代替微商，即通常取Δt为常数τ,就得到由第k步的值到第k+1步的值之间的关系式4、改良的Euler方法以一维情况为例,对问题Euler迭代格式是其中而改良的Euler迭代格式那么是其中5、仿真方法如果建立微分方程很困难，或者微分方程很复杂而较难作出数值处理，常常可以用仿真方法。所谓仿真方法，顾名思义，指的是模仿真实行为和过程的方法。在这个具体问题中，就是一步步地模拟导弹追踪敌艇的实际过程。而计算机仿真，那么是在计算机上通过相应的程序和软件来实现对事件运行的实际过程的模拟。此题仿真格式为仍然可以如前那样,当时,仿真停止;或者事先给定误差界,当时,仿真停止,这时，三、实验任务1.应用数学软件或编制计算程序对问题(3.12)~(3.14)进行数值计算，先运用Euler法，与表3.2以及表3.3的数据比拟，并以更小的步长计算结果；再用改良的Euler法计算〔步长与Euler法相同〕。2．在本实验介绍的计算过程中，我们是计算到即停止，然后取,这样做法可能会有不小的误差。有时甚至会出现整体步长改小而结果却未必能改良的情况。由于Euler法或改良的Euler法的计算格式中每一步值的取得仅仅依赖上一步的值，因此在计算过程中改变步长是可行的，即当计算到而y远大于H时，可缩小步长〔例如为原来的十分之一〕以xy作为新起点继续进行迭代。试用这种变步长方法来改良在任务１中得到的结果。3．如果当基地发射导弹的同时，敌艇立即由仪器觉察。假定敌艇为一高速快艇，它即刻一135km/h的速度与导弹方向垂直的方向逃逸，问导弹何时何地击中快艇？试建立数学模型并求解。4、试对上一问题运用仿真方法进行计算。四、实验程序与结果欧拉法的程序：function[x,y]=euler(h)H=120;Ve=90;Vw=450;x(1)=0;y(1)=0;t(1)=0;k=1;whiley(k)<H&k<1000%K<1000防止死循环k=k+1;x(k)=x(k-1)+Vw*h/sqrt(1+((H-y(k-1))/(Ve*t(k-1)-x(k-1)))^2);y(k)=y(k-1)+Vw*h/sqrt(1+((Ve*t(k-1)-x(k-1))/(H-y(k-1)))^2);t(k)=(k-1)*h;end;改良欧拉法的程序：function[x,y]=euler2(h)H=120;Ve=90;Vw=450;x(1)=0;y(1)=22.5;k=1;whiley(k)<H&k<1000k=k+1;x(k)=x(k-1)+Vw*h*((k-1)*Ve*h-x(k-1))/sqrt(((k-1)*Ve*h-x(k-1))^2+(H-y(k-1))^2);y(k)=y(k-1)+Vw*h*(H-y(k-1))/sqrt(((k-1)*Ve*h-x(k-1))^2+(H-y(k-1))^2);t(k)=(k-1)*h;end;改变步长欧拉法的程序：function[x,y]=euler3(h)H=120;Ve=90;Vw=450;x(1)=0;y(1)=0;t(1)=0;k=1;whiley(k)<H&k<1000k=k+1;x(k)=x(k-1)+Vw*h/sqrt(1+((H-y(k-1))/(Ve*t(k-1)-x(k-1)))^2);y(k)=y(k-1)+Vw*h/sqrt(1+((Ve*t(k-1)-x(k-1))/(H-y(k-1)))^2);t(k)=(k-1)*h;end;ify(k)-H>10x(1)=x(k);y(1)=y(k);t(1)=t(k);whiley(k)<H&k<1000k=k+1;x(k)=x(k-1)+Vw*h/sqrt(1+((H-y(k-1))/(Ve*t(k-1)-x(k-1)))^2);y(k)=y(k-1)+Vw*h/sqrt(1+((Ve*t(k-1)-x(k-1))/(H-y(k-1)))^2);t(k)=(k-1)*h\10;endend;总程序：[x,y]=euler(0.05)%不同的欧拉法，不同的精度[x,y]=euler(0.005)[x,y]=euler2(0.05)[x,y]=euler2(0.005)[x,y]=euler3(0.05)[x,y]=euler3(0.005)H=120;Ve=135;Vw=450;h=0.001;%垂直躲避策略t=0;s=0;k=1;Xw(1)=0;Yw(1)=0;Xe(1)=0;Ye(1)=H;while(Xw(k)-Xe(k))^2+(Yw(k)-Ye(k))^2>0.4Xw(k+1)=Xw(k)+Vw*h*(Xe(k)-Xw(k))/sqrt((Xe(k)-Xw(k))^2+(Ye(k)-Yw(k))^2);Yw(k+1)=Yw(k)+Vw*h/sqrt(1+((Xe(k)-Xw(k))/(Ye(k)-Yw(k)))^2);Xe(k+1)=Xe(k)+Ve*h/sqrt(1+((Xe(k)-Xw(k))/(Ye(k)-Yw(k)))^2);Ye(k+1)=Ye(k)-Ve*h*(Xe(k)-Xw(k))/sqrt((Xe(k)-Xw(k))^2+(Ye(k)-Yw(k))^2);s=s+sqrt((Xe(k+1)-Xe(k))^2+(Ye(k+1)-Ye(k))^2);k=k+1;t=t+h;endsprintf('k=%d,t=%8.4f\n',k-1,t)ans=sprintf('X=%8.4f,Y=%8.4f,T=%8.4f\n',Xe(k),Ye(k),s/Ve)x(1)=0;y(1)=0;a(1)=0;b(1)=120;i=0;f(1)=pi/2;t=0.01;%仿真法fori=1:10000 s(i)=sqrt((x(i)-a(i))^2+(y(i)-b(i))^2); x(i+1)=x(i)+450*((a(i)-x(i))/s(i))*t; y(i+1)=y(i)+450*((b(i)-y(i))/s(i))*t; a(i+1)=a(i)+135*((b(i)-y(i))/s(i))*t; b(i+1)=b(i)-135*((a(i)-x(i))/s(i))*t; ifx(i)<a(i)&&x(i+1)>a(i+1)break;endendsprintf('k=%d,t=%7.4f\n',i-1,i*t)运行结果：x=Columns1through6001.03743.41217.646214.8679Column729.1948y=Columns1through6022.500044.976167.350489.4484110.7580Column7128.1070x=Columns1through6000.00860.02600.05230.0877Columns7through120.13240.18660.25060.32440.40840.5027Columns13through180.60770.72340.85030.98861.13851.3004Columns19through241.47461.66131.86112.07412.30092.5418Columns25through302.79723.06763.35353.65543.97394.3095Columns31through364.66295.03485.42585.83696.26876.7223Columns37through427.19867.69888.22418.77589.35549.9645Columns43through4810.605011.279111.989112.737813.528514.3650Columns49through5415.252016.195117.202318.280519.444920.7132Columns55through5722.115323.709025.6673y=Columns1through602.25004.50006.74998.999811.2495Columns7through1213.499015.748417.997520.246322.494724.7427Columns13through1826.990329.237331.483733.729535.974538.2186Columns19through2440.461942.704144.945247.185149.423751.6607Columns25through3053.896256.129958.361660.591362.818665.0435Columns31through3667.265569.484671.700373.912576.120778.3245Columns37through4280.523582.717284.905087.086389.260491.4264Columns43through4893.583395.729997.865099.9867102.0932104.1819Columns49through54106.2497108.2925110.3049112.2793114.2046116.0630Columns55through57117.8228119.4110120.5190x=01.03743.41217.646214.867929.1948y=22.500044.976167.350489.4484110.7580128.1070x=Columns1through600.01040.03140.06330.10630.1608Columns7through120.22700.30530.39590.49920.61550.7453Columns13through180.88901.04691.21951.40731.61081.8305Columns19through242.06702.32102.59302.88383.19423.5250Columns25through303.87714.25154.64925.07145.51945.9947Columns31through366.49887.03357.60098.20328.84309.5236Columns37through4210.248311.021611.848612.735813.691614.7268Columns43through4715.856917.105118.511020.161922.3864y=Columns1through622.500024.750026.999929.249731.499233.7486Columns7through1235.997638.246240.494442.742044.989047.2353Columns13through1849.480751.725253.968556.210758.451560.6907Columns19through2462.928265.163967.397369.628571.857074.0825Columns25through3076.304878.523480.738082.948085.153087.3522Columns31through3689.545091.730593.907896.075798.2328100.3774Columns37through42102.5075104.6205106.7130108.7807110.8176112.8153Columns43through47114.7609116.6329118.3896119.9184120.2558x=Columns1through6001.03743.41217.646214.8679Column729.1948y=Columns1through6022.500044.976167.350489.4484110.7580Column7128.1070x=Columns1through6000.00860.02600.05230.0877Columns7through120.13240.18660.25060.32440.40840.5027Columns13through180.60770.72340.85030.98861.13851.3004Columns19through241.47461.66131.86112.07412.30092.5418Columns25through302.79723.06763.35353.65543.97394.3095Columns31through364.66295.03485.42585.83696.26876.7223Columns37through427.19867.69888.22418.77589.35549.9645Columns43through4810.605011.279111.989112.737813.528514.3650Columns49through5415.252016.195117.202318.280519.444920.7132Columns55through5722.115323.709025.6673y=Columns1through602.25004.50006.74998.999811.2495Columns7through1213.499015.748417.997520.246322.494724.7427Columns13through1826.990329.237331.483733.729535.974538.2186Columns19through2440.461942.704144.945247.185149.423751.6607Columns25through3053.896256.129958.361660.591362.818665.0435Columns31through3667.265569.484671.700373.912576.120778.3245Columns37through4280.523582.717284.905087.086389.260491.4264Columns43through4893.583395.729997.865099.9867102.0932104.1819Columns49through54106.2497108.2925110.3049112.2793114.2046116.0630Columns55through57117.8228119.4110120.5190ans=k=266,t=0.2660ans=X=33.1005,Y=110.1920,T=0.2660ans=k=26,t=0.2700五、实验总结与扩展1、追击问题：这个模型是一个追击问题，追击问题由易到难有以下几种情况：直线追击，这方面比拟有价值的模型是考虑加速度的，甚至加速度的导数非直线追击，这方面主要考虑追击策略与躲避策略〔如成角度躲避，蛇形曲线等〕，也就是方向调整的问题，对于指向性的追击策略，从模型上来说只要追击速度大于逃逸速度在两者方向上的分量，那么就可以击中。解决问题的关键是考虑二者之间的距离关于时间的函数，以及它的导数，以下讨论2、实际问题：对于实际问题来说，还要考虑舰艇的大小，这个因素对最后的结果至关重要，因为导弹调整方向的能力有限。3、追击问题的求解与求“何地〞击中不同，求“何时〞击中，是不需要采用任何解析或数值方法的。易证，定理：两物体之间的距离d(t)关于t的导数是二者速