针对scipy计算的完整的运动方程，选取不同初始条件 ， 研究不同初始条件时的单摆运动图像。用python写出具体代码

为了研究单摆运动的图像，我们可以使用`scipy.integrate.solve_ivp`来解决常微分方程，并结合`matplotlib`来绘制图形。这里我们假设单摆的运动方程是简谐振动方程，即：

\[
\ddot{\theta} + \frac{g}{L}\sin(\theta) = 0
\]

其中，\( \theta \) 是摆角，\( g \) 是重力加速度，\( L \) 是摆长。

以下是使用Python的示例代码：

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
from math import sin, cos, pi

# 定义单摆的运动方程
def single_pendulum(theta, t, g, L):
    thdot = theta[1]
    thetaddot = -g/L * sin(theta[0])
    return [thdot, thetaddot]

# 参数设置
g = 9.81  # 重力加速度
L = 1.0  # 摆长
initial_angles = [(pi/6, 0), (pi/4, 0), (pi/3, 0)]  # 不同的初始摆角和初始速度
timespan = (0, 10)  # 观察时间段（秒）

# 解决每个初始条件下的微分方程
solns = []
for angle in initial_angles:
    init_theta = angle[0]
    init_thetadot = angle[1]
    sol = spi.solve_ivp(single_pendulum, timespan, [init_theta, init_thetadot], args=(g, L))
    solns.append(sol)

# 绘制结果
fig, axs = plt.subplots(len(initial_angles), 2, figsize=(10, len(initial_angles)*5), sharex=True, sharey=True)
for i, (angle, ax_row) in enumerate(zip(initial_angles, axs)):
    ax_angle, ax_speed = ax_row
    ax_angle.set_title(f"Initial Angle: {angle[0]} rad, Initial Speed: {angle[1]} rad/s")
    
    # 绘制角度
    time = solns[i].t
    theta = solns[i].y[:, 0]
    ax_angle.plot(time, theta, label=f'Angle ({angle[0]:.2f}, {angle[1]:.2f})')

    # 绘制速度
    thetadot = solns[i].y[:, 1]
    ax_speed.plot(time, thetadot, label=f'Speed ({angle[1]:.2f})')

    ax_angle.legend()
    ax_angle.grid(True)
    ax_speed.legend()
    ax_speed.grid(True)

plt.xlabel("Time (s)")
plt.tight_layout()
plt.show()