课后作业2

Q1

问题

假设空中目标距离炮弹发射点 $S=6000m$ ,高度 $H=1450m$ ,做匀速运动,以速度 $v_1=50m/s$ 向发射台运动.炮弹质量为 $1kg$ ,初速度 $v_0=500\sqrt{2}m/s$ ,空气阻尼系数为 $k = 0.1$ ,计算出击中时间以及击中位置。

程序

% --- 主脚本 ---

% 1. 设置求解器的选项
options = optimoptions('fsolve', ...
    'Display', 'iter', ...          % 显示迭代过程
    'TolFun', 1e-8, ...             % 设置容忍度
    'Algorithm', 'levenberg-marquardt'); % 使用L-M算法

% 2. 提供初始猜测值 [t_guess, theta_guess_rad]
initial_guess = [10, 0.44]; % t=10秒, theta=0.44弧度(约25度)

% 3. 调用 fsolve 求解
%    @projectile_system 是我们之前创建的函数的句柄
[solution, fval] = fsolve(@projectile_system, initial_guess, options);

% 4. 显示结果
if norm(fval) < 1e-4 % 检查解是否收敛
    time_to_impact = solution(1);
    angle_rad = solution(2);
    angle_deg = rad2deg(angle_rad);

    fprintf('\n--- 求解结果 ---\n');
    fprintf('炮弹击中目标所需时间 (t): %.4f 秒\n', time_to_impact);
    fprintf('炮弹发射角度 (θ): %.4f 度\n', angle_deg);

    % 5. 计算击中位置
    S = 6000;
    v1 = 50;
    H = 1450;
    x_impact = S - v1 * time_to_impact;
    y_impact = H; % 击中高度即为目标高度

    fprintf('击中位置坐标 (x, y): (%.2f m, %.2f m)\n', x_impact, y_impact);
    
    % 6. 绘制轨迹图
    plot_trajectories(time_to_impact, angle_rad, x_impact, y_impact);
else
    fprintf('\n求解失败，求解器未能收敛到一个解。\n');
    disp('最后的函数值 (误差):');
    disp(fval);
end

function F = projectile_system(vars)
    % vars(1) 是时间 t
    % vars(2) 是发射角 theta (单位：弧度)

    % 已知参数
    m = 1;          % kg
    k = 0.1;        % 空气阻力系数
    v0 = 500*sqrt(2); % m/s
    g = 9.8;        % m/s^2
    S = 6000;       % m
    H = 1450;       % m
    v1 = 50;        % m/s

    % 从输入向量中获取 t 和 theta
    t = vars(1);
    theta = vars(2);

    % 方程1: 水平方向
    x_projectile = (m/k) * v0 * cos(theta) * (1 - exp(-k*t/m));
    x_target = S - v1*t;
    F(1) = x_projectile - x_target;

    % 方程2: 竖直方向 (使用整理后的形式)
    % 原始形式: (m/k)*v0*sin(theta)*(1-exp(-k*t/m)) - (m*g/k)*t - (m^2*g/k^2)*(exp(-k*t/m)-1) - H
    y_projectile_term1 = (m/k) * v0 * sin(theta) * (1 - exp(-k*t/m));
    y_projectile_term2 = (m*g/k) * t;
    y_projectile_term3 = (m^2*g/k^2) * (exp(-k*t/m) - 1);
    
    F(2) = y_projectile_term1 - y_projectile_term2 - y_projectile_term3 - H;
end

