Matlab风选模拟：Euler法正常但RK4法轨迹异常排查求助

我用Matlab开发风选模拟器，Euler数值解法结果符合预期，但RK4法输出异常：谷物轨迹无法显示，谷壳呈现自由下落状态。以下是问题排查和修复方案：

核心问题定位与修复 #

1. RK4中间步状态变量传递错误 #

在rungeKutta4函数中，计算k2、k3、k4时，未正确更新完整状态向量（包含y、v、x、w），仅使用初始的v和w值，导致中间步导数计算完全错误，无法正确模拟粒子运动。

错误代码片段：

k2 = func(t(i) + dt/2, [y(i), v, x(i), w] + dt/2 * k1', yc, denf, u0, h, denp, dp, visf, g, Vp);
k3 = func(t(i) + dt/2, [y(i), v, x(i), w] + dt/2 * k2', yc, denf, u0, h, denp, dp, visf, g, Vp);
k4 = func(t(i) + dt, [y(i), v, x(i), w] + dt * k3', yc, denf, u0, h, denp, dp, visf, g, Vp);

修正方法：
维护完整状态向量，每次计算中间步时，基于当前状态和k值更新得到完整中间状态，再传入函数计算导数；同时k是列向量，无需转置，直接做向量加法。

修正后代码：

% 维护完整状态向量
current_state = [y(i); v; x(i); w];

k1 = func(t(i), current_state, yc, denf, u0, h, denp, dp, visf, g, Vp);
k2 = func(t(i) + dt/2, current_state + dt/2 * k1, yc, denf, u0, h, denp, dp, visf, g, Vp);
k3 = func(t(i) + dt/2, current_state + dt/2 * k2, yc, denf, u0, h, denp, dp, visf, g, Vp);
k4 = func(t(i) + dt, current_state + dt * k3, yc, denf, u0, h, denp, dp, visf, g, Vp);

2. RK4垂直阻力项符号错误 #

在computeTrajectory函数中，垂直方向阻力计算与Euler版本不一致：Euler中使用(velocity_final - v_vel(i))（即0 - 粒子垂直速度），但RK4代码直接使用v_particle，导致阻力方向完全相反，破坏浮力与阻力平衡，谷壳呈现自由下落状态。

错误代码片段：

Rep = denf * sqrt((uf - w_particle)^2 + v_particle^2) * dp / visf;
% ...
Fdrag_y = (0.5 * pi * denf * dp^2 * Cd * sqrt((uf - w_particle)^2 + v_particle^2) * v_particle);

修正方法：
将v_particle替换为(0 - v_particle)，与Euler版本保持一致：

修正后代码：

Rep = denf * sqrt((uf - w_particle)^2 + (0 - v_particle)^2) * dp / visf;
% ...
Fdrag_y = (0.5 * pi * denf * dp^2 * Cd * sqrt((uf - w_particle)^2 + (0 - v_particle)^2) * (0 - v_particle));

3. RK4函数状态变量管理优化 #

原代码中v和w作为单独变量维护易出错，建议统一用状态向量管理，同时确保循环结束后数组截断正确。

修正后的rungeKutta4函数完整代码：

function [y, x] = rungeKutta4(func, t, initial_conditions, yc, denf, u0, h, denp, dp, visf, g, Vp)
    dt = t(2) - t(1);
    n = length(t);
    % 初始化完整状态数组：y, v, x, w
    state = zeros(4, n);
    state(:, 1) = initial_conditions;

    for i = 1:n-1
        if state(1, i) < yc
            break; % Stop loop when y < yc
        end

        current_state = state(:, i);
        k1 = func(t(i), current_state, yc, denf, u0, h, denp, dp, visf, g, Vp);
        k2 = func(t(i) + dt/2, current_state + dt/2 * k1, yc, denf, u0, h, denp, dp, visf, g, Vp);
        k3 = func(t(i) + dt/2, current_state + dt/2 * k2, yc, denf, u0, h, denp, dp, visf, g, Vp);
        k4 = func(t(i) + dt, current_state + dt * k3, yc, denf, u0, h, denp, dp, visf, g, Vp);

        state(:, i+1) = current_state + dt/6 * (k1 + 2*k2 + 2*k3 + k4);
    end

    % 提取结果并截断
    y = state(1, 1:i);
    x = state(3, 1:i);
end

其他调试建议 #

• 移除全局变量：将global u_0改为函数参数传递，减少代码耦合。
• 添加中间值输出：在computeTrajectory中输出uf、Rep、Fdrag_x、Fdrag_y等关键变量，对比Euler和RK4的计算结果，验证一致性。
• 统一变量命名：保持Euler和RK4代码中变量名一致，便于对比排查。

内容的提问来源于stack exchange，提问作者user22975163