Solving an ODE Numerically with SciPy

Question

I'm trying to find a numerical solution and eventually graph, the Gyllenberg-Webb model (cancer cell growth model). This model looks like:

Where β is the reproduction rate of proliferating cells, µp is the death rate of proliferating cells, µq is the death rate of quiescent cells, and r0 and ri are functions (transition rates) of N(t). Also N(t) = P(t)+Q(t).

For my purposes here I defined r_0(N) = bN and r_i(N) = aN to make things more simple.

My problem is when I try and plot my solution with pyplot I get

ValueError: x and y must have same first dimension

which I guess is self-explanatory, but I'm not sure how to go about fixing it without breaking everything else.

My code, which I've done only for the first equation so far, is:

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

def fun(P,t, params):
    beta, mp,b,N,Q = params
    return(beta-mp-(b*N))*P+(a*N)*Q

params = (0.5,0.6,0.7,0.8,0.9)

tvec = np.arange(0,6,0.1)
s1 = scipy.integrate.odeint(
    fun,
    y0 = 1,
    t = tvec,
    args = (params,))

#print(s1)

plt.plot(fun,tvec)

Answer 1

You are currently plotting fun vs. tvec. What you actually want is to plot tvec vs s1. You will also have to define the parameter a in fun; I set it to 1 in the code below:

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


def fun(P, t, params):
    beta, mp, b, N, Q = params
    return (beta-mp-(b*N))*P + (1.0 * N)*Q

params = (0.5, 0.6, 0.7, 0.8, 0.9)

tvec = np.arange(0, 6, 0.1)
s1 = scipy.integrate.odeint(
    fun,
    y0=1.,
    t=tvec,
    args=(params,))

plt.plot(tvec, s1)
plt.show()

This will plot: