SciPy solve_ivp的解包含一阶常微分方程组的振荡

Question

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

其中，该示例的Tf被视为常数，并且给出Q（t）。Q（t）的绘图如下所示。用于创建时间vs Q打印的数据文件可在here处使用。

我的Python代码用于解决给定Q（t）（指定为^{cd1>}）的系统是：

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

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

Answer 1

正如已经在评论中所说的，建议插值Q。当试图用一个显式方法（如RK45（求解μivp的标准）求解刚性常微分方程组时，通常会出现振荡。因为Kurada的方法更像是一个隐式的方法

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

通过将雅可比矩阵提供给求解器，最终得到了求解常微分方程组的合理解。请参阅下面我的工作解决方案。在

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并不能改善结果。在