使用Python实现3D矢量场插值时出现报错的问题咨询

• scipy.interpolate.RegularGridInterpolator要求传入的网格轴参数为一维且严格递增的坐标序列，你代码中传入的old_points是np.meshgrid生成的三维网格数组，不符合参数要求，因此触发维度0数值非严格递增的报错。
• np.meshgrid默认使用indexing='xy'规则，生成的x,y,z数组形状为(Ny, Nx, Nz)，你强行将u,v,wreshape为(Nx, Ny, Nz)会打乱坐标与矢量分量的对应关系，即使解决了轴的问题也会得到错误插值结果。

以下是可直接运行的完整代码，核心修改点为单独提取一维网格轴传入插值器，统一网格索引规则避免坐标错位：

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import RegularGridInterpolator

# ---------------------- 1. 生成示例原始矢量场 ----------------------
Nx = 8
Ny = 8
Nz = 3
# 先定义各轴的一维严格递增坐标序列，后续插值要用到
x_axis = np.linspace(-1, 1, Nx)
y_axis = np.linspace(-1, 1, Ny)
z_axis = np.linspace(-1, 1, Nz)
# 用indexing='ij'保证数组形状为(Nx, Ny, Nz)，和轴顺序对应
x, y, z = np.meshgrid(x_axis, y_axis, z_axis, indexing='ij')

# 构造u,v,w分量
u = np.sin(np.pi * x) * np.cos(np.pi * y) * np.cos(np.pi * z)
v = -np.cos(np.pi * x) * np.sin(np.pi * y) * np.cos(np.pi * z)
w = (np.sqrt(2.0 / 3.0) * np.cos(np.pi * x) * np.cos(np.pi * y) *
     np.sin(np.pi * z))

# 可视化原始场
fig = plt.figure(dpi=150)
ax = fig.gca(projection='3d')
ax.quiver(x, y, z, u, v, w, length=0.2)
plt.title("原始矢量场")
plt.show()

# ---------------------- 2. 定义插值函数 ----------------------
def interpolate_field(old_axes, u, v, w, new_points):
    # old_axes: 原始网格各轴的一维坐标，格式为(x_axis, y_axis, z_axis)
    # u,v,w: 原始网格对应的矢量分量，形状为(Nx, Ny, Nz)
    # new_points: 待插值点坐标，形状为(N, 3)，每行是(x,y,z)
    u_interp = RegularGridInterpolator(old_axes, u)(new_points)
    v_interp = RegularGridInterpolator(old_axes, v)(new_points)
    w_interp = RegularGridInterpolator(old_axes, w)(new_points)
    return u_interp, v_interp, w_interp

# ---------------------- 3. 执行新网格插值 ----------------------
# 定义新网格参数
NNx = 20
NNy = 20
NNz = 3
new_x_axis = np.linspace(-1, 1, NNx)
new_y_axis = np.linspace(-1, 1, NNy)
new_z_axis = np.linspace(-1, 1, NNz)
new_x, new_y, new_z = np.meshgrid(new_x_axis, new_y_axis, new_z_axis, indexing='ij')

# 整理待插值点格式
new_points = np.vstack([new_x.ravel(), new_y.ravel(), new_z.ravel()]).T

# 调用插值函数，原始轴直接传一维的x_axis/y_axis/z_axis即可
u_int, v_int, w_int = interpolate_field(
    old_axes=(x_axis, y_axis, z_axis),
    u=u, v=v, w=w,
    new_points=new_points
)

# 插值结果恢复为网格形状，方便后续可视化或计算
u_int_grid = u_int.reshape(NNx, NNy, NNz)
v_int_grid = v_int.reshape(NNx, NNy, NNz)
w_int_grid = w_int.reshape(NNx, NNy, NNz)

# 可视化插值后的场
fig = plt.figure(dpi=150)
ax = fig.gca(projection='3d')
ax.quiver(new_x, new_y, new_z, u_int_grid, v_int_grid, w_int_grid, length=0.2)
plt.title("插值后矢量场")
plt.show()

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