双摆RK4

Question

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

由于这是一个一阶微分方程组，我们准备用RK4求解。然而，我真的迷失了用 python 代码编写上面的系统。

这里有一个相关的问题：System of second order ODEs Runge Kutta 4th order which writes the system of ODES for a different physical system into Python code and also uses RK4 to solve it. However, I have not been able to use that to make my own code work. There is another question 这更接近我想要的，但它并没有通过使用数组来处理这个过程，而是通过写几个方程来代替。

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

相反。

双摆RK4

Double Pendulum RK4

python

simulation

ode

differential-equations

runge-kutta