SciPy solve_ivp 的解决方案包含一阶 ODE 系统的振荡

Posted 2023-03-12

【中文标题】SciPy solve_ivp 的解决方案包含一阶 ODE 系统的振荡

【英文标题】：Solution from SciPy solve_ivp contains oscillations for a system of first-order ODEs

【发布时间】：2019-12-23 06:13:01

【问题描述】：

我正在尝试求解耦合的一阶 ODE 系统：

这里的 Tf 被认为是常数，Q(t) 是给定的。 Q(t) 的图如下所示。用于创建时间与 Q 图的数据文件可在here 获得。

对于给定的 Q(t)（指定为 qheat），我用于解决此系统的 Python 代码是：

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

# Data

time, qheat = np.loadtxt('timeq.txt', unpack=True)

# Calculate Temperatures    

def tc_dt(t, tc, ts, q):
    rc = 1.94
    cc = 62.7
    return ((ts - tc) / (rc * cc)) + q / cc

def ts_dt(t, tc, ts):
    rc = 1.94
    ru = 3.08
    cs = 4.5
    tf = 298.15
    return ((tf - ts) / (ru * cs)) - ((ts - tc) / (rc * cs))

def func(t, y):
    idx = np.abs(time - t).argmin()
    q = qheat[idx]

    tcdt = tc_dt(t, y[0], y[1], q)
    tsdt = ts_dt(t, y[0], y[1])
    return tcdt, tsdt

t0 = time[0]
tf = time[-1]
sol = solve_ivp(func, (t0, tf), (298.15, 298.15), t_eval=time)

# Plot

fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0], label='tc')
ax.plot(sol.t, sol.y[1], label='ts')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Temperature [K]')
ax.legend(loc='best')

plt.show()

这会产生如下所示的图，但不幸的是，结果中出现了几个振荡。有没有更好的方法来解决这个 ODE 耦合系统？

【问题讨论】：

您对q 的近似值可能是问题所在。我建议为此创建一个interpolation function。

我没有在我的示例中展示它，但我尝试使用 interp1d 表示 q 但它对结果中的波动没有帮助。

【参考方案1】：

就像在 cmets 中已经说过的那样，建议对 Q 进行插值。当尝试使用 RK45（solve_ivp 的标准）等显式方法求解刚性 ODE 系统时，通常会发生振荡。由于您的 ODE 系统似乎是一个僵硬的系统，因此进一步建议使用像“Radau”这样的隐式 Runge-Kutta 方法：

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d

# Data
time, qheat = np.loadtxt('timeq.txt', unpack=True)

# Interpolate Q
Q = interp1d(time, qheat) # linear spline

# Calculate Temperatures

def tc_dt(t, tc, ts, q):
    rc = 1.94
    cc = 62.7
    return ((ts - tc) / (rc * cc)) + q / cc

def ts_dt(t, tc, ts):
    rc = 1.94
    ru = 3.08
    cs = 4.5
    tf = 298.15
    return ((tf - ts) / (ru * cs)) - ((ts - tc) / (rc * cs))

def func(t, y):
    idx = np.abs(time - t).argmin()

    tcdt = tc_dt(t, y[0], y[1], Q(t))
    tsdt = ts_dt(t, y[0], y[1])
    return tcdt, tsdt

t0 = time[0]
tf = time[-1]
# Note the passed values for rtol and atol.
sol = solve_ivp(func, (t0, tf), (298.15, 298.15), method="Radau", t_eval=time, atol=1e-8, rtol=1e-8)

# Plot

fig, ax = plt.subplots()
ax.plot(sol.t, sol.y[0], label='tc')
ax.plot(sol.t, sol.y[1], label='ts')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Temperature [K]')
ax.legend(loc='best')

plt.show()

给我：

【讨论】：

您的答案错过了应该出现的高峰。 Q函数有10个峰

@MatthewFlamm 我编辑了我的答案。我认为没有必要提供雅可比，而是降低误差估计的绝对容差（atol）和相对容差（rtol）会产生相同的合理解决方案。对于面临同样问题且使用更复杂的 ode-rhs 并且不想手动提供 jacobian 的用户来说，了解这一点可能很有用。

我同意不需要提供 jacobian（尽管推荐）。除了隐式方法之外，降低公差或修改 max_step 大小应该可以捕获所有峰值。您可能想说明为什么不使用默认公差，答案中同时包含这两种情况实际上非常有用。

【参考方案2】：

通过向求解器提供雅可比矩阵，我终于得到了 ODE 系统的合理解。请参阅下面的我的工作解决方案。

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import interp1d

# Data

time, qheat = np.loadtxt('timeq.txt', unpack=True)

# Calculate Temperatures

interp_qheat = interp1d(time, qheat)

def tc_dt(t, tc, ts, q):
    """
    dTc/dt = (Ts-Tc)/(Rc*Cc) + Q/Cc
    """
    rc = 1.94
    cc = 62.7
    return ((ts - tc) / (rc * cc)) + q / cc

def ts_dt(t, tc, ts):
    """
    dTs/dt = (Tf-Ts)/(Ru*Cs) - (Ts-Tc)/(Rc*Cs)
    """
    rc = 1.94
    ru = 3.08
    cs = 4.5
    tf = 298.15
    return ((tf - ts) / (ru * cs)) - ((ts - tc) / (rc * cs))

def jacobian(t, y):
    """
    Given the following system of ODEs

    dTc/dt = (Ts-Tc)/(Rc*Cc) + Q/Cc
    dTs/dt = (Tf-Ts)/(Ru*Cs) - (Ts-Tc)/(Rc*Cs)

    determine the Jacobian matrix of the right-hand side as

    Jacobian matrix = [df1/dTc, df2/dTc]
                      [df1/dTs, df2/dTs]

    where

    f1 = (Ts-Tc)/(Rc*Cc) + Q/Cc
    f2 = (Tf-Ts)/(Ru*Cs) - (Ts-Tc)/(Rc*Cs)
    """
    cc = 62.7
    cs = 4.5
    rc = 1.94
    ru = 3.08
    jc = np.array([
        [-1 / (cc * rc), 1 / (cs * rc)],
        [1 / (cc * rc), -1 / (cs * ru) - 1 / (cs * rc)]
    ])
    return jc

def func(t, y):
    """
    Right-hand side of the system of ODEs.
    """
    q = interp_qheat(t)
    tcdt = tc_dt(t, y[0], y[1], q)
    tsdt = ts_dt(t, y[0], y[1])
    return tcdt, tsdt

t0 = time[0]
tf = time[-1]
sol = solve_ivp(func, (t0, tf), (298.15, 298.15), method='BDF', t_eval=time, jac=jacobian)

# Plot

fig, ax = plt.subplots(tight_layout=True)
ax.plot(sol.t, sol.y[0], label='tc')
ax.plot(sol.t, sol.y[1], label='ts')
ax.set_xlabel('Time [s]')
ax.set_ylabel('Temperature [K]')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), frameon=False)

plt.show()

生成的图如下所示。

插入 Q 的唯一好处是通过删除 main 函数中的 argmin() 来加快代码的执行。否则，插值 Q 不会改善结果。

【讨论】：

通过切换到 BDF（一种隐式求解器）而不是默认求解器 RK45（一种显式求解器），可能会改进解决方案。

您可能还想考虑限制最大步长，即使您在这里没问题，在某些情况下您也可以轻松跨过尖峰。