如何基于现有Python Runge-Kutta代码实现耦合微分方程组及n维通用求解？

没问题！咱们先从你现有的单方程代码入手，修改适配你给出的耦合方程组，再拓展到适用于任意n维耦合系统的通用实现。

首先得理清变量名：你原来的代码里自变量是x，但现在的耦合方程组自变量是t，因变量是x(t)和y(t)——为了避免混淆，我会在代码里明确区分自变量t和因变量x。

一、针对你的耦合方程组的修改实现 #

你的耦合方程组是：

dy/dt = (x - y²)/2
dx/dt = 3x + 3y

核心修改点是把单个因变量换成状态向量（用numpy数组存储所有耦合变量），让导数函数返回对应长度的导数向量，再调整Runge-Kutta的计算逻辑适配向量运算。修改后的完整代码如下：

import numpy as np
import matplotlib.pyplot as plt
import numba
import time

start_time = time.perf_counter()  # 替代已弃用的time.clock()，保证兼容性

@numba.jit(nopython=True)  # nopython模式能最大化加速效果
def derivatives(t, state):
    # state是状态向量：[x(t), y(t)]
    x, y = state
    dxdt = 3 * x + 3 * y
    dydt = (x - y**2) / 2
    return np.array([dxdt, dydt])

def runge_kutta(t0, state0, t_target, h):
    # state0是初始状态向量：[x0, y0]
    n = int((t_target - t0) / h)
    state = np.copy(state0)
    current_t = t0
    for _ in range(n):
        # 计算Runge-Kutta的四个k值（向量形式）
        k1 = h * derivatives(current_t, state)
        k2 = h * derivatives(current_t + 0.5*h, state + 0.5*k1)
        k3 = h * derivatives(current_t + 0.5*h, state + 0.5*k2)
        k4 = h * derivatives(current_t + h, state + k3)
        # 更新状态向量
        state += (k1 + 2*k2 + 2*k3 + k4) / 6
        current_t += h
    return state

def plot_coupled_system(t_start, t_end, steps, initial_state, h):
    t_values = np.linspace(t_start, t_end, steps)
    x_values = []
    y_values = []
    for t in t_values:
        x, y = runge_kutta(t_start, initial_state, t, h)
        x_values.append(x)
        y_values.append(y)
    
    plt.figure(figsize=(10, 5))
    plt.subplot(121)
    plt.plot(t_values, x_values, label='x(t)')
    plt.plot(t_values, y_values, label='y(t)')
    plt.xlabel('t')
    plt.ylabel('Value')
    plt.legend()
    plt.title('Time Series of x(t) and y(t)')
    
    plt.subplot(122)
    plt.plot(x_values, y_values)
    plt.xlabel('x')
    plt.ylabel('y')
    plt.title('Phase Portrait (x vs y)')
    
    plt.tight_layout()
    print(f"Execution time: {time.perf_counter() - start_time:.4f} seconds")
    plt.show()

# 参数设置
t_start, t_end = 0, 10
steps = t_end * 100
initial_state = np.array([1.0, 1.0])  # 初始x(0)=1, y(0)=1
h = 0.002

plot_coupled_system(t_start, t_end, steps, initial_state, h)

关键修改说明：

• 把原来的dydx改成derivatives，接收自变量t和状态向量state，返回包含所有变量导数的数组
• Runge-Kutta函数现在处理向量而非单个数值，所有k值都是向量，状态更新也是向量运算
• 绘图部分新增了相图（x vs y），更直观展示耦合系统的动态行为
• 替换了已弃用的time.clock()为time.perf_counter()，保证代码在新版本Python中正常运行

二、适用于n维耦合微分方程组的通用实现 #

上面的代码其实已经是通用的了！只要做以下两点，就能适配任意n维耦合系统：

1. 定义导数函数：接收自变量和n维状态向量，返回长度为n的导数数组（每个元素对应一个变量的导数）
2. 设置初始状态：用n维numpy数组传入各变量的初始值

比如，如果要解3维耦合系统：

dx/dt = ax + by
dy/dt = cy + z²
dz/dt = dx - e*z

只需要修改derivatives函数：

@numba.jit(nopython=True)
def derivatives(t, state):
    x, y, z = state
    dxdt = a*x + b*y
    dydt = c*y + z**2
    dzdt = d*x - e*z
    return np.array([dxdt, dydt, dzdt])

然后初始状态改成np.array([x0, y0, z0])，其他代码完全不用改——Runge-Kutta的向量运算会自动处理n维的情况。

额外优化建议：

• 用numba.jit(nopython=True)可以大幅加速迭代计算，尤其是n较大或步长很小的时候
• 如果需要更高的精度，可以考虑用自适应步长的Runge-Kutta方法（比如RK45），但基础的RK4已经能满足大多数场景需求

内容的提问来源于stack exchange，提问作者Amirhossein Rezaei