python - SciPysolve_ivp 的解决方案包含一阶 ODE 系统的振荡

Question

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

其中本例中的 Tf 被视为常数，并且给出了 Q(t)。 Q(t) 的图如下所示。用于创建时间 vs 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

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

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

给我: