耦合非线性ODE组RK4求解非预期发散 算法修正求助

问题描述 #

现有2个二阶非线性常微分方程，经状态空间法转换得到4个耦合的一阶ODE，使用四阶龙格-库塔（RK4）方法求解时输出曲线出现非预期发散，经排查为RK4实现逻辑错误。方程中L、Fa项包含状态空间变量，不影响本次问题定位。
方程示意图

原始状态空间方程（存在定义疏漏） #

f1 = @(x2) x2; % = x1'

f2 = @(x1, x2, x3) K/m*(l_0/sqrt((X_d(t, t1, t2, a_x0, X_d0, X_d0_tag)-x1).^2+(Z_d(t, t0, a_z0, Z_d0, Z_d0_tag)-x3).^2)-1)*(x1-X_d(t, t1, t2, a_x0, X_d0, X_d0_tag)) ...
    - 0.5*Rho*A*C_d*(x2-Interpolation(Z_d(t, t0, a_z0, Z_d0, Z_d0_tag), data)).^2*sgn(x2, Z_d(t, t0, a_z0, Z_d0, Z_d0_tag), data)/m;
% = x2'

f3 = @(x1) x1; % = x3' （此处存在定义错误）

f4 = @(x1, x3) K/m*(l_0/sqrt((X_d(t, t1, t2, a_x0, X_d0, X_d0_tag)-x1).^2+(Z_d(t, t0, a_z0, Z_d0, Z_d0_tag)-x3).^2)-1)*(x3-Z_d(t, t0, a_z0, Z_d0, Z_d0_tag))-g;
% x4'

原始错误RK4实现 #

h=0.2; % step size
t_array = 0:h:10;

w = zeros(1,length(t_array)); 
x = zeros(1,length(t_array)); 
y = zeros(1,length(t_array)); 
z = zeros(1,length(t_array));

for i=1:(length(t_array)-1) % calculation loop
    t = 0 +h*i; % A parameter needed for the interpolation in f2
    
    k_1 = f1(x(i));
    k_2 = f1(x(i)+0.5*h*k_1);
    k_3 = f1(x(i)+0.5*h*k_2);
    k_4 = f1(x(i)+k_3*h);
    x(i+1) = x(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  
    disp(x(i+1));
    
    m_1 = f3(z(i));
    m_2 = f3(z(i)+0.5*h*k_1);
    m_3 = f3(z(i)+0.5*h*k_2);
    m_4 = f3(z(i)+k_3*h);
    z(i+1) = z(i) + (1/6)*(m_1+2*m_2+2*m_3+m_4)*h;  

    n_1 = f2(x(i), z(i), w(i));
    n_2 = f2(x(i), z(i) ,w(i)+0.5*h*k_1);
    n_3 = f2(x(i), z(i) ,w(i)+0.5*h*k_2);
    n_4 = f2(x(i), z(i) ,w(i)+k_3*h);
    w(i+1) = w(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  

    l_1 = f4(x(i), z(i));
    l_2 = f4(x(i), z(i));
    l_3 = f4(x(i), z(i));
    l_4 = f4(x(i), z(i));
    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  
    
end

核心错误点 #

• 状态导数定义错误：z向位移x3的导数应为z向速度x4，原f3错误写为@(x1) x1，完全不符合物理关系；f2、f4的参数列表未完整包含耦合的状态变量，导数计算时无法获取正确的状态输入。
• RK4迭代逻辑错误：耦合方程组的RK4不能逐变量单独计算斜率，所有状态的斜率必须在同一子步（k1/k2/k3/k4阶段）基于同一组预测状态同步计算，否则会完全丢失变量间的耦合关系。原代码逐变量单独算k值，且速度项更新时错误复用位移项的k值，z向速度的l系列斜率甚至全程使用相同初始状态计算，完全没有实现RK4的子步预测逻辑。
• 时变参数取值错误：计算k2/k3/k4子步斜率时，时间t未对应调整为中点、终点时刻，导致时变的X_d、Z_d、插值项取值错误。

修正后实现代码 #

% 1. 修正统一导数函数，输入为(当前时间t, 状态向量state)
% 状态向量顺序：state(1)=x1(x向位移), state(2)=x2(x向速度), state(3)=x3(z向位移), state(4)=x4(z向速度)
ode_fun = @(t, state) [
    state(2); % x1' = x2 对应原f1
    K/m*(l_0/sqrt((X_d(t, t1, t2, a_x0, X_d0, X_d0_tag)-state(1)).^2 + (Z_d(t, t0, a_z0, Z_d0, Z_d0_tag)-state(3)).^2)-1)*(state(1)-X_d(t, t1, t2, a_x0, X_d0, X_d0_tag)) ...
    - 0.5*Rho*A*C_d*(state(2)-Interpolation(Z_d(t, t0, a_z0, Z_d0, Z_d0_tag), data)).^2*sgn(state(2), Z_d(t, t0, a_z0, Z_d0, Z_d0_tag), data)/m; % x2' 对应原f2
    state(4); % x3' = x4 修正原f3的定义错误
    K/m*(l_0/sqrt((X_d(t, t1, t2, a_x0, X_d0, X_d0_tag)-state(1)).^2 + (Z_d(t, t0, a_z0, Z_d0, Z_d0_tag)-state(3)).^2)-1)*(state(3)-Z_d(t, t0, a_z0, Z_d0, Z_d0_tag)) - g; % x4' 对应原f4
];

% 2. RK4求解配置
h = 0.2; % 步长，若仍存在数值震荡可缩小至0.05~0.1
t_array = 0:h:10;
state_hist = zeros(4, length(t_array));
% 填入自定义初始条件：state_hist(:,1) = [x1初始值; x2初始值; x3初始值; x4初始值];

% 3. 同步迭代RK4（所有状态斜率同子步计算）
for i = 1:(length(t_array)-1)
    t_cur = t_array(i);
    s_cur = state_hist(:,i);
    
    k1 = ode_fun(t_cur, s_cur);
    k2 = ode_fun(t_cur + 0.5*h, s_cur + 0.5*h*k1);
    k3 = ode_fun(t_cur + 0.5*h, s_cur + 0.5*h*k2);
    k4 = ode_fun(t_cur + h, s_cur + h*k3);
    
    state_hist(:,i+1) = s_cur + (h/6)*(k1 + 2*k2 + 2*k3 + k4);
end

% 4. 提取结果
x_disp = state_hist(1,:); % x向位移
x_vel = state_hist(2,:); % x向速度
z_disp = state_hist(3,:); % z向位移
z_vel = state_hist(4,:); % z向速度

排查补充 #

• 修正后若仍有发散，优先检查sgn函数、Interpolation函数的返回值逻辑，尤其是阻力项的符号是否符合物理预期，符号反向会直接导致加速度方向错误引发发散；
• 可使用Matlab自带ode45传入相同的ode_fun求解，将结果与自写RK4结果对比，快速验证实现正确性。

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