233数学建模

1——（1）

竖直：

根据前面的求解过程，可以得到竖直方向速度与时间的关系为:

𝑣𝑦(𝑡) = (𝑚𝑔/𝑘) * (1 - exp(-𝑘/𝑚 * 𝑡))

其中，𝑚、𝑔和𝑘为已知参数，可以直接赋值，然后利用上面的代码进行计算。

我们可以使用MATLAB的ode45函数求解竖直方向上的运动方程，得到物资在竖直方向上的速度随时间的变化关系 𝑣𝑦(𝑡)。

% 设置参数
m = 50;     % 物资质量，单位：kg
r = 0.2;    % 物资半径，单位：m
A = pi*r^2; % 物资受力面积，单位：m^2
g = 9.8;    % 重力加速度，单位：m/s^2
k = 0.1;    % 空气阻力系数
v_x0 = 100; % 水平速度，单位：m/s
v_y0 = 0;   % 竖直速度，单位：m/s
t_span = [0 10]; % 时间区间

% 竖直方向运动方程
f = @(t, v_y) (-m*g - k*A*v_y*abs(v_y))/m;

% 初始条件
v0 = [v_y0];

% 使用ode45函数求解微分方程
[t, v] = ode45(f, t_span, v0);

% 绘制速度-时间图像
plot(t, v);
xlabel('时间（s）');
ylabel('竖直方向速度（m/s）');
title('竖直方向速度-时间图像');

 水平：

根据上述公式，vx(t)可以表示为：

vx(t) = v0 * exp(-k*t/m) * cos(theta)

以下是用MATLAB代码求解水平方向上的运动方程得到水平方向速度与时间的关系：

% 定义参数
m = 50; % 物资质量，单位kg
r = 0.2; % 物资半径，单位m
rho = 1.2; % 空气密度，单位kg/m^3
A = pi*r^2; % 物资表面积，单位m^2
Cd = 0.5; % 物资阻力系数
v0 = 30; % 无人机飞行速度，单位m/s
h0 = 200; % 无人机飞行高度，单位m
theta = 45; % 投放角度，单位度

g = 9.8; % 重力加速度，单位m/s^2
v0x = v0*cosd(theta); % 投放时刻物资的水平速度，单位m/s
v0y = v0*sind(theta); % 投放时刻物资的竖直速度，单位m/s

% 定义水平方向上的运动方程
f = @(t, v) (-0.5*rho*Cd*A/m)*(v^2);

% 求解水平方向速度与时间的关系
[t, v] = ode45(f, [0, 100], v0x);

% 绘制速度-时间曲线
figure();
plot(t, v);
xlabel('Time (s)');
ylabel('Horizontal Velocity (m/s)');
title('Horizontal Velocity vs Time');

根据投放物资的水平运动方程，有

dx(t) = vx(t) * dt

将vx(t)代入上式，得到

dx(t) = v0 * cos(theta) * exp(-k*t/m) * dt

将上式两端从0到d积分，得到

d = integral(v0 * cos(theta) * exp(-k*t/m), 0, t)

将t代入上式即可得到投放距离d。最终的表达式为：

d = v0 * cos(theta) * m/k * (1 - exp(-k*t/m))

1——（2）

求d：

% 定义常量
m = 50; % 物资质量
r = 0.2; % 物资半径
g = 9.8; % 重力加速度
k = 0.1; % 空气阻力系数

% 定义初始条件
h = 300; % 飞行高度
v0 = 300; % 飞行速度
theta = pi/4; % 投放角度    0   pi/2  pi  3*pi/2

% 计算水平初速度和竖直初速度
vx0 = v0*cos(theta);
vy0 = v0*sin(theta);

% 求解竖直方向运动方程
t = linspace(0, 20, 1000); % 时间变化范围
v_y = -(m*g/k)*ones(1, length(t)) + ((m*g)/k + vy0)*exp(-k*t/m); % 竖直方向速度随时间的变化关系
h = -(m/(2*k))*log(1-(k/m)*v_y); % 高度随时间的变化关系

% 计算落地时间
idx = find(h<=0, 1, 'first'); % 第一个小于等于0的索引
t_ground = t(idx); % 落地时间

% 求解水平方向运动方程
x = zeros(1, length(t)); % 水平方向位移
for i = 2:length(t)
    x(i) = x(i-1) + vx0*(t(i)-t(i-1)) - (k/m)*vx0^2/g*(t(i)-t(i-1)) + ...
           (k/m)*(vx0-g*(t(i)-t(i-1))/2)^2/g*(exp(-2*k*(t(i)-t(i-1))/m)-1);
end

% 计算投放距离
d = x(idx); % 投放距离

