PAYATHON 2026

移动热源模拟中，如何处理时间与空间数组长度不同的插值问题？

数据阿明

时间数组和空间数组不需要等长。100 个空间点与 9000 个时间点可以共同构成温度解。当前代码的问题是把时间到位置的映射、空间索引查找和温度快照提取混在了一起。

wf_pos(t) 也需要修改。solve_ivp 调用 ODE 函数时传入的是当前时刻的标量 t，而 interp1d(t, position) 需要成对的采样数据。这里很可能触发了 TypeError: object of type 'float' has no len()，但仍应以完整 traceback 为准。

1. 先让位置函数返回当前位置

如果热源轨迹已有解析表达式，可以直接计算位置。保留原有表达式时，可以先把函数改成：

def wf_pos(t):
    return np.sqrt(t * dt)

这个函数输入一个时刻，返回对应的位置。原代码返回的是插值函数对象，后续却把它当作位置传给 lindex_func。

还需要统一 t 的含义：

  • 如果 t 是采样序号，t * dt 可以表示经过的时间。
  • 如果 t 已经是物理时间，就不应再次把它当作采样序号进行换算。

solve_ivp 中的 t 使用 (t_start, t_end) 定义的时间坐标。如果绘图时使用 t[...] / dt，积分时却直接使用 t，同一个位置函数收到的输入就有不同含义。

位置还必须具有长度单位。如果输入时间的单位是秒，平方根轨迹通常需要一个具有相应单位的比例系数。原材料没有提供该系数，因此无法确认原公式的量纲，也无法确认热源是否始终位于杆内。

2. 只有离散轨迹数据才需要构造时间插值

如果热源位置来自离散采样，应准备两组一一对应的数据：time_samples 和 position_samples。它们需要等长，与空间网格 x_points 的长度无关。

from scipy.interpolate import interp1d

# time_samples 和 position_samples 是已知的一维轨迹数据。
# 两者长度相同，time_samples 严格递增。
position_interp = interp1d(
    time_samples,
    position_samples,
    kind='linear',
    bounds_error=True,
)


def wf_pos(t):
    return float(position_interp(t))

插值器应在积分前创建一次，ODE 中只需查询当前位置。轨迹数据需要覆盖整个 (t_start, t_end)，因为积分器会在内部选择时间点评估导数，不能只考虑 t_eval 中的输出时刻。

下面的写法不适用于当前问题：

interp1d(time_samples, x_points)

x_points 是杆上的固定空间坐标，不是各个采样时刻对应的热源位置。即使强行把两组数组改成相同长度，也不能建立所需的物理关系。

3. 在 ODE 中区分空间网格和热源位置

原代码中的 x = wf_pos(t) 会覆盖参数 x 的局部名称。建议使用含义明确的变量名：

source_position = wf_pos(t)
lindex = lindex_func(source_position)
Q_val = Qi_func(lindex)

还需要检查 Qi_func 的返回值。

对于长度为 Nx 的温度状态 T，移动局部热源通常需要一个形状为 (Nx,) 的源项数组，用来描述当前时刻各空间点的发热量。如果 Q_val 是标量，下面的操作会通过广播给所有空间点加上相同的值：

Tp += Q_val

这未必符合移动热源的设定。由于没有提供 lindex_func 和 Qi_func，目前只能确认需要检查源项的形状和空间分布，不能直接给出完整的替代实现。

热源落在两个网格点之间是正常情况。如果源项的离散方式要求把热量分配给相邻节点，可以根据位置计算权重；如果要求取最近节点，则可以查找最近索引。具体采用哪种方式，应与热源宽度、源项单位和离散模型保持一致。

4. 温度快照应从解矩阵中按时间取列

对于当前默认的调用方式，如果初始温度 T0 有 Nx 个分量，成功积分后：

sol.t.shape == (输出时间点数,)
sol.y.shape == (Nx, 输出时间点数)

某个输出时刻对应的完整温度分布可以这样取出：

j = np.argmin(np.abs(sol.t - snapshot_time))
ax.plot(x_points, sol.y[:, j])

这种方式取到的是最接近 snapshot_time 的已存储时刻，绘图时应使用 sol.t[j] 标注实际时间。

如果希望在指定时间准确查询积分器的连续近似解，可以启用 dense_output=True：

sol = integrate.solve_ivp(
    ode,
    (t_start, t_end),
    T0,
    args=args,
    t_eval=t_eval,
    rtol=1e-6,
    atol=1e-6,
    method='RK45',
    dense_output=True,
)

if not sol.success:
    raise RuntimeError(sol.message)

下面的绘图示例假定时间单位是秒，并且积分区间覆盖从起始时刻算起的第 3、6、9、12、15 分钟：

import matplotlib.pyplot as plt
import numpy as np

snapshot_times = t_start + 180.0 * np.arange(1, 6)

if np.any(snapshot_times < sol.t[0]) or np.any(snapshot_times > sol.t[-1]):
    raise ValueError('快照时间超出已求解区间')

fig, axes = plt.subplots(5, 1, figsize=(7, 12))

for ax, snapshot_time in zip(axes, snapshot_times):
    temperature = sol.sol(snapshot_time)
    ax.plot(x_points, temperature)
    ax.set_xlabel('Position (m)')
    ax.set_ylabel('Temperature')
    ax.set_title(f't = {snapshot_time:g} s')
    ax.grid(True)

fig.tight_layout()
plt.show()

每个快照的横坐标和纵坐标都有 Nx 个元素，无需截取一段时间数组来匹配空间坐标。

如果要绘制热源经过某个位置时的温度，或某个时间段的平均温度，则需要先明确对应的物理量，再采用相应的数据提取方式。原材料中 u 的定义不完整，不能直接判断应该怎样切片。

5. 截取空间区间时，不要依赖浮点数精确相等

以下写法要求热源位置恰好等于某个网格坐标：

np.where(x_points == wf_start)

任意连续轨迹通常不满足这个条件。即使数学上相等，浮点计算也可能带来微小差异。

如果只想显示两个位置之间的温度，可以直接建立空间掩码：

position_start = wf_pos(segment_start)
position_end = wf_pos(segment_end)
left, right = sorted((position_start, position_end))

mask = (x_points >= left) & (x_points <= right)
temperature = sol.sol(snapshot_time)

ax.plot(x_points[mask], temperature[mask])

这样选出的是两端位置之间已有的网格点，不能保证曲线包含精确的区间端点。如果区间窄于网格间距，选出的点也可能为空。

如果需要绘制任意空间坐标处的温度，应对同一时刻的空间温度分布进行插值：

# x_query 必须位于 x_points 覆盖的范围内。
temperature_query = np.interp(
    x_query,
    x_points,
    temperature,
)

参与配对的是 x_points 和 temperature，两者长度相同；查询位置 x_query 可以有不同数量的元素。

对于往返移动的轨迹，仅凭一个时间段的起点和终点，不能确定热源在这段时间内经过的全部位置，还需要考虑轨迹在区间内部的变化。

6. 同时检查网格和边界条件

Nx = 100 在当前代码中表示空间点数，空间步长是 dx。原来的 np.arange(0, Lx, dx) 不包含右端点 Lx。如果模型要求 100 个节点同时包含杆的两个端点，可以改为：

Nx = 100
Lx = 0.15
x_points = np.linspace(0.0, Lx, Nx)
dx = x_points[1] - x_points[0]

这会改变原来的网格定义，需要同步检查边界离散方式。如果保留原网格，空间插值范围也应以实际的 x_points[0] 和 x_points[-1] 为准。

原代码中的 Tp[0] = Tp[1] 和 Tp[-1] = Tp[-2] 设置的是温度时间导数，不能仅凭这两行认定已经正确实现了绝热边界。边界节点的离散公式需要根据实际边界条件推导。

另外，range(1, 5) 只会循环四次。如果要生成五个快照，需要调整循环范围。绘图循环也建议使用 snapshot_index 等名称，避免顶层变量 k 覆盖 ODE 中作为导热系数使用的 k。

备注：内容仅供参考。