I’m trying to generate a phase portrait for a system of coupled ordinary differential equations (ODEs) using Python’s scipy.integrate.solve_ivp. The system models the frequency of cooperators (x) and the multiplication of cooperators (r_c) in an evolutionary game theory context with a time-dependent environment w(t) .However, my output graph does not match the expected phase portrait, particularly in terms of convergence and the shape of the cyclic region.

Equations and Parameters

It's a coupled ode equation of the system

where : w(t) = 1 - 0.5 * sin ( a*t + delta)

parameters : N = 4 , r_d = 0.6 , epchilon = 6 , alpha = 1.5 , beta = 3.5 , theta = 0.1 , delta = - pie / 2

Expected Output Expected graph

Actual Output My output graph

Code

Here’s my complete code:

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

# Parameters
N = 4
rd = 0.6
epsilon = 6
alpha = 1.5
beta = 3.5
theta = 0.5
a = 0.1
delta = -np.pi / 2

# Time-dependent w2 function
def w2(t):
    return 1 - 0.5 * np.sin(a * t + delta)

# Define the system of ODEs
def coupled_system(t, state):
    x, rc = state
    w = w2(t)  # Time-dependent w2
    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))

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

    return [dx_dt, drc_dt]

# Time span
t_span = (0, 50)
t_eval = np.linspace(*t_span, 1000)

# Increase initial conditions density
x0_vals = np.linspace(0.00, 1.00, 30)  # More starting points
rc0_vals = np.linspace(1.50, 3.50, 30)

# Store trajectories for different initial conditions
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-6,atol=1e-9)
        trajectories.append(sol)

# Plot trajectories with improved color coding and smoothing
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]

    # Color classification based on the final state
    if final_x < 0.001:  # Extinct cooperators
        color = '#1034a6'
    elif final_x > 0.80 and abs(final_rc - alpha) < 0.3 :  # Full cooperation
        color = '#FF1493'  # Deep Pink
    else:  # Cyclic behavior
        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()

Questions

Why is my phase portrait not matching the expected graph? Specifically:

1.Why do the pink trajectories converge to r_c = 2.0 - 2.5 instead of r_c = 1.5

2.Why is the cyclic region (orange trajectories) too elongated and not compact?