求解具有间断输入/强迫d的常微分方程

2024-03-29 10:10:41 发布

您现在位置:Python中文网/ 问答频道 /正文

我试图用Python解决耦合的一阶ODEs系统。我是新手,但是Zombie Apocalypse example来自SciPy.org网站到目前为止帮助很大。在

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

简化示例

考虑一个底部有洞的水桶(这种“线性水库”是许多水文模型的基本组成部分)。到铲斗的输入流量为R,孔的输出流量为QQ假设与桶中的水量成比例,V。比例常数通常写为formula1,其中T是商店的“停留时间”。这是一个简单的形式颂歌

fromula2

Simple linear reservoir

实际上,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值:

^{pr2}$

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

Result

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

我的问题

  1. 战略1是最好的方法,还是有更好的方法?在
  2. 战略2是个坏主意吗?为什么这么慢?在

谢谢你!在


Tags: oftheto方法dffor系统as
1条回答
网友
1楼 · 发布于 2024-03-29 10:10:41

1.)是的

2.)是的

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

对于具有自动步长调整的实现来说尤其如此。无论何时接近一个不满足微分顺序的点,解算器都会看到一个僵硬的系统,并将步长推向0,这将导致解算器的速度减慢。在

如果使用固定步长的解算器,且步长仅为1天的一小部分,则可以组合策略1和策略2。然后,在一天的转折点作为(隐式)新常量的重新开始点。在

相关问题 更多 >