python-2.7 - 具有多个参数的 ODEINT(时间相关)

Question

我正在尝试使用 ODEINT 求解单个一阶 ODE。以下是代码。我希望在 3 个时间点获得 3 个 y 值。我遇到的问题是能够传递 mt 和 nt 的第 n 个值来计算 dydt。我认为 ODEINT 传递了 mt 和 nt 的所有 3 个值，而不是仅传递第 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

这个方程代表了我正在尝试建模的小鼠内分泌系统的一部分。该系统类似于双 jar 系统，其中第一个 jar 接收特定激素 [以未知速率] 但我们的传感器将在特定时间间隔(1 秒)检测该水平 (mt) .然后，这个水箱会进入第二个水箱，在第二个水箱中，另一个传感器会检测到这种激素 (y) 的水平。我使用单独的变量标记级别，因为检测级别的传感器彼此独立，并且未相互校准。 “c2”可被视为显示两个级别之间相关性的系数。此外，这种激素从 jar 1 到 jar 2 的转移是扩散驱动的。这种激素被生化过程进一步消耗(类似于第二个水箱的排水阀)。目前，尚不清楚哪些参数会影响消耗；然而，另一个传感器可以检测在特定时间间隔(在本例中也是 1 秒)消耗的激素 (nt) 量。

因此，mt 和 nt 是特定时间点的激素浓度/水平。虽然代码中的长度只有 4 个元素，但在我的研究中这些数组要长得多。所有传感器均以 1 秒为间隔报告浓度 - 因此 tspan 由间隔 1 秒的时间点组成。

目标是通过数学方法确定第二个 jar 中这种激素的浓度 (y)，然后根据实验数据优化这些系数的值。我能够将这些数组 mt 和 nt 传递给定义的 ODE，并在 MATLAB 中使用 ODE45 求解，没有任何问题。在尝试用 Python 复制代码时，我遇到了这个 RunTimeError。

Answer 1

正如我在评论中提到的，如果您想使用常微分方程对该系统建模，您必须对 m 和 n 的值做出假设采样时间之间。一种可能的模型是使用线性插值。这是一个使用 scipy.interpolate.interp1d 的脚本根据示例 mt 和 nt 创建函数 mfunc(t) 和 nfunc(t)。

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)生成解决方案，而是将 tspan 设置为一组更密集的点.这显示了模型在采样时间之间的行为。