使用 Runge Kutta 的简单摆周期 - python

Question

我正在尝试使用Python中的Runge-Kutta方法对一个简单的摆进行建模，周期应该是2秒。不仅如此，我的角速度和角位移值也不正确。我不确定我的代码中有什么错误。

这是我尝试过的代码 - v1 是速度步长，t1 是角位移步长。

g = 9.81
pi = math.pi
l = g/pi

t0 = 0
t = t0
t_list = [t0]

v0 = 0
v = v0
v_list = [v0]

theta0 = 0.1
theta = theta0
theta_list = [theta0]

dt = 0.02
tfinal = 20

def theta2(t, v, theta):
  return -(g/l)*math.sin(theta)

def theta1(t, v, theta):
  return v

while t <= tfinal:
  if t == 0:
    v1 = ((dt)/2)*theta2(t, v, theta)
  else:
    v1 = theta1(t-(dt)/2, v, theta) + (dt)*theta2(t, v, theta)
  t1 = theta + (dt)*v1
  
  v = v1
  theta = t1
  t = t + dt

  t_list.append(round(t, 2))
  v_list.append(v)
  theta_list.append(theta)

print(t_list)
print(v_list)
print(theta_list)
ax = plt.plot(t_list, theta_list)
ax = plt.xlabel('time')
ax = plt.ylabel('theta')
#ax = plt.xlim(0,3.5)
ax = plt.show()

这是结果图。周期应为 2 秒。

Answer 1

您有一些问题。主要原因之一是您使用的摆锤长度不正确。要获得摆长度，您需要求解：

T = 2 * pi * sqrt(l/g) ==> l=g/pi^2

自从我参加数值动力系统课程以来也有一段时间了，但这看起来不像龙格-库塔，所以我编写了四阶龙格-库塔方法。这似乎提供了正确的结果。

import math
import matplotlib.pyplot as plt

# Constants
g = 9.81
pi = math.pi
l = g / (pi ** 2)  # Corrected length for a period of 2 seconds

# Initial conditions
t0 = 0
v0 = 0
theta0 = 0.1  # radians

dt = 0.02
tfinal = 20

# Initialize variables
t = t0
v = v0
theta = theta0

# Initialize lists to store data
t_list = [t0]
v_list = [v0]
theta_list = [theta0]

# Define derivative functions
def theta_dot(t, theta, v):
    return v

def v_dot(t, theta, v):
    return - (g / l) * math.sin(theta)

# Runge-Kutta 4th order method
while t <= tfinal:
    # Compute k1
    k1_theta = dt * theta_dot(t, theta, v)
    k1_v = dt * v_dot(t, theta, v)
    
    # Compute k2
    k2_theta = dt * theta_dot(t + dt/2, theta + k1_theta/2, v + k1_v/2)
    k2_v = dt * v_dot(t + dt/2, theta + k1_theta/2, v + k1_v/2)
    
    # Compute k3
    k3_theta = dt * theta_dot(t + dt/2, theta + k2_theta/2, v + k2_v/2)
    k3_v = dt * v_dot(t + dt/2, theta + k2_theta/2, v + k2_v/2)
    
    # Compute k4
    k4_theta = dt * theta_dot(t + dt, theta + k3_theta, v + k3_v)
    k4_v = dt * v_dot(t + dt, theta + k3_theta, v + k3_v)
    
    # Update theta and v
    theta += (k1_theta + 2*k2_theta + 2*k3_theta + k4_theta) / 6
    v += (k1_v + 2*k2_v + 2*k3_v + k4_v) / 6
    
    # Update time
    t += dt
    
    # Append data to lists
    t_list.append(t)
    theta_list.append(theta)
    v_list.append(v)

# Plotting after simulation
plt.plot(t_list, theta_list)
plt.xlabel('Time (s)')
plt.ylabel('Theta (rad)')
plt.title('Simple Pendulum Motion')
plt.show()