解常微分方程组

Euler法

function euler_method

    % 参数设置

    v_missile = 450; % 导弹速度 km/h

    v_enemy = 90; % 敌艇速度 km/h

    % 初始条件

    x0 = 0; % 导弹初始位置 x

    y0 = 0; % 导弹初始位置 y

    xe0 = 120; % 敌艇初始位置 y

    t0 = 0; % 初始时间

    % 时间步长和总时间

    dt = 0.01; % 时间步长

    t_final = 0.28; % 终止时间

    % 初始化变量

    t = t0:dt:t_final;

    n = length(t);

    x = zeros(1, n);

    y = zeros(1, n);

    xe = 90 * t;

    ye = xe0 * ones(1, n);

    x(1) = x0;

    y(1) = y0;

    % 欧拉法

    for i = 1:n-1

        dx = xe(i) - x(i);

        dy = ye(i) - y(i);

        distance = sqrt(dx^2 + dy^2);

        x(i+1) = x(i) + 450 * dx / distance * dt;

        y(i+1) = y(i) + 450 * dy / distance * dt;

        if y(i+1)>=xe0% 判定击中条件

            fprintf('欧拉法: 击中时间: %.2f 小时\n', t(i));

            fprintf('欧拉法: 击中位置: (%.2f, %.2f)\n', x(i), y(i));

            break;

        end

    end

[t;x;y]'

    % 绘图

    figure;

    plot(x, y, 'r', xe, ye, 'b');

    legend('导弹轨迹', '敌艇轨迹');

    xlabel('x (km)');

    ylabel('y (km)');

    title('导弹追击敌艇轨迹 - 欧拉法');

    grid on;

end

Matlab运行结果界面：

改进欧拉法：

function improved_euler_method

    % 参数设置

    v_missile = 450; % 导弹速度 km/h

    v_enemy = 90; % 敌艇速度 km/h

    % 初始条件

    x0 = 0; % 导弹初始位置 x

    y0 = 0; % 导弹初始位置 y

    xe0 = 120; % 敌艇初始位置 y

    t0 = 0; % 初始时间

    % 时间步长和总时间

    dt = 0.01; % 时间步长

t_final = 0.30; % 终止时间

    % 初始化变量

    t = t0:dt:t_final;

    n = length(t);

    x = zeros(1, n);

    y = zeros(1, n);

    xe = 90 * t;

    ye = xe0 * ones(1, n);

    x(1) = x0;

y(1) = y0;

    % 改进欧拉法

    for i = 1:n-1

        dx = xe(i) - x(i);

        dy = ye(i) - y(i);

        distance = sqrt(dx^2 + dy^2);

        x_star = x(i) + 450 * dx / distance * dt;

        y_star = y(i) + 450 * dy / distance * dt;

        dx_star = xe(i+1) - x_star;

        dy_star = ye(i+1) - y_star;

        distance_star = sqrt(dx_star^2 + dy_star^2);

        x(i+1) = x(i) + 0.5 * 450 * (dx / distance + dx_star / distance_star) * dt;

        y(i+1) = y(i) + 0.5 * 450 * (dy / distance + dy_star / distance_star) * dt;

        if y(i)>=xe0 % 判定击中条件

            fprintf('改进欧拉法: 击中时间: %.2f 小时\n', t(i));

            fprintf('改进欧拉法: 击中位置: (%.2f, %.2f)\n', x(i), y(i));

            break;

        end

    end

[t;x;y]'

    % 绘图

    figure;

    plot(x, y, 'r', xe, ye, 'b');

    legend('导弹轨迹', '敌艇轨迹');

    xlabel('x (km)');

    ylabel('y (km)');

    title('导弹追击敌艇轨迹 - 改进欧拉法');

    grid on;

end

龙格库塔法：

