我试图用Python解决耦合的一阶ODEs系统。我是新手,但是Zombie Apocalypse example来自SciPy.org网站到目前为止帮助很大。在
在我的案例中,一个重要的区别是,用于“驱动”我的ODEs系统的输入数据在不同的时间点突然发生变化,我不确定如何最好地处理这一问题。下面的代码是我能想到的最简单的例子来说明我的问题。我很欣赏这个例子有一个直接的解析解,但是我的实际ODEs系统更复杂,这就是为什么我试图理解数值方法的基本原理。在
简化示例
考虑一个底部有洞的水桶(这种“线性水库”是许多水文模型的基本组成部分)。到铲斗的输入流量为R,孔的输出流量为Q。Q假设与桶中的水量成比例,V。比例常数通常写为,其中T是商店的“停留时间”。这是一个简单的形式颂歌
实际上,R是观测到的日降雨量总量的时间序列。在每天的内,假设降雨率是恒定的,但是在天之间,降雨率突然变化(即R是一个不连续的时间函数)。我试图理解这对解决我的颂歌的意义。在
策略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值:
这两种方法给出的答案相同,如下所示:
然而,第二种策略虽然在代码方面更紧凑,但它比第一种策略慢得多。我想这和R中的不连续性有关,导致odeint
的问题?在
我的问题
谢谢你!在
1.)是的
2.)是的
两者的原因:Runge-Kutta解算器期望ODE函数的可微阶数至少与解算器的阶数一样高。这是必要的,以便泰勒展开,给出预期的误差项存在。这意味着即使是1阶欧拉方法也需要一个可微的ODE函数。因此不允许跳跃,在1阶中可以容忍扭结,但在高阶解算器中则不能。在
对于具有自动步长调整的实现来说尤其如此。无论何时接近一个不满足微分顺序的点,解算器都会看到一个僵硬的系统,并将步长推向0,这将导致解算器的速度减慢。在
如果使用固定步长的解算器,且步长仅为1天的一小部分,则可以组合策略1和策略2。然后,在一天的转折点作为(隐式)新常量的重新开始点。在
相关问题 更多 >
编程相关推荐