二阶系统RK4仿真代码

下面是一个使用RK4方法对二阶系统进行仿真的Python代码示例：

import numpy as np
import matplotlib.pyplot as plt

# 定义二阶系统的微分方程
def system(t, x):
    m = 1.0  # 质量
    k = 2.0  # 弹性系数
    c = 0.5  # 阻尼系数
    x1_dot = x[1]
    x2_dot = (-k/m)*x[0] - (c/m)*x[1]
    return np.array([x1_dot, x2_dot])

# 定义RK4方法
def rk4(t, x, f, h):
    k1 = h * f(t, x)
    k2 = h * f(t + h/2, x + k1/2)
    k3 = h * f(t + h/2, x + k2/2)
    k4 = h * f(t + h, x + k3)
    return x + (1/6)*(k1 + 2*k2 + 2*k3 + k4)

# 初始化参数
t0 = 0.0
tf = 10.0
h = 0.01
t = np.arange(t0, tf + h, h)
x = np.zeros((len(t), 2))
x[0] = np.array([1.0, 0.0])  # 初始状态

# RK4仿真
for i in range(len(t) - 1):
    x[i+1] = rk4(t[i], x[i], system, h)

# 绘制结果
plt.figure()
plt.plot(t, x[:, 0], 'b-', label='x')
plt.plot(t, x[:, 1], 'r-', label='v')
plt.xlabel('Time (s)')
plt.ylabel('Position (m) / Velocity (m/s)')
plt.legend(loc='best')
plt.show()

在这个例子中，我们定义了一个二阶系统的微分方程，并使用RK4方法对系统进行仿真。我们使用numpy和matplotlib库进行计算和绘图。定义RK4方法时，我们使用步长h，将当前状态x和微分方程f作为参数输入，并返回下一个状态。在主循环中，我们使用RK4方法进行仿真，并将结果保存在数组x中。最后，我们绘制x和v（速度）关于时间的曲线。