% 输出结果
disp(['落地时间：', num2str(t_ground), '秒']);
disp(['投放距离：', num2str(d), '米']);

问题2

2——（1）

建模：

对于无人机发射距离与无人机的飞行高度、飞行速度、俯冲角度及发射速度等因素之间的关系，可以按照以下步骤建立数学模型：

1. 设无人机的初始位置为坐标原点，发射点的坐标为 (d, h)。

2. 根据无人机的俯冲角度 𝜃, 可以得到爆炸物的初始水平速度和竖直速度分量： 𝑣𝑥0 = 𝑣0 × cos(𝜃) 𝑣𝑦0 = 𝑣0 × sin(𝜃) − 𝑣1

3. 对于爆炸物的竖直运动，可以建立如下的运动方程： h + 𝑣𝑦0 × t − 1/2 × g × t^2 = 0 其中 g 为重力加速度，取值为 9.8 m/s^2。

解这个方程可以得到爆炸物落地的时间 t，即： t = (𝑣𝑦0 + √(𝑣𝑦0^2 + 2 × g × h)) / g

1. 对于爆炸物的水平运动，可以建立如下的运动方程： d + 𝑣𝑥0 × t = 0

解这个方程可以得到爆炸物的发射距离 d，即： d = −𝑣𝑥0 × t

1. 将步骤 3 中求得的 t 代入步骤 4 中，可以得到发射距离 d 与无人机的飞行高度、飞行速度、俯冲角度及发射速度等因素之间的关系。

综上所述，无人机发射距离与无人机的飞行高度、飞行速度、俯冲角度及发射速度等因素之间的关系可以用以下的数学模型表示： d = −(𝑣0 × cos(𝜃)) × (𝑣0 × sin(𝜃) − 𝑣1 + √((𝑣0 × sin(𝜃) − 𝑣1)^2 + 2 × g × h)) / g

竖直：

我们可以通过竖直方向上的运动方程求解爆炸物的竖直速度与时间的关系。考虑到在发射时刻，爆炸物在竖直方向上的速度为负的发射速度 𝑣1，因此可以将爆炸物的竖直速度表示为：

𝑣𝑦(𝑡) = −𝑣1 + 𝑔𝑡

其中 g 是重力加速度，t 表示爆炸物的运动时间。

我们可以将上式表示为标准的二次函数形式：

𝑣𝑦(𝑡) = − 0.5 𝑔𝑡^2 − 𝑣1𝑡

现在我们用 MATLAB 求解 𝑣𝑦(𝑡) 与 t 的关系，假设发射速度为 600 km/h，转换为 m/s 为：

𝑣1 = 600 × 1000/3600 = 166.7 m/s

重力加速度为 9.8 m/s^2。

可以用以下 Matlab 代码求解：

g = 9.8; % 重力加速度
v1 = 166.7; % 爆炸物发射速度
t = linspace(0, 50, 1000); % 设定时间范围和时间点数
vy = -0.5 * g * t.^2 - v1 * t; % 爆炸物竖直速度与时间的关系
plot(t, vy); % 绘制竖直速度随时间的变化曲线
xlabel('时间 (s)'); ylabel('竖直速度 (m/s)'); title('竖直速度随时间的变化');

vy​(t)=v0​sinθ−gt

求t

syms t g v0y h0
h = h0 + v0y*t - 0.5*g*t^2;
eqn = h == 0;
t = solve(eqn, t);

 

水平：

 

% 输入参数
v0 = 300/3.6; % 飞行速度，单位：m/s
theta = 45; % 俯冲角度，单位：度
v1 = 600/3.6; % 发射速度，单位：m/s
dmin = 1000; % 最小发射距离，单位：m
dmax = 3000; % 最大发射距离，单位：m
h0 = 800; % 飞行高度，单位：m
g = 9.8; % 重力加速度，单位：m/s^2
r = 0.08; % 爆炸物半径，单位：m
m = 5; % 爆炸物重量，单位：kg

% 计算竖直方向速度分量
vx0 = v0*cosd(theta); % 水平方向初速度，单位：m/s
vy0 = v0*sind(theta) - v1; % 竖直方向初速度，单位：m/s

% 计算落地时间
tf = (-vy0 + sqrt(vy0^2 + 2*g*h0))/g; % 落地时间，单位：s

% 计算水平方向的运动
t = 0:0.01:tf; % 时间范围，步长：0.01s
dx = vx0*t; % 水平方向位移，单位：m

% 判断发射距离是否满足要求
d = dx(end);
if d < dmin || d > dmax
    disp('发射距离不在指定范围内');
else
    disp(['发射距离为 ', num2str(d), ' m']);
end

结果：

发射距离为 1641.077 m

2——（2）