Transcript 导弹跟踪

纯粹数学与应用数学是理解世
界及其发展的一把主要钥匙。
―
里约热内卢宣言
学习Euler,是认识数学的最好途
径。
― 高斯
数学实验
导弹跟踪问题
上海交通大学数学科学学院
实验目的
• 微分方程数学建模
• 回顾微分方程解析解法
• 介绍Euler方法
• 介绍改进Euler方法
• 介绍仿真方法
实际问题
某军的一导弹基地发现正北方向120 km
处海面上有敌艇一艘以90 km/h的速度向正东方
向行驶. 该基地立即发射导弹跟踪追击敌艇,
导弹速度为450 km/h，自动导航系统使导弹在
任一时刻都能对准敌艇.试问导弹在何时何处
击中敌艇？
（追踪问题）
建立坐标系
120
当t =0时,导弹位于原点O,
敌艇位于(0,H)点, H =120 (km).
P(x, y)
当时刻t ,导弹位于P(x(t),y(t))
敌艇位于(90t, H)点
O
x
数学模型
时刻t 导弹的位置为P(x(t),y(t)),其
速度可由水平分速度与垂直分速度合成
2
2
 dx 
 dy 
2


v




w
 dt 
 dt 
导弹方向指向敌艇,故有
dy
Hy

dx ve t  x

dy dx H  y

(
)
dt
dt ve t  x
Vw = 450 (km/h) Ve = 90 (km/h)
导出一阶微分方程组








dx

dt
vw
Hy 2
1 (
)
vet  x
dy

dt
vw
vet  x 2
1 (
)
Hy
x(0)  0
y (0)  0
解析方法
技巧性较强: 消去 t ,化为二阶方程
d 2 x ( H  y)

2
2
d y ( ddyx )  1
连同初始条件
dx
设 p
dy
x y 0
ve
(  )
vw
dx
 0,
dy
0
y 0
降阶得到一阶可分离变量方程
dp  p 2  1

dy
Hy

1
H  Hy 
p  [(
) (
) ]
2 Hy
H
进而解出
1 ( H  y) 1 H  ( H  y)1
H
x [ 

]
2 H (  1)
1 
1  2
导弹击中艇时y=H, 得到此时其位置和时刻
H
L
xL
 25 (km), t  T   0.2778
2
1 
ve
数值方法 (Euler方法)
Euler 方法十分简单, 就是用差商代替微商
即采用
dx
x xk 1  xk dy
y yk 1  yk


,


dt
 t tk 1  tk
dt
t tk 1  tk
通常取Δt为常数τ,就得到由第 k 步的值到
第k+1步的值之间的关系式
伟大的欧拉（Leonhard Euler，1707～1783
）
 数学天才：约翰·伯努利指引
13岁进巴塞尔大学
 著作辉煌：涉及几乎每一数学领域
886本书和论文 ，汇成100巨册

毅力惊人：失明17年，完成书和约400篇论文
 非凡记忆和心算力：前100个素数的前六次幂
 风格高尚：支持年青的拉格朗日
 热爱生活
 瑞士法郎上的欧拉
 邮票上的欧拉
解析法降阶后方程的Euler格式
由方程
dx
p
dy
此例告诉
我们高阶方
程可化为一
阶方程组
dp  p  1

dy
Hy
2
和初始条件 x y 0  0,
p y 0
dx

dy
0
y 0
 Euler 迭代格式
xk 1  xk  hpk ,
x0  0, p0  0
pk 1  pk  h
 1  pk2
H  kh
h=H/n为步长
Matlab
文件m4_1
原模型的 Euler 迭代格式
设 t = tk+1 时导弹位置为(xk+1 ，yk+1)
vw
 xk 1  xk 
H  yk 2
1 (
)

vet k  xk


vw
 yk 1  yk 
vet k  xk 2

1 (
)
H  yk

 x  0, y  0
