用一级，二级欧拉法求解单摆摆角和时间关系曲线，已知摆长1.016米，重力加速度9.81米每秒的平方，摆球质量0.05千克，阻力与运动速度成正比，比例系数为0.02，时间在25秒内。用python求解，写出相关代码

一级和二级欧拉法通常用于数值解算物理方程，特别是在处理连续系统的微分方程时。对于单摆的问题，我们可以使用简谐振动的模型，即θ'' = -g/l * sin(θ) + d(v)/dt，其中θ是摆角，v是摆动速率，l是摆长，g是重力加速度。由于有阻力的影响，我们需要考虑其对速度的影响。

首先，我们忽略空气阻力，这是一个理想化的假设。然后，引入阻力修正的简化模型，将阻力简化为一个线性项，并将其纳入动力学方程。在实际应用中，一级和二级欧拉法分别对应于简单更新（Forward Euler）和改进更新（RK2，如Midpoint Method），这里我们先展示一级欧拉法的代码：

import numpy as np

# 定义常量
l = 1.016  # 摆长，单位米
m = 0.05   # 摆球质量，单位千克
g = 9.81   # 重力加速度，单位 m/s^2
k = 0.02   # 阻力比例系数

# 初始条件（例如摆角0度，初始速度0）
theta_0 = 0.0
omega_0 = 0.0

# 时间步长和总时间
delta_t = 0.01  # 时间步长，单位秒
t_max = 25       # 总时间，单位秒

# 使用一级欧拉法求解
def euler_step(theta, omega, t):
    alpha = -(g / l) * np.sin(theta)
    dv_dt = -k * omega
    omega_new = omega + delta_t * alpha
    theta_new = theta + delta_t * omega
    return theta_new, omega_new

# 创建数组存储结果
times = np.arange(0, t_max, delta_t)
theta_euler = [theta_0]
omega_euler = [omega_0]

for t in times[:-1]:
    theta_euler.append(euler_step(theta_euler[-1], omega_euler[-1], t)[0])
    omega_euler.append(euler_step(theta_euler[-1], omega_euler[-1], t)[1])

# 结果是列表形式，可以画图显示摆角随时间变化
theta_euler_array = np.array(theta_euler)