graphing nonlinear systems of ODEs in python

Question

the model Im trying to graph this model however I am just not sure whether my code is leading me to a reasonable solution:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
K = 17.3
epsilon = 4.9
r = 2.7
mu = 0.9
gamma = 3.2 
alpha = 4 / (K - epsilon) ** 2 
def mosquito_odes(t, Y, alpha, r, K, epsilon, mu, gamma):
 N, M = Y 
 
 dNdt = alpha * r * (N**2 / (N + M + 1e-6)) * (K - N) * (N - epsilon) - mu * N
 dMdt = gamma * N / (1 + N) - mu * M
 
 return [dNdt, dMdt]
N0 = K/1.5 
M0 = 5.6667 
Y0 = [N0, M0]
T = 5 
t_eval = np.linspace(0, T, 100) 
sol_ode = solve_ivp(
 mosquito_odes,
 [0, T],
 Y0,
 args=(alpha, r, K, epsilon, mu, gamma),
 method="RK45",
 t_eval=t_eval
)
t_vals = sol_ode.t
N_vals, M_vals = sol_ode.y
plt.figure(figsize=(8, 5))
plt.plot(t_vals, N_vals / K, label="Wild Population (N)", linestyle="-", color="b")
plt.plot(t_vals, M_vals / K, label="Sterile Population (M)", linestyle="--", color="r")
plt.axhline(1.9 / K, linestyle="dashed", color="brown", label="Threshold γc=1.9") # Show threshold
plt.xlabel("Time")
plt.ylabel("Fraction of Population")
plt.title("Non-Spatial Mosquito Population Dynamics")
plt.legend()
plt.grid()
plt.show()

the outcome is not what I expect (for these parameters I expect both ODEs to go to 0. however it seems that only one is.

Answer 1

TL;DR.

It is always worthy to draw a phase diagram whenever you can:

enter image description here

We see there is three fixed states: two are sinks, one is a source.

The reason why both curves does not goes to zero is because your initial condition leads to a path that is not attracted towards the (0,0) fixed point.

How to proceed?

With your model and setup:

import numpy as np
from scipy import integrate, optimize
import matplotlib.pyplot as plt
K = 17.3
epsilon = 4.9
r = 2.7
mu = 0.9
gamma = 3.2 
alpha = 4 / (K - epsilon) ** 2 
def model(t, y, alpha, r, K, epsilon, mu, gamma):
 return np.array([
 gamma * y[1] / (1 + y[1]) - mu * y[0],
 alpha * r * (y[1] ** 2 / (y[1] + y[0])) * (K - y[1]) * (y[1] - epsilon) - mu * y[1]
 ])
M0 = 5.6667
N0 = K / 1.5

We can draw such a diagram, first we create grids:

Mlin = np.linspace(-1, 20, 200)
Nlin = np.linspace(-1, 20, 200)
M, N = np.meshgrid(Mlin, Nlin)
t = np.linspace(2, 10, 200)

Then we compute vector field over the grid:

p0 = (alpha, r, K, epsilon, mu, gamma)
U, V = model(t, [M, N], *p0)
C = np.sqrt(U ** 2 + V ** 2)

We compute roots of the field:

def proxy(x, *args):
 return model(None, x, *args)
sol0 = optimize.root(proxy, x0=[0.025, 0.01], args=p0, tol=1e-3)
sol1 = optimize.root(proxy, x0=[5, 7.5], args=p0)
sol2 = optimize.root(proxy, x0=[2.5, 17.5], args=p0)

We solve some IVP for specific inital values:

y0 = (M0, N0)
y1 = (4, 6)
y2 = (10, 2.5)
ivp0 = integrate.solve_ivp(model, [t.min(), t.max()], args=p0, y0=y0, t_eval=t)
ivp1 = integrate.solve_ivp(model, [t.min(), t.max()], args=p0, y0=y1, t_eval=t)
ivp2 = integrate.solve_ivp(model, [t.min(), t.max()], args=p0, y0=y2, t_eval=t)

And we draw results:

prop_cycle = plt.rcParams['axes.prop_cycle']
colors = prop_cycle.by_key()['color']
fig, axe = plt.subplots()
for ivp, color in zip([ivp0, ivp1, ivp2], colors[:3]):
 axe.plot(t, ivp.y[0], "-", color=color)
 axe.plot(t, ivp.y[1], "--", color=color)
axe.grid()
fig, axe = plt.subplots()
axe.streamplot(M, N, U, V, color="grey")
for sol in [sol0, sol1, sol2]:
 axe.scatter(*sol.x, marker="o", color="blue")
for y, ivp, color in zip([y0, y1, y2], [ivp0, ivp1, ivp2], colors[:3]):
 axe.scatter(*y, color=color)
 axe.plot(*ivp.y, color=color)
axe.grid()

We can confirm we can reach the (0,0) fixed point by choosing the ad hoc initial value:

enter image description here

Bonus

Also notice, the (0,0) fixed point is a bit special and prevent optimize.root to converge. Situation around its neighborhood is about:

enter image description here

Answer 2

Your system is nonlinear and, consequently, very dependent on the initial boundary conditions. If you set N0 = K/2.5 for example, rather than N0 = K/1.5 you will indeed get a solution with both populations tending to 0.

You should investigate the possible steady-state solutions to your problem; i.e. those with the RHS=0 in both equations. The second equation gives M in terms of N; if you substitute that in the RHS of the first equation you will get (after removing the trivial solution M=N=0) a cubic equation for N. Solving this gives possible non-zero steady states of N, and thence of M.

I had a go at solving the steady-state equation below. You will see that, for your coefficients, it finds the non-trivial roots. In particular, with your current values for N0 and M0 your timestepping allows a particular steady solution

N/K, M/K = 0.9184 0.1934

which is the one shown on your graph (at the bottom of this post).

Root-finder:

import numpy as np
import matplotlib.pyplot as plt
from numpy.polynomial import Polynomial
K = 17.3
epsilon = 4.9
r = 2.7
mu = 0.9
gamma = 3.2 
alpha = 4 / (K - epsilon) ** 2 
ar = alpha * r
c0 = -ar * K * epsilon - mu - gamma
c1 = ar * ( K + epsilon - K * epsilon ) - mu
c2 = ar * ( K + epsilon - 1 )
c3 = -ar
p = Polynomial( [ c0, c1, c2, c3 ] )
for N in p.roots():
 M = ( gamma / mu ) * N / ( 1 + N )
 print( 'N/K, M/K =', N/K, M/K )

Output of root-finder

N/K, M/K = -0.07815057938043678 0.7893885266055993
N/K, M/K = 0.38517782284686214 0.17870522690567683
N/K, M/K = 0.918406282545136 0.19335395971294628

The original post gives time-dependent graph as below. It is heading toward a steady state corresponding to the last of these roots. enter image description here