0
0
计算到 yk－1 <H, yk≥H 停止，取L ≈ xk
利用Matlab
参见教材中介绍的m文件并运行
例如 m4_3,当τ = 0.05
k
1
2
3
4
5
6
则
tk
0.05
0.10
0.15
0.20
0.25
0.30
xk
0.000 00
1.037 36
3.412 05
7.646 15
14.867 90
29.194 80
yk
22.500 00
44.976 07
67.350 41
89.448 43
110.757 96
128.107 02
L≈x6= 29.195 T ≈0.3244
当τ = 0.001
L≈x278=25.049 T ≈0.27833
表1是对应不同的τ ，用Euler法所得相应的步长推
进次数和计算结果
表1
τ
n
L
T
0.1
3
0.05
6
0.005
56
0.0001
278
22.67495 29.19480 25.66731 25.04935
0.25194 0.32439 0.28519 0.27833
关于微分方程数值方法
大多数微分方程－即使形式十分简单，
也可能无法用解析方法求解，因此数值方法
极其重要。
改进Euler法（预测－校正法）
一维例子 方程与初始条件
dx
 f ( x, t ),
dt
Euler 格式
x t 0  x0
xk 1  xk   f ( xk , tk ) (tk  k )
改进Euler 格式
*
x
预测值 k 1  xk   f ( xk , t k )

*
x

x

[
f
(
x
,
t
)

f
(
x
k
k k
k 1 , t k 1 )]
校正值 k 1
2
改进Euler法的几何解释
dx
 f ( x, t ) 的解为x(t) , 代入方程从tk→tk+1
若
dt
f
积分导出
xk 1  xk  
t k 
tk
f =f (x(t),t)
f ( x(t ), t )dt
Euler法
O
xk+1- xk=τf (xk,tk)
改进Euler法
xk+1-
*
f
(
x
xk=τ[f (xk,tk)+
k 1 , t k 1 ) ]/2
tk
tk+1
t
仿真方法 模仿真实事件行为和过程
在这个问题上，就是一步步地模拟导弹追踪敌
舰的实际过程
在t = 0 时,敌艇在M0(0,H),
M0
M1
导弹在原点P0指向M0
在t =τ时,敌艇的位置为
M1(veτ, H)， 导弹的位置为
P1 (x1,y1), 而 x1=0, y1= vwτ,
P1 θ1
导弹飞行方向的倾角为
P0
θ1 =arctan[(H - vwτ)/veτ]
M0 M1
M2
在t =2τ时,敌艇的位置
为M2(2τve, H)， 导弹的
θ2
P2
P1
位置为P2(x2, y2),
x2= x1+ vwτcosθ1
y2= y1+ vwτsinθ1
θ1
P0
由θ1表达式可写出cosθ1和
sinθ1的表达式.
此时导弹飞行方向的倾角为
依此类推…
θ2= arctan[(H-y2)/(2veτ-x2)]
仿真迭代格式
在t=(k+1)τ时，导弹位于Pk+1(xk+1, yk+1)





xk 1  xk  vw cos  k
yk 1  yk  vw sin  k
其中
cos  k 
sin  k 
kve  xk
( kve  xk ) 2  ( H  yk ) 2
H  yk
( kve  xk ) 2  ( H  y k ) 2
使用Matlab的计算结果 （m4_4）
例如τ = 0.05
k
1
2
tk
0.05
0.10
xk
0.000 00
1.037 36
yk
22.500 00
44.976 07
3
4
5
0.15
0.20
0.25
3.412 05
7.646 15
14.867 90
67.350 41
89.448 43
110.757 96
6
0.30
29.194 80
128.107 02
未用微分方程得到同样结果!
实验任务
参见教材
2.在对由模型直接得到的微分方程组用Euler法计
算时，我们计算到yk<H,yk+1≥H即停止，但这样
的做法可能会有不小误差，即使步长减小而结果
仍可能不能改善结果，因此可以采取变步长法：
当计算到yk<H<yk+1，两端距离H较远时，以xk, yk
为起点,缩小步长继续计算，将改善结果.
试用这个方法改进授课中计算的结果.
实验任务
参见教材
3.如果当基地发射导弹的同时,敌艇立即由仪器发
觉. 假定敌艇为一高速快艇，它即刻以135 km/h
的速度与导弹方向垂直的方向逃逸,问导弹何时
何地击中敌舰?
任务5是成其它角度
补充任务
追踪问题有时也以猎犬追兔形式出现，
若猎犬位于一个半径300m的圆形场地中央，
发现兔子在场地边缘，即以 800m/min
速度追逐，同时兔子以 650m/min 的速度
逃逸，问猎犬需要多久才能捕获兔子
（首先尝试建立数学模型，再用数值方法
求解;其次考虑仿真方法）.
谢谢各位！