用一级，二级欧拉法求解双摆摆角和时间关系，并绘制曲线，用python求解，写出相关代码

在Python中，我们可以使用一阶和二阶欧拉方法来近似解决双摆系统的微分方程，双摆系统通常描述两个质量通过轻质杆相连的物理模型。这里，我们将假设使用theta1和theta2分别表示两个摆角，而dt代表时间步长。

首先，我们需要一些基本的数学库，如NumPy用于数值计算和matplotlib用于绘图。以下是一级和二级欧拉法的代码示例：

**一阶欧拉法（Euler's Method）：**

import numpy as np
import matplotlib.pyplot as plt

def double_pendulum(theta1, theta2, L1, L2, m1, m2, g, dt):
    dtheta1_dt = (m2 * g * L2 * np.sin(theta2) - m1 * g * L1 * np.sin(theta1)) / (L1 ** 2 + L2 ** 2 * np.cos(theta1 - theta2) ** 2)
    dtheta2_dt = (m1 * g * L1 * np.sin(theta1) * np.cos(theta1 - theta2) - m2 * g * L2 * np.sin(theta2)) / L2 / np.cos(theta2)

    theta1_new = theta1 + dtheta1_dt * dt
    theta2_new = theta2 + dtheta2_dt * dt

    return theta1_new, theta2_new

# 初始化参数
L1, L2 = 1, 1
m1, m2 = 1, 1
g = 9.81
theta1_0, theta2_0 = np.pi/4, np.pi/2
t = np.linspace(0, 10, num=1000)
dt = t[1] - t[0]

theta1_euler = [theta1_0]
theta2_euler = [theta2_0]

for t_i in t[:-1]:
    theta1_euler.append(double_pendulum(*theta1_euler[-1], L1, L2, m1, m2, g, dt)[0])
    theta2_euler.append(double_pendulum(*theta2_euler[-1], L1, L2, m1, m2, g, dt)[1])

# 绘制曲线
plt.figure()
plt.plot(t, theta1_euler, label='First Pendulum')
plt.plot(t, theta2_euler, label='Second Pendulum')
plt.xlabel('Time (s)')
plt.ylabel('Angle (rad)')
plt.legend()
plt.show()

**二级欧拉法（Midpoint Method）精度更高：**

def double_pendulum_midpoint(theta1, theta2, L1, L2, m1, m2, g, dt):
    # 使用midpoint rule for derivatives
    dtheta1_dt = (m2 * g * L2 * np.sin(theta2) - m1 * g * L1 * np.sin(theta1)) / (L1 ** 2 + L2 ** 2 * np.cos(theta1 - theta2) ** 2)
    dtheta2_dt = (m1 * g * L1 * np.sin(theta1) * np.cos(theta1 - theta2) - m2 * g * L2 * np.sin(theta2)) / L2 / np.cos(theta2)

    midpoint_theta1 = theta1 + dtheta1_dt * dt / 2
    midpoint_theta2 = theta2 + dtheta2_dt * dt / 2

    theta1_new = double_pendulum(midpoint_theta1, midpoint_theta2, L1, L2, m1, m2, g, dt)[0]
    theta2_new = double_pendulum(midpoint_theta1, midpoint_theta2, L1, L2, m1, m2, g, dt)[1]

    return theta1_new, theta2_new

theta1_second_order = [theta1_0]
theta2_second_order = [theta2_0]

for t_i in t[:-1]:
    theta1_second_order.append(double_pendulum_midpoint(*theta1_second_order[-1], L1, L2, m1, m2, g, dt)[0])
    theta2_second_order.append(double_pendulum_midpoint(*theta2_second_order[-1], L1, L2, m1, m2, g, dt)[1])

# 继续绘制曲线
plt.figure()
plt.plot(t, theta1_second_order, label='First Pendulum (Midpoint)')
plt.plot(t, theta2_second_order, label='Second Pendulum (Midpoint)')
plt.xlabel('Time (s)')
plt.ylabel('Angle (rad)')
plt.legend()
plt.show()

这两个版本的代码都只包含了一个简单的一维摆动，实际上双摆系统需要考虑四个状态变量。在实际应用中，你需要对每个角度分别应用上述方法。