龙格库塔方法求解单摆微分方程python

以下是使用龙格库塔方法求解单摆微分方程的Python代码：

import numpy as np
import matplotlib.pyplot as plt

# 定义微分方程
def pend(y, t, b, c):
    theta, omega = y
    dydt = [omega, -b*omega - c*np.sin(theta)]
    return dydt

# 定义龙格库塔方法
def rk4(y, t, dt, derivs, *args):
    k1 = dt*derivs(y, t, *args)
    k2 = dt*derivs(y + 0.5*k1, t + 0.5*dt, *args)
    k3 = dt*derivs(y + 0.5*k2, t + 0.5*dt, *args)
    k4 = dt*derivs(y + k3, t + dt, *args)
    y_next = y + (k1 + 2*k2 + 2*k3 + k4)/6
    return y_next

# 定义初始条件
theta0 = np.pi/4
omega0 = 0
y0 = [theta0, omega0]

# 定义常数
b = 0.25
c = 5

# 定义时间步长和时间点
dt = 0.01
t = np.arange(0, 20, dt)

# 使用龙格库塔方法求解微分方程
y = np.zeros((len(t), 2))
y[0, :] = y0
for i in range(len(t)-1):
    y[i+1, :] = rk4(y[i, :], t[i], dt, pend, b, c)

# 绘制图像
plt.plot(t, y[:, 0], 'b', label='theta(t)')
plt.plot(t, y[:, 1], 'g', label='omega(t)')
plt.legend(loc='best')
plt.xlabel('t')
plt.grid()
plt.show()