用solve_ivp解微分方程时返回的奇数结果

Question

我正在尝试解决一系列微分方程，如下面的代码所示。但是当我模拟这些微分方程时，我得到了一个不准确的结果。在这里，我有4个变量，我首先重塑了我的向量，然后将这些值连接起来，希望每个变量的向量长度相等。这样，我将初始条件设置为变量的数量*N (N表示我想要模拟的单元数量)。代码中所示的for循环用于模拟细胞群体，而不是单个细胞，因此我将for循环的值从1设置为N。然而，我遇到的问题是，当我返回结果并绘制变量(X，Y，Z，V)时，我的X，Z和V值为0，而实际上我应该看到振荡。此外，表示平均场的F值也应该显示振荡，但当我使用print语句时，F值为负值(对于这些方程，我不会得到任何负值)。我已经检查过我的代码很多次了，我不确定为什么当所有的方程都正确列出时，我没有得到任何振荡。我认为我没有正确地将变量传递到求解器中，或者如果我需要使用不同的方法来获得相等长度的向量。有没有人知道为什么我没有得到X，Z，V的振荡，却得到了Y的振荡？我该如何纠正这个错误呢？任何帮助都是非常感谢的(代码和图在下面)。

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline     
from scipy.integrate import solve_ivp

def testing(t,y):
    X,Y,Z,V = y.reshape([4,-1])
    N = 1000; v1 = 6.8355; K1 = 2.7266; n = 5.6645; v2 = 8.4297; K2 = 0.2910 
    k3 = 0.1177; v4 = 1.0841; K4 = 8.1343; k5 = 0.3352; v6 = 4.6645 
    K6 = 9.9849; k7 = 0.2282; v8 = 3.5216; K8 = 7.4519; vc = 6.7924
    Kc = 4.8283; K = 1; L = 0 
    dXdt = np.zeros(N)
    dYdt = np.zeros(N)
    dZdt = np.zeros(N)
    dVdt = np.zeros(N)
    for i in range(1,N+1):
        F = 1/N * sum([V[i]])
        dXdt[i] = (v1*(K1**n))/((K1**n)+(Z[i]**n))-((v2*X[i])/(K2+X[i])) + ((vc*K*F)/(Kc+(K*F))) + L
        dYdt[i] = (k3*X[i]) - ((v4*Y[i])/(K4+Y[i]))
        dZdt[i] = (k5*Y[i]) - ((v6*Z[i])/(K6+Z[i]))
        dVdt[i] = (k7*X[i]) - ((v8*V[i])/(K8+V[i]))
        return np.concatenate([dXdt,dYdt,dZdt,dVdt])
t_span = (0, 96)
t = np.linspace(t_span[0], t_span[1], 3000) 
N = 1000
y0 = np.zeros(4*N)
soln = solve_ivp(testing, t_span, y0, t_eval=t)
t = soln.t; X = soln.y[0]; Y = soln.y[1]; Z = soln.y[2]; V = soln.y[3]

结果图：CLICK TO SEE PLOT

Answer 1

ODE函数testing中有几个问题

• 返回语句的缩进是错误的，它在循环中，所以在第一次运行时触发它

• F = 1/N * sum([V[i]])显然是错误的，您希望将F= 1/N * sum(V)放在循环之外，以便对所有i。

计算一次。

• N应该是接收到的数组的实际长度，以便平均值计算实际计算平均值N=len(V).

• 求解器尝试以某个负值计算Z的分数幂，可能是在RK45方法的某个中间阶段。通过绝对值来保证它的安全性，使用Z*abs(Z)**(n-1).

来保持一些标志

在您对解决方案的解释中，索引也是不正确的。

• X = soln.y[0]; Y = soln.y[1]; Z = soln.y[2]; V = soln.y[3]与状态组件布局不对应。第一个单元格的组件是X,Y,Z,V = soln.y[0::N].

或X = soln.y[0]; Y = soln.y[N]; Z = soln.y[2*N]; V = soln.y[3*N]

使用numpy数组运算的更正版本是

def testing(t,y):
    X,Y,Z,V = y.reshape([4,-1])
    N = len(X); v1 = 6.8355; K1 = 2.7266; n = 5.6645; v2 = 8.4297; K2 = 0.2910 
    k3 = 0.1177; v4 = 1.0841; K4 = 8.1343; k5 = 0.3352; v6 = 4.6645 
    K6 = 9.9849; k7 = 0.2282; v8 = 3.5216; K8 = 7.4519; vc = 6.7924
    Kc = 4.8283; K = 1; L = 0 
    F = 1/N * sum(V)
    dXdt = (v1*(K1**n))/((K1**n)+(Z*abs(Z)**(n-1)))-((v2*X)/(K2+X)) + ((vc*K*F)/(Kc+(K*F))) + L
    dYdt = (k3*X) - ((v4*Y)/(K4+Y))
    dZdt = (k5*Y) - ((v6*Z)/(K6+Z))
    dVdt = (k7*X) - ((v8*V)/(K8+V))
    return np.concatenate([dXdt,dYdt,dZdt,dVdt])

并给出了，即使是对于N=10，良好的振荡

​

​