用四个人微分方程，画出双摆中第二个摆球的运动定律，并写出代码

在物理中，双摆系统通常涉及两个摆球通过一根轻质杆相连，它们各自受到重力影响并相互作用。要描述这种系统的运动，我们可以使用简化的模型，即每个摆球都由四个自由度（两个角度θ1和θ2代表两根摆线相对于竖直方向的偏转，以及各自的摆角速度ω1和ω2）来表示。这两个摆球之间的相互作用可以简化为拉伸力，假设杆是刚性的。

对于第二个摆球（设为球2），其运动方程可以用类似于这样的微分方程组来表示：

1. 对于角度θ2的运动：
\( \ddot{\theta}_2 = -\frac{g}{l_2}\sin(\theta_2) - \frac{k}{m_2 l_2} (\theta_2 - \theta_1) \)

2. 对角速度ω2的变化：
\( \ddot{\omega}_2 = \frac{g}{l_2}\cos(\theta_2)\dot{\theta}_2 - \frac{k}{I_2 + m_2 l_2^2} (\omega_2 - \omega_1) \)

这里，\( g \)是重力加速度，\( l_2 \)是第二个摆长，\( k \)是连接杆对摆动的影响系数，\( m_2 \)是第二球的质量，\( I_2 \)是第二球的转动惯量，\( \theta_1 \)是第一个摆球的角度。

为了模拟这个系统，你可以使用数值积分方法（如欧拉法、龙格-库塔法等）编写程序，例如Python中的`scipy.integrate.solve_ivp`函数。以下是使用Python的简单示例代码框架：

import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt

# 定义常数和初始条件
g = 9.81
l1 = l2 = 1  # 摆长假设相等
k = ...  # 杆的作用系数
m2 = 1
I2 = ...  # 第二球的转动惯量
theta1_0, theta2_0, omega1_0, omega2_0 = ...

def pendulum_equations(t, state):
    theta1, theta2, omega1, omega2 = state
    dtheta1dt = omega1
    dtheta2dt = omega2
    domega1dt = -g / l1 * np.sin(theta1)
    domega2dt = -g / l2 * np.sin(theta2) - k / (m2 * l2**2) * (theta2 - theta1)
    
    return [dtheta1dt, dtheta2dt, domega1dt, domega2dt]

tspan = (0, T)  # 时间范围
initial_state = [theta1_0, theta2_0, omega1_0, omega2_0]
sol = solve_ivp(pendulum_equations, tspan, initial_state, method='RK45')

# 绘制结果
plt.plot(sol.t, sol.y[1], label="θ2")
plt.plot(sol.t, sol.y[3], label="ω2")
plt.xlabel("时间(s)")
plt.ylabel("角度/角速度")
plt.legend()

注意：实际应用中，你需要提供所有未知常数的具体值，并可能需要调整数值积分方法的精度。此外，由于这是一个复杂的问题，可能还需要考虑非线性和离散化的时间步长。