function plot_trajectories(time_to_impact, angle_rad, x_impact, y_impact)
    % 参数定义
    m = 1;          % kg
    k = 0.1;        % 空气阻力系数
    v0 = 500*sqrt(2); % m/s
    g = 9.8;        % m/s^2
    S = 6000;       % m
    H = 1450;       % m
    v1 = 50;        % m/s
    
    % 创建时间向量
    t_vector = linspace(0, time_to_impact, 100);
    
    % 计算炮弹轨迹
    x_projectile = (m/k) * v0 * cos(angle_rad) * (1 - exp(-k*t_vector/m));
    y_projectile = (m/k) * v0 * sin(angle_rad) * (1 - exp(-k*t_vector/m)) ...
                   - (m*g/k) * t_vector ...
                   - (m^2*g/k^2) * (exp(-k*t_vector/m) - 1);
    
    % 计算目标物轨迹
    x_target = S - v1 * t_vector;
    y_target = H * ones(size(t_vector));
    
    % 创建图形
    figure('Position', [100, 100, 700, 500]);
    
    % 子图1: 整体轨迹
    plot(x_projectile, y_projectile, 'b-', 'LineWidth', 2, 'DisplayName', '炮弹轨迹');
    hold on;
    plot(x_target, y_target, 'r--', 'LineWidth', 2, 'DisplayName', '目标轨迹');
    plot(x_impact, y_impact, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'red', ...
         'DisplayName', sprintf('击中点 (%.1f, %.1f)', x_impact, y_impact));
    
    % 标记起点
    plot(0, 0, 'go', 'MarkerSize', 8, 'MarkerFaceColor', 'green', 'DisplayName', '发射点');
    plot(S, H, 'mo', 'MarkerSize', 8, 'MarkerFaceColor', 'magenta', 'DisplayName', '目标起点');
    
    xlabel('水平距离 (m)');
    ylabel('高度 (m)');
    title('炮弹与目标物轨迹');
    legend('Location', 'best');
    grid on;
    axis equal;
    
    % 添加信息文本框
    info_str = sprintf(['求解结果:\n' ...
                       '发射角度: %.2f°\n' ...
                       '飞行时间: %.2f s\n' ...
                       '击中位置: (%.1f, %.1f)\n' ...
                       '目标速度: %.1f m/s\n' ...
                       '炮弹初速: %.1f m/s'], ...
                       rad2deg(angle_rad), time_to_impact, ...
                       x_impact, y_impact, v1, v0);
    
    annotation('textbox', [0.02, 0.02, 0.3, 0.2], 'String', info_str, ...
               'BackgroundColor', 'white', 'EdgeColor', 'black', ...
               'FontSize', 10, 'VerticalAlignment', 'bottom');
end

输出

>> q1

                                        First-order                     Norm of 
 Iteration  Func-count   ||f(x)||^2      optimality       Lambda           step
     0           3      2.12856e+06        3.15e+06         0.01
     1           6            81596        6.76e+05        0.001        4.57997
     2           9          534.633        4.88e+04       0.0001        1.39244
     3          12         0.019458            21.4        1e-05       0.135578
     4          15      4.49331e-11         0.00676        1e-06    0.000893721
     5          18      8.78879e-25        3.22e-09        1e-07    4.22753e-08

方程已解。

fsolve 已完成，因为按照函数容差的值衡量，
函数值向量接近于零，并且按照梯度的值衡量，
问题似乎为正则问题。

<停止条件详细信息>

--- 求解结果 ---
炮弹击中目标所需时间 (t): 16.1075 秒
炮弹发射角度 (θ): 23.3662 度
击中位置坐标 (x, y): (5194.63 m, 1450.00 m)

分析

理论推导

题目中给出了炮弹击中空中目标的方程组，即：
$\begin{cases} -\frac{m}{k}v_0\cos \theta e^{-\frac{k}{m}t}+\frac{m}{k}v_0\cos \theta =S-\left( v_1t+\frac{1}{2}\alpha t^2 \right)\\ -\frac{m}{k}v_0\sin \theta e^{-\frac{k}{m}t}-\frac{mg}{k}\left( t+\frac{m}{k}e^{-\frac{k}{m}t} \right) +\frac{m}{k}v_0\sin \theta +\frac{m^2g}{k^2}=H\\ \end{cases}$
炮弹击中空中目标需要满足的条件为：
$\begin{cases} x_t\left( t,\theta \right) =S-\left( v_1t+\frac{1}{2}\alpha t^2 \right)\\ y_t\left( t,\theta \right) =H\\ \end{cases}$
题干中提到目标做匀速运动，则 $\alpha=0$ ，则水平方向只需满足
$x_t\left( t,\theta \right) =S-v_1t$
对方程组进行整理可得：
$\begin{cases} x_p=\frac{m}{k}v_0\cos \theta (1-e^{-\frac{k}{m}t})\\ y_p=\frac{m}{k}v_0\sin \theta \left( 1-e^{-\frac{k}{m}t} \right) -\frac{mg}{k}t-\frac{m^2g}{k^2}\left( e^{-\frac{k}{m}t}-1 \right)\\ \end{cases}$
则需要求解的方程组为：
$\begin{cases} x_p\left( t,\theta \right) =x_t\left( t \right)\\ y_p\left( t,\theta \right) =y_t\\ \end{cases}$
即
$\begin{cases} F\left( 1 \right) =x_p\left( t,\theta \right) -x_t\left( t \right) =0\\ F\left( 2 \right) =y_p\left( t,\theta \right) -y_t=0\\ \end{cases}$

将上述理论推导在Matlab中实现，获得 $F(1)$ 和 $F(2)$ ；

再根据物理直觉取一个合理的初始值，这里取initial_guess = [10, 0.44]即预测在 $t=10s,\theta =0.44rad\left( \approx 25\degree \right)$ 时能够击中目标；

最后使用Matlab提供的fsolve函数对方程组进行求解。

‍

数据分析课后作业2

课后作业2

Q1

问题

程序

输出

分析