python - 使用不连续的输入/强制数据求解 ODE

我正在尝试在 Python 中解决耦合的一阶 ODE 系统。我是新手，但 SciPy.org 的 Zombie Apocalypse example 到目前为止帮了大忙。

就我而言，一个重要的区别是，用于“驱动”我的 ODE 系统的输入数据在不同的时间点突然改变，我不确定如何最好地处理这个问题。下面的代码是我能想到的最简单的例子来说明我的问题。我很欣赏这个例子有一个简单的解析解，但我的实际 ODE 系统更复杂，这就是为什么我试图理解数值方法的基础知识。

简化示例

考虑一个底部有洞的水桶(这种“线性水库”是许多水文模型的基本构建块)。桶的输入流速为 R ，孔的输出流速为 Q 。 Q 假定与桶中的水体积 V 成正比。比例常数通常写为，其中 T 是商店的“驻留时间”。这给出了一个简单的 ODE 形式

实际上， R 是观测到的每日降雨总量的时间序列。在每一天中，降雨率被假定为常数，但在几天之间，降雨率突然变化(即 R 是 时间的不连续 函数)。我试图了解这对解决我的 ODE 的影响。

策略 1

最明显的策略(至少对我而言)是在每个降雨时间步长内分别应用 SciPy 的 odeint 函数。这意味着我可以将 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]

这两种方法都给出了相同的答案，如下所示:

然而第二种策略虽然在代码方面更紧凑，但它 比第一种慢很多 。我想这与 中的不连续性有关 R 导致 odeint 出现问题？

我的问题

策略 1 是这里最好的方法，还是有 更好的方法 ？

策略 2 是个坏主意吗？为什么它这么慢？

谢谢!

最佳答案

1.) 是的

2.) 是的

两者的原因:Runge-Kutta 求解器期望 ODE 函数的可微阶数至少与求解器的阶数一样高。这是必需的，以便给出预期误差项的泰勒展开式存在。这意味着即使是 1 阶欧拉方法也需要可微的 ODE 函数。因此不允许跳跃，可以容忍 1 阶的扭结，但不能容忍更高阶的求解器。

对于具有自动步长调整的实现尤其如此。每当接近不满足微分阶数的点时，求解器就会看到一个刚性系统并将步长推向 0，这会导致求解器变慢。

如果您使用具有固定步长且步长为 1 天的一小部分的求解器，则可以组合策略 1 和 2。然后，日轮的采样点作为(隐式)重新开始点，具有新的常数。

关于python - 使用不连续的输入/强制数据求解 ODE，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/32225472/