ODEINT с несколькими параметрами (зависит от времени)

Question

Я пытаюсь решить один ODE первого порядка, используя ODEINT. Ниже приведен код. Я ожидаю получить 3 значения y для 3 временных точек. Проблема, с которой я борюсь, - это способность передавать n-е значение mt и nt для расчета dydt. Я думаю, что ODEINT передает все 3 значения mt и nt, вместо этого всего 0, 1 или 2, в зависимости от итерации. Из-за этого, я получаю эту ошибку:

RuntimeError: The size of the array returned by func (4) does not match the size of y0 (1).

Интересно, если я заменю начальное условие, которое (и должно быть) одно значение как: a0 = [2] * 4, код работает, но дает мне матрицу 4X4 в качестве решения, что кажется неправильным.

mt = np.array([3,7,4,2]) # Array of constants 
nt = np.array([5,1,9,3]) # Array of constants 
c1,c2,c3 = [-0.3,1.4,-0.5] # co-efficients 
para = [mt,nt] # Packing parameters 

#Test ODE function 
def test (y,t,extra): 
    m,n = extra 
    dydt = c1*c2*m - c1*y - c3*n 
    return dydt 

a0= [2] # Initial Condition 
tspan = range(len(mt)) # Define tspan 

#Solving the ODE 
yt= odeint(test, a0,tspan,args=(para,)) 
#Plotting the ODE 
plt.plot(tspan,yt,'g') 
plt.title('Multiple Parameters Test') 
plt.xlabel('Time') 
plt.ylabel('Magnitude') 

Дифференциальное уравнение порядка первого является:

dy/dt = c1*(c2*mt-y(t)) - c3*nt

Это уравнение представляет собой часть мышиного эндокринной системы, что я пытаюсь моделировать. Система аналогична системе с двумя резервуарами, где первый резервуар получает определенный гормон [по неизвестной скорости], но наш датчик обнаружит этот уровень (mt) с определенными временными интервалами (1 секунда). Затем этот бак поступает во второй резервуар, где уровень этого гормона (y) обнаруживается другим датчиком. Я обозначил уровни с использованием отдельных переменных, потому что датчики, которые обнаруживают уровни, независимы друг от друга и не откалиброваны друг с другом. «c2» можно рассматривать как коэффициент эффективности, который показывает корреляцию между двумя уровнями. Кроме того, передача этого гормона из резервуара 1 в резервуар 2 осуществляется с помощью диффузии. Этот гормон дополнительно потребляется биохимическим процессом (аналогично сливному клапану для второго резервуара). На данный момент неясно, какие параметры влияют на потребление; однако другой датчик может обнаруживать количество гормона (nt), потребляемого с определенным временным интервалом (1 секунда, в этом случае тоже).

Таким образом, mt и nt являются концентрациями/уровнями гормона в определенные моменты времени. хотя только 4-элементный код в коде, эти массивы намного больше в моем исследовании. Все датчики сообщают о концентрациях с интервалом в 1 секунду - следовательно, tspan состоит из временных точек, разделенных на 1 секунду.

Цель состоит в том, чтобы определить концентрацию этого гормона во втором резервуаре (y) математически, а затем оптимизировать значения этих коэффициентов на основе экспериментальных данных. Мне удалось передать эти массивы mt и nt в определенное ОДУ и решить проблему с помощью ODE45 в MATLAB. Я работал в этом RunTimeError, пытаясь реплицировать код в Python.

Answer 1

Как я уже упоминал в комментарии, если вы хотите смоделировать эту систему с использованием обычного дифференциального уравнения, вы должны сделать предположение о значениях m и n между временами выборки. Одной из возможных моделей является использование линейной интерполяции. Вот сценарий, который использует scipy.interpolate.interp1d для создания функций mfunc(t) и nfunc(t) на основе образцов mt и nt.

import numpy as np 
from scipy.integrate import odeint 
from scipy.interpolate import interp1d 
import matplotlib.pyplot as plt 

mt = np.array([3,7,4,2]) # Array of constants 
nt = np.array([5,1,9,3]) # Array of constants 

c1, c2, c3 = [-0.3, 1.4, -0.5] # co-efficients 

# Create linear interpolators for m(t) and n(t). 
sample_times = np.arange(len(mt)) 
mfunc = interp1d(sample_times, mt, bounds_error=False, fill_value="extrapolate") 
nfunc = interp1d(sample_times, nt, bounds_error=False, fill_value="extrapolate") 

# Test ODE function 
def test (y, t): 
    dydt = c1*c2*mfunc(t) - c1*y - c3*nfunc(t) 
    return dydt 

a0 = [2] # Initial Condition 
tspan = np.linspace(0, sample_times.max(), 8*len(sample_times)+1) 
#tspan = sample_times 

# Solving the ODE 
yt = odeint(test, a0, tspan) 

# Plotting the ODE 
plt.plot(tspan, yt, 'g') 
plt.title('Multiple Parameters Test') 
plt.xlabel('Time') 
plt.ylabel('Magnitude') 
plt.show() 

Вот сюжет, созданный сценарий:

Заметим, что вместо того, чтобы генерировать решение только на sample_times (т.е. время от времени 0, 1, 2 и 3), I установите tspan на более плотный набор точек. Это показывает поведение модели между временами выборки.