双摆RK4

Question

我正在尝试实施经典的 RK4 算法来求解控制运动方程的微分方程组。但是，我有一段时间遇到了几个问题。运动方程可写为：

由于这是一个一阶微分方程系统，我们准备使用 RK4 来解决它。 但是，对于用 python 代码编写上面的系统，我真的很迷茫。

这里有一个相关的问题：System of second order ODEs Runge Kutta 4th order 将不同物理系统的 ODES 系统写入 Python 代码，并使用 RK4 解决。但是，我无法使用它来使我自己的代码工作。还有另一个问题I want to have the pendulum blob in my double pendulum 更接近我想要的，但它并没有通过使用数组来处理这个过程，而是通过编写几个方程来代替。

import matplotlib.pyplot as plt
import numpy as np


def RK4(t, y, h, f):  
    
    #Runge Kutta standard calculations
    k1 = f(t, y)
    k2 = f(t + h/2, y + h/2 * k1)
    k3 = f(t + h/2, y + h/2 * k2)
    k4 = f(t + h, y + h * k3)    

    return 1/6*(k1 + 2 * k2 + 2 * k3 + k4) 

    
def RHS(t, y):
    
    ################################
    #Critical physical Parameters
    g = 9.8
    l1 = 1
    l2 = 1
    m1 = 1
    m2 = 1
    w1 = w2 = 0
    theta_1 = theta_2 = np.pi/4
    delta_theta = theta_1 - theta_2
    ################################
    
    #Writing the system of ODES.
    f0 = w1
    f1 = w2

    f2 = (m2*l1*w1**2*np.sin(2*delta_theta) + 2*m2*l2*w2**2*np.sin(delta_theta) + 2*g*m2*np.cos(theta_2)*np.sin(delta_theta)
    + 2*g*m1*np.sin(theta_1))/(-2*l1*(m1 + m2*np.sin(delta_theta)**2))
    
    f3 = (m2*l2*w2**2*np.sin(2*delta_theta) + 2*(m1 + m2)*l1*w1**2*np.sin(delta_theta) 
    + 2*g*(m1 + m2)*np.cos(theta_1)*np.sin(delta_theta))/(2*l2*(m1 + m2*np.sin(delta_theta)**2))
    
    return np.array([f0, f1, f2, f3])


def main():
    #Specify time interval and number of domain subintervals
    t0 = 0 
    T = 100
    n = 2048
    # initial conditions
    y0 = [np.pi/4, np.pi/4, 0, 0]                  
    
    #Domain discretization
    t_n = np.linspace(t0, T) 
    y_n = [np.array(y0)]  
    
    #Step size
    h  = (T - t0)/n        
    
    while t_n[-1] < T:
        #Keep going until T is reached.
        #Store numerical solutions into y_n
        y_n.append(y_n[-1] + h * RK4(t_n[-1], y_n[-1], h, RHS))
        t_n.append(t_n[-1] + h)
    
    
    print(y_n)
 

main()

但是，终端给了我以下输出：

[array([0.78539816, 0.78539816, 0.        , 0.        ])]

这可能表明方程组没有被求解。这是为什么？我想我也没有适当地通过初始条件，我只是不知道如何正确地做到这一点。

我真的很难做到这一点。有人可以帮我修复代码以适当地集成 ODES 系统吗？

提前致谢，卢卡斯

Answer 1

您没有在导数函数中使用传递的状态变量，您只是将状态设置为常量。你应该有这样一行

theta1, theta2, w1, w2 = y

改为。