Python（SciPy）耦合ODE相图与预期结果不匹配的问题排查

嘿，我来帮你一步步揪出相图和预期不符的问题，从代码里的小细节到方程实现，咱们逐一梳理：

一、先搞定几个容易忽略的小错误 #

首先看你的代码，有几个明显的问题，可能直接影响了结果：

1. 忘了导入numpy：代码里用到了np.pi、np.linspace这些，但开头没写import numpy as np——你能运行起来估计是之前在环境里导入过，但最好加上，不然换个环境就报错啦。
2. theta参数写错了！：题目里给定的theta是0.1，但你代码里写成了theta = 0.5，这个参数直接影响drc_dt的计算，绝对是导致收敛位置不对的关键原因之一！

二、检查ODE方程的实现是否准确 #

你的耦合系统方程实现可能存在推导偏差，尤其是那两个公共项：
你写的：

common_term1 = ((w * (w * x - x + 1) ** (N - 1)) - 1) / (N * (w - 1))
common_term2 = (((w * x - x + 1) ** (N - 1)) - 1) / (N * (w - 1))

其实w*x -x +1可以简化成x*(w-1)+1，这部分没问题，但要仔细核对原方程里的这两个项的分子是否和你写的一致——比如common_term1的分子是w*(...)^N-1还是w*(...)^(N-1)-1？你这里用的是N-1，要确保和原方程完全匹配。

再看dx_dt和drc_dt的表达式：

• dx_dt里的(rc * common_term1) - 1 - (rd * common_term2)，是不是应该写成rc*common_term1 - rd*common_term2 -1？虽然数学上等价，但要确认原方程的顺序和符号有没有错。
• drc_dt里的-x * ((rc * common_term1) - 1)其实和x*(1 - rc*common_term1)是一样的，这部分没问题，但因为theta参数错了，结果还是会跑偏。

三、针对你两个具体疑问的解答 #

1. 为啥粉色轨迹收敛到r_c=2.0-2.5而不是1.5？ #

• 首要原因是theta参数错误：你把0.1写成了0.5，这会大幅改变r_c的演化驱动力，直接让平衡点偏移。
• 收敛判断条件太松：你写的abs(final_rc - alpha) < 0.3，alpha是1.5，也就是1.2到1.8之间才算粉色，但因为参数错了，轨迹实际收敛到2.0-2.5，这部分其实不该归为粉色，等修正参数后，再把阈值调小一点（比如0.1），就能准确识别收敛到1.5的轨迹了。
• 时间跨度不够：你用的t_span=(0,50)，对于时变系统来说，可能还没来得及收敛到稳定状态，试试把时间拉长到200甚至300，看看轨迹会不会最终落到1.5附近。

2. 为啥循环区域（橙色）太细长不紧凑？ #

• 初始条件和分类条件的问题：你设置的初始rc范围是1.5到3.5，x是0到1，有些初始条件可能根本不属于循环区域，但被你的分类条件归为橙色了。等修正参数后，循环区域的形状会自然更接近预期。
• 求解精度和平滑度：虽然你设置了rtol=1e-6,atol=1e-9，但可以试试换个更精准的求解器，比如method='DOP853'，同时把t_eval的点数从1000增加到2000，让轨迹更平滑，看起来就不会那么细长了。
• 方程实现错误：前面的参数错误和方程推导偏差，直接导致了循环区域变形，修正这些后，形状会紧凑很多。

四、修正后的代码参考 #

我把几个关键问题改好了，你可以试试运行这个版本：

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp

# 完全匹配题目给定的参数
N = 4
rd = 0.6
epsilon = 6
alpha = 1.5
beta = 3.5
theta = 0.1  # 修正了这里！
a = 0.1
delta = -np.pi / 2

# 时变环境函数
def w2(t):
    return 1 - 0.5 * np.sin(a * t + delta)

# 耦合ODE系统
def coupled_system(t, state):
    x, rc = state
    w = w2(t)
    x_w_term = x*(w-1) + 1  # 简化重复计算的项
    pow_term = x_w_term ** (N-1)
    common_term1 = (w * pow_term - 1) / (N * (w - 1))
    common_term2 = (pow_term - 1) / (N * (w - 1))

    dx_dt = x * (1 - x) * (rc * common_term1 - rd * common_term2 - 1)
    drc_dt = epsilon * (rc - alpha) * (beta - rc) * (x*(1 - rc*common_term1) + rd * theta * (1 - x) * common_term2)

    return [dx_dt, drc_dt]

# 延长时间跨度，让轨迹充分演化
t_span = (0, 200)
t_eval = np.linspace(*t_span, 2000)

# 调整初始条件，避免极端边界，同时适当减少密度提升运行速度
x0_vals = np.linspace(0.05, 0.95, 20)
rc0_vals = np.linspace(1.6, 3.4, 20)

trajectories = []

for x0 in x0_vals:
    for rc0 in rc0_vals:
        # 用更精准的求解器，提高精度
        sol = solve_ivp(coupled_system, t_span, [x0, rc0], t_eval=t_eval, rtol=1e-8, atol=1e-11, method='DOP853')
        trajectories.append(sol)

plt.figure(figsize=(8, 6))

for sol in trajectories:
    x_traj, rc_traj = sol.y
    final_x, final_rc = x_traj[-1], rc_traj[-1]

    # 收紧颜色判断条件，更精准
    if final_x < 0.01:  # 合作者灭绝
        color = '#1034a6'
    elif final_x > 0.9 and abs(final_rc - alpha) < 0.1:  # 收敛到全合作状态
        color = '#FF1493'
    else:  # 循环行为
        color = 'orange'

    plt.plot(x_traj, rc_traj, color=color, linewidth=0.5)

plt.xlabel("Frequency of Cooperators, x")
plt.ylabel("Multiplication of Cooperators, r_c")
plt.title("Phase Portrait of the System")
plt.xlim([0.0, 1.0])
plt.ylim([1.5, 3.5])
plt.show()

五、额外小建议 #

• 先单独跑一两个初始条件（比如x=0.5, rc=1.5），看看轨迹演化是否符合预期，再批量跑所有初始条件，这样更容易排查问题。
• 可以先画一下w(t)随时间的变化，确认时变环境的函数是对的，避免sin函数的参数（比如a*t的系数、delta的符号）出错。
• 一定要仔细核对原方程的每一个项，符号、系数、指数都不能错，这是耦合ODE相图出错的最常见原因！

备注：内容来源于stack exchange，提问作者Emon Hossain