我正在尝试用 Python 求解耦合一阶 ODE 系统。我对此很陌生，但是僵尸启示录示例 http://wiki.scipy.org/Cookbook/Zombie_Apocalypse_ODEINT到目前为止，SciPy.org 提供了很大的帮助。
在我的例子中，一个重要的区别是用于“驱动”我的 ODE 系统的输入数据发生了变化abruptly在不同的时间点，我不知道如何最好地处理这个问题。下面的代码是我能想到的用来说明我的问题的最简单的示例。我很欣赏这个例子有一个简单的解析解，但我的实际 ODE 系统更复杂，这就是为什么我试图理解数值方法的基础知识。
简化示例
Consider a bucket with a hole in the bottom (this kind of "linear reservoir" is the basic building block of many hydrological models). The input flow rate to the bucket is R and the output from the hole is Q. Q is assumed to be proportional to the volume of water in the bucket, V. The constant of proportionality is usually written as , where T is the "residence time" of the store. This gives a simple ODE of the form
事实上，R是观测到的每日降雨量时间序列。Within假设每天的降雨量是恒定的，但是between汇率突然变化的天数（即R is a 不连续的时间的函数）。我试图理解这对于解决我的常微分方程的影响。
策略1
最明显的策略（至少对我来说）是应用 SciPyodeint在每个降雨时间步长内分别起作用。这意味着我可以治疗R作为常数。像这样的东西：
import numpy as np, pandas as pd, matplotlib.pyplot as plt, seaborn as sn
from scipy.integrate import odeint
np.random.seed(seed=17)
def f(y, t, R_t):
""" Function to integrate.
"""
# Unpack parameters
Q_t = y[0]
# ODE to solve
dQ_dt = (R_t - Q_t)/T
return dQ_dt
# #############################################################################
# User input
T = 10 # Time constant (days)
Q0 = 0. # Initial condition for outflow rate (mm/day)
days = 300 # Number of days to simulate
# #############################################################################
# Create a fake daily time series for R
# Generale random values from uniform dist
df = pd.DataFrame({'R':np.random.uniform(low=0, high=5, size=days+20)},
index=range(days+20))
# Smooth with a moving window to make more sensible
df['R'] = pd.rolling_mean(df['R'], window=20)
# Chop off the NoData at the start due to moving window
df = df[20:].reset_index(drop=True)
# List to store results
Q_vals = []
# Vector of initial conditions
y0 = [Q0, ]
# Loop over each day in the R dataset
for step in range(days):
# We want to find the value of Q at the end of this time step
t = [0, 1]
# Get R for this step
R_t = float(df.ix[step])
# Solve the ODEs
soln = odeint(f, y0, t, args=(R_t,))
# Extract flow at end of step from soln
Q = float(soln[1])
# Append result
Q_vals.append(Q)
# Update initial condition for next step
y0 = [Q, ]
# Add results to df
df['Q'] = Q_vals
策略2
第二种方法是将所有内容简单地提供给odeint并让它处理不连续性。使用相同的参数和R值如上：
def f(y, t):
""" Function used integrate.
"""
# Unpack incremental values for S and D
Q_t = y[0]
# Get the value for R at this t
idx = df.index.get_loc(t, method='ffill')
R_t = float(df.ix[idx])
# ODE to solve
dQ_dt = (R_t - Q_t)/T
return dQ_dt
# Vector of initial parameter values
y0 = [Q0, ]
# Time grid
t = np.arange(0, days, 1)
# solve the ODEs
soln = odeint(f, y0, t)
# Add result to df
df['Q'] = soln[:, 0]
这两种方法给出了相同的答案，如下所示：
然而，第二种策略虽然代码更紧凑，但它much比第一个慢。我想这与不连续性有关R造成问题odeint?
我的问题
谢谢你！
1.) Yes
2.) Yes
两者的原因：Runge-Kutta 求解器期望 ODE 函数的可微性阶数至少与求解器的阶数一样高。这是必要的，以便给出预期误差项的泰勒展开式存在。这意味着即使是 1 阶欧拉方法也需要可微的 ODE 函数。因此不允许跳跃，在 1 阶中可以容忍扭结，但在更高阶的求解器中则不允许。
对于具有自动步长调整的实现尤其如此。每当接近微分阶数不满足的点时，求解器就会看到一个刚性系统，并将步长驱动到 0，这会导致求解器速度减慢。
如果您使用具有固定步长且步长为 1 天的一小部分的求解器，则可以组合策略 1 和 2。然后，日间循环的采样点将用作具有新常数的（隐式）重启点。