使用scipy中的函数求解完整的单摆运动方程(不用小角度近似)。画出 摆角与时间关系的运动图像。用python写出具体代码

在Python中，可以利用`scipy.integrate.solve_ivp`函数来求解非线性的单摆运动方程，这是一个全阶的动力学问题。这里我们将考虑哈密顿函数描述下的完整单摆方程：

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

其中，$\theta$是摆角，$g$是重力加速度，$l$是摆长。

首先，我们需要导入所需的库并定义初始条件和函数：

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

# 参数值
g = 9.81  # 重力加速度 (m/s^2)
l = 1     # 摆长 (m)

def single_pendulum(t, y):
    theta, omega = y  # y是一个包含当前摆角和角速度的向量
    dtheta_dt = omega  # 角速度对时间的导数
    domega_dt = -g/l * np.sin(theta)  # 角加速度的计算
    return [dtheta_dt, domega_dt]

# 初始条件：假设初始角度很小，但不是小角度近似
initial_angle = np.pi / 4  # 摆角 (rad)
initial_velocity = 0        # 初始角速度 (rad/s)
y0 = [initial_angle, initial_velocity]

# 时间范围
t_span = (0, 10)  # 计算时间范围，单位为秒
t_eval = np.linspace(*t_span, num=1000)  # 等间距的时间点用于绘制图形

# 求解动力学方程
solution = solve_ivp(single_pendulum, t_span, y0, method='LSODA')

# 绘制摆角随时间的变化图
plt.plot(solution.t, solution.y[0], label="摆角")
plt.xlabel('Time (s)')
plt.ylabel('Angle (rad)')
plt.title('完整单摆摆角与时间的关系')
plt.legend()
plt.show()