function runge_kutta_method

    % 参数设置

    v_missile = 450; % 导弹速度 km/h

    v_enemy = 90; % 敌艇速度 km/h

    % 初始条件

    x0 = 0; % 导弹初始位置 x

    y0 = 0; % 导弹初始位置 y

    xe0 = 120; % 敌艇初始位置 y

    t0 = 0; % 初始时间

    % 时间步长和总时间

    dt = 0.01; % 时间步长

    t_final = 0.3; % 终止时间

    % 初始化变量

    t = t0:dt:t_final;

    n = length(t);

    x = zeros(1, n);

    y = zeros(1, n);

    xe = 90 * t;

    ye = xe0 * ones(1, n);

    x(1) = x0;

    y(1) = y0;

    % 龙格库塔法

    for i = 1:n-1

        k1x = 450 * (xe(i) - x(i)) / sqrt((xe(i) - x(i))^2 + (ye(i) - y(i))^2);

        k1y = 450 * (ye(i) - y(i)) / sqrt((xe(i) - x(i))^2 + (ye(i) - y(i))^2);

        k2x = 450 * (xe(i) - (x(i) + 0.5 * dt * k1x)) / sqrt((xe(i) - (x(i) + 0.5 * dt * k1x))^2 + (ye(i) - (y(i) + 0.5 * dt * k1y))^2);

        k2y = 450 * (ye(i) - (y(i) + 0.5 * dt * k1y)) / sqrt((xe(i) - (x(i) + 0.5 * dt * k1x))^2 + (ye(i) - (y(i) + 0.5 * dt * k1y))^2);

        k3x = 450 * (xe(i) - (x(i) + 0.5 * dt * k2x)) / sqrt((xe(i) - (x(i) + 0.5 * dt * k2x))^2 + (ye(i) - (y(i) + 0.5 * dt * k2y))^2);

        k3y = 450 * (ye(i) - (y(i) + 0.5 * dt * k2y)) / sqrt((xe(i) - (x(i) + 0.5 * dt * k2x))^2 + (ye(i) - (y(i) + 0.5 * dt * k2y))^2);

        k4x = 450 * (xe(i) - (x(i) + dt * k3x)) / sqrt((xe(i) - (x(i) + dt * k3x))^2 + (ye(i) - (y(i) + dt * k3y))^2);

        k4y = 450 * (ye(i) - (y(i) + dt * k3y)) / sqrt((xe(i) - (x(i) + dt * k3x))^2 + (ye(i) - (y(i) + dt * k3y))^2);

        x(i+1) = x(i) + (1/6) * dt * (k1x + 2*k2x + 2*k3x + k4x);

        y(i+1) = y(i) + (1/6) * dt * (k1y + 2*k2y + 2*k3y + k4y);

        distance = sqrt((xe(i+1) - x(i+1))^2 + (ye(i+1) - y(i+1))^2);

        if  y(i+1)-y(i)<0.001% 判定击中条件

            fprintf('龙格库塔法: 击中时间: %.2f 小时\n', t(i));

            fprintf('龙格库塔法: 击中位置: (%.2f, %.2f)\n', x(i), y(i));

            break;

        end

    end

[t;x;y]'

    % 绘图

    figure;

    plot(x, y, 'r', xe, ye, 'b');

    legend('导弹轨迹', '敌艇轨迹');

    xlabel('x (km)');

    ylabel('y (km)');

    title('导弹追击敌艇轨迹 - 龙格库塔法');

    grid on;

end

仿真算法：

function fangzhen(h)

L = 120;Vs = 90;Vm = 450;

x(1) = 0;y(1) = 0;

for i = 1:1e6

    M = sqrt(((i - 1)*Vs*h - x(i)).^2 + (L - y(i)).^2);

    ctheta = ((i - 1)*Vs*h - x(i))./M;

    stheta = (L - y(i))./M;

    x(i+1) = x(i) + Vm*h*ctheta;

    y(i+1) = y(i) + Vm*h*stheta;

    if y(i+1) >= L

        break;

    end

end

plot(x,y,x(1):0.01:x(end),L)

x = x(end)

y = y(end)

t = x / Vs

通过本次实验我学会了对实际问题建立模型，然后借助数值计算的方法：EULER法、预报校正法—改进欧拉法、龙格库塔法来求解模型进而解决实际问题。

掌握了欧拉方法、改进欧拉法、龙格库塔法的算法迭代过程，并编写matlab程序并运行程序，这期间我也手工迭代了三次，将两者结果进行比对后发现结果一致，证实迭代过程没有问题；训练了我对实际问题的建模能力和使用数值计算方法解决实际问题的能力。

发现了步长对于结果的影响，此外还发现了自己在建立模型过程中和计算求解时的一些不足。matlab程序编写时的能力需要提高，尤其是提高代码的简洁性且输出想要的结果的能力。今后我的思考问题的能力和给实际问题建立数学模型的能力也有待提高。