CODEWITHSUNDEEP
pythonmatplotlibscipydifferential-equations
dejsdukes
dejsdukes

Reputation: 121

Solving an ODE Numerically with SciPy

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

enter image description here

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)

Upvotes: 1

Views: 108

Answers (2)

Lutz Lehmann
Lutz Lehmann

Reputation: 26040

In the end you will want to solve the coupled system. This is not complicated, just make the state object a vector and return the derivatives in the correct order.

def fun(state,t, params):
    P, Q = state
    beta, mp, mq, a, b = params
    N = P+Q
    r0N, riN = b*N, a*N
    return [ (beta-mp-r0N)*P + riN*Q, r0N*P - (riN+mq)*Q ]

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

tsol = np.arange(0,6,0.1)
sol = odeint( fun, y0 = [ 1, 0],  t = tsol,  args = (params,))

Psol, Qsol = sol.T; plt.plot(tsol, Psol, tsol, Qsol)

enter image description here

Upvotes: 2

Cleb
Cleb

Reputation: 26039

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:

enter image description here

Upvotes: 1

Related Questions