Численное решение ОДУ с помощью SciPy

Question

Я пытаюсь найти численное решение и, в конечном итоге, график, модель Гилленберга-Уэбба (модель роста раковых клеток). Эта модель выглядит так:

Где β - скорость размножения пролиферирующих клеток, µp - смертность пролиферирующих клеток, µq - смертность покоящихся клеток, а r0 и ri - функции (скорости перехода ) из N(t). Также N(t) = P(t)+Q(t).

Для моих целей здесь я определил r_0(N) = bN и r_i(N) = aN, чтобы упростить задачу.

Моя проблема в том, что когда я пытаюсь построить свое решение с помощью pyplot, я получаю

ValueError: x и y должны иметь одинаковое первое измерение

что само собой разумеется, но я не уверен, как это исправить, не ломая все остальное.

Мой код, который я сделал только для первого уравнения:

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

В конце концов вы захотите решить связанную систему. Это не сложно, просто сделайте объект состояния вектором и верните производные в правильном порядке.

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)

Answer 2

Вы в настоящее время планируете fun против tvec. То, что вы на самом деле хотите, это построить tvec против s1. Вам также нужно будет определить параметр a в fun; Я установил его на 1 в коде ниже:

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()

Это будет сюжет: