基于Scipy的一维瞬态热传导方程求解需求（采用solve_ivp与odeint方法）

嘿各位大佬！我已经在MATLAB里搞定了一个一维瞬态热传导的模型，现在想转到Python用Scipy的solve_ivp和odeint两个方法来复现它，目前整理了问题定义和部分代码，想请教怎么完善剩下的部分～

问题详情 #

控制方程 #

我需要求解的一维瞬态热传导方程如下：

dT/dt = (k/(rho*cp))*d²T/dy²

几何与初始/边界条件 #

• 几何定义：y方向的一维直线结构，被划分为若干计算单元
• 初始条件（I.C.）：在t=0时刻，所有单元的温度均为25°C
• 边界条件（B.C.）：
  • a) 底部单元与热辊接触，满足：-k*dT/dy = h*(T_roller - T)
  • b) 上部单元暴露在大气中，满足：-k*dT/dy = h*(T - T_ambient)

已编写的solve_ivp代码片段 #

目前我写了部分初始化代码，还有一些参数没补全，也还没完成微分方程的定义和求解部分：

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

# Constants
Ny = 51            # Number of spatial grid points
Nt = 5000          # Number of time steps
velocity_x = 0.2   # Line speed in m/s

# Thermal properties of different layers
layer1_thickness = 0.001  # Thickness of layer 1 (meters)
layer1_k = 0.15           # Thermal conductivity of layer 1 (W/m-K)

rho_1 = ...  # 待填充的层1密度值
cp_1 = ...   # 待填充的层1比热容值

另外，我还想尝试用Scipy的odeint方法来实现同样的求解，希望能得到这两种方法的完整实现思路或代码示例～

备注：内容来源于stack exchange，提问作者Adya