将龙格库塔应用于耦合方程

Question

所以我有 2 个二阶非线性 ODE，应用状态空间定理后我有 4 个一阶 ODE。 我正在尝试应用 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

正如我所说，我的图表正在发生变化（它们不应该发生变化），所以我怀疑我的代码是错误的。 请帮我修复算法。 非常感谢！

Answer 1

XD和ZD是什么？为什么他们没有与之相关的微分方程？您能提供更多详细信息吗？

此外，如果您使用向量作为状态变量，并且使用一个函数句柄来生成向量输出，这也会有所帮助。这减少了您编写的代码，并允许您将结果与 ode45( ) 进行比较，因为它可以采用相同的函数句柄作为输入。

代码中的根本缺陷似乎是，尽管您声明方程是耦合的，但您试图将它们零碎地积分。例如，您使用 k_1、k_2、k_3、k_4 对一个变量执行 RK4 方案，将 x(i) 传播到 x(i+1)。但在此过程中，所有其他耦合变量 z(i) 、 w(i) 和 y(i) 在代码中保持静态。这是一个缺陷。所有耦合变量需要同时通过这些中间计算进行传播。即，您需要首先生成所有 k_1、m_1、n_1 和 l_1，然后使用这些结果计算所有 k_2、m_2、n_2 和 l_2。然后使用这些结果计算所有 k_3、m_3、n_3 和 l_3。然后使用所有这些结果计算 k_4、m_4、n_4 和 l_4。最后使用所有这些将所有变量向前传播一步。这就是矢量函数句柄可以为您提供很大帮助的地方。通过创建一个接受向量输入（向量的每个元素代表变量之一）并返回向量输出的函数句柄，您可以将代码简化为仅编写一组 RK4 方程，该方程自动将所有变量向前传播同时，因为它们都是同一向量的一部分。这也将使您的代码更易于调试。

最后，您将混合变量和导数。 x 应该与 k 一起使用，z 应该与 m 一起使用，w 应该与 n 一起使用，y 应该与 l 一起使用。特别是，l 甚至没有实施 RK4 方案。