用python中的四阶Runge Kutta求解三个耦合的非线性微分方程

import numpy as np

# 定义三个耦合的非线性微分方程
def coupled_ode_system(t, y):
    y1, y2, y3 = y
    dy1dt = -y1 + y2 * y3
    dy2dt = -y2 + y1 * y3
    dy3dt = -y3 + y1 * y2
    return np.array([dy1dt, dy2dt, dy3dt])

# 四阶Runge-Kutta方法
def rk4(func, t_span, y0, N):
    t = np.linspace(t_span[0], t_span[1], N)
    y = np.zeros((N, len(y0)))
    y[0] = y0
    
    for i in range(N - 1):
        h = t[i+1] - t[i]
        k1 = func(t[i], y[i])
        k2 = func(t[i] + h/2, y[i] + h*k1/2)
        k3 = func(t[i] + h/2, y[i] + h*k2/2)
        k4 = func(t[i] + h, y[i] + h*k3)
        y[i+1] = y[i] + h * (k1 + 2*k2 + 2*k3 + k4) / 6
    
    return t, y

# 初始条件和时间范围
y0 = np.array([1.0, 1.0, 1.0])
t_span = [0, 10]
N = 1000

# 求解微分方程
t, y = rk4(coupled_ode_system, t_span, y0, N)

# 打印结果
print("Time:", t)
print("Solution:", y)