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

Question

我正在尝试求解耦合的一阶 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 耦合系统？

Answer 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()

给我：

Answer 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 不会改善结果。