[Python] Phase Portrait of Coupled ODEs Not Matching Expected Graph in Python (SciPy)

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?

Continue reading...