使用 scipy odeint 求解简单 ODE 在 0 处给出直线

Question

我正在尝试解决一个简单的 ODE： dN/dt = N*(rho(t)-beta)/lambda

Rho 是时间的函数，我使用 linspace 生成了它。该代码适用于其他方程，但不知何故在 0 处给出了一条平坦的直线。（您可以在图中看到它）。有关如何纠正它的任何指导方针？

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

def model2(N, t, rho):
    beta_val = 0.0065
    lambda_val = 0.00002
    k = (rho - beta_val) / lambda_val
    dNdt = k*N
    print(rho)
    return dNdt

# initial condition
N0 = [0]

# number of time points
n = 200

# time points
t = np.linspace(0,200,n)

rho = np.linspace(6,9,n)
#rho =np.array([6,6.1,6.2,6.3,6.4,6.5,6.6,6.7,6.8,6.9,7.0,7.1,7.2,7.3,7.4,7.5,7.6,7.7,7.8,7.9])  # Array of constants

# store solution
NSol = np.empty_like(t)
# record initial conditions
NSol[0] = N0[0]

# solve ODE
for i in range(1,n):
    # span for next time step
    tspan = [t[i-1],t[i]]
    # solve for next step
    N = odeint(model2,N0,tspan,args=(rho[i],))
    print(N)
    # store solution for plotting
    NSol[i] = N[0][0]
    # next initial condition
    #z0 = N0[0]

# plot results
plt.plot(t,rho,'g:',label='rho(t)')
plt.plot(t,NSol,'b-',label='NSol(t)')
plt.ylabel('values')
plt.xlabel('time')
plt.legend(loc='best')
plt.show()

This is the graph I get after running this code

Answer 1

我修改了您的代码（和系数）以使其正常工作。 当系数也依赖于t时，它们必须是导数函数调用的python函数：

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

# Define
def model2(N, t, rho):
    beta_val = 0.0065
    lambda_val = 0.02
    k = ( rho(t) - beta_val )/lambda_val
    dNdt = k*N
    return dNdt

def rho(t):
    return .001 + .003/20*t


# Solve
tspan = np.linspace(0, 20, 10)
N0 = .01
N = odeint(model2, N0 , tspan, args=(rho,))

# Plot
plt.plot(tspan, N, label='NS;ol(t)');
plt.ylabel('N');
plt.xlabel('time'); plt.legend(loc='best');