如何将给定Td、Ta离散递推系统转换为ODE并使用Python求解

一、离散转ODE推导 #

你给出的离散递推本质是向前欧拉差分格式，当时间步长Δt趋近于0时，差分可以直接转化为微分形式，推导过程如下：

1. 对于Td的递推式：

Td(t+Δt) = Td(t) + Δt * S/(MC) * (H(Ta(t)-Td(t)) - L*P)
两边减Td(t)后除以Δt，取Δt→0的极限，得到微分方程：
dTd/dt = (S/(M*C)) * (H*(Ta - Td) - L*P)

2. 对于Ta的递推式：

Ta(t+Δt) = Ta(t) + Δt * HS/(AP) * (Td(t+Δt) - Ta(t))
当Δt趋近于0时，Td(t+Δt)与Td(t)的差值为高阶无穷小，可以忽略，同理转化为：
dTa/dt = (H*S/(A*P)) * (Td - Ta)

注意：你当前给出的系统中S为常数，没有对应的演化规则，因此S随时间的变化规律为恒定值。如果你需要S随时间变化，需要补充S的表达式或对应的微分方程，加入到方程组中即可。

二、Python求解代码 #

我们使用scipy.integrate.solve_ivp求解上述ODE方程组，示例代码如下：
首先安装依赖库：

pip install numpy scipy matplotlib

求解代码：

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

# 定义常数（可根据你的实际参数替换）
M = 1.0       # 示例值
C = 4186.0    # 示例值，水的比热容
H = 10.0      # 示例值
L = 0.01      # 示例值
P = 101325.0  # 示例值，标准大气压
A = 0.1       # 示例值
S = 1000.0    # 示例值，固定常数

# 定义ODE方程组
def ode_system(t, y):
    Td, Ta = y
    dTd_dt = (S / (M * C)) * (H * (Ta - Td) - L * P)
    dTa_dt = (H * S / (A * P)) * (Td - Ta)
    return [dTd_dt, dTa_dt]

# 初始条件
Td0 = 20.0  # Td初始值，替换为你的实际值
Ta0 = 25.0  # Ta初始值，替换为你的实际值
y0 = [Td0, Ta0]

# 求解时间范围（可调整）
t_start = 0
t_end = 100
t_eval = np.linspace(t_start, t_end, 1000)

# 求解ODE
sol = solve_ivp(ode_system, [t_start, t_end], y0, t_eval=t_eval)

# 提取结果
t = sol.t
Td_res = sol.y[0]
Ta_res = sol.y[1]
S_res = np.full_like(t, S)  # S为常数

# 绘图查看结果
plt.rcParams['font.sans-serif'] = ['SimHei']  # 解决中文显示问题
plt.figure(figsize=(10,6))
plt.plot(t, Td_res, label='Td')
plt.plot(t, Ta_res, label='Ta')
plt.plot(t, S_res, label='S')
plt.xlabel('时间')
plt.ylabel('值')
plt.legend()
plt.grid(True)
plt.show()

三、结果说明 #

• 求解得到的Td_res、Ta_res、S_res分别对应三个参数在各个时间点的取值
• 如果你需要S随时间变化，只需要修改ode_system函数，加入S的微分方程，同时将S加入状态变量y中即可

内容的提问来源于stack exchange，提问作者Shay