计算插值后 3D 规则网格数据的导数答案

【问题标题】：Compute the derivative of a 3D regular grid data after interpolation计算插值后 3D 规则网格数据的导数
【发布时间】：2019-06-27 18:59:46
【问题描述】：

我的数据是一个常规的 3D 网格，作为一个最简单的例子，我在这里给出了一个简单立方体的简单 h5 文件的代码，我使用RegularGridInterpolator 函数进行插值。但是，我想知道如何更改我的代码，以便我可以从这个插值函数中获得导数？为了您的友好信息，我在这里给出了我的代码：

生成h5文件的代码：

import numpy as np
import h5py

def f(x,y,z):
   return 2 * x**3 + 3 * y**2 - z

x = np.linspace(-1, 1, 2)
y = np.linspace(-1, 1, 2)
z = np.linspace(-1, 1, 2)
mesh_data = f(*np.meshgrid(x, y, z, indexing='ij', sparse=True))

h5file = h5py.File('cube.h5', 'w')
h5file.create_dataset('/x', data=x)
h5file.create_dataset('/y', data=y)
h5file.create_dataset('/z', data=z)
h5file.create_dataset('/mesh_data', data=mesh_data)

h5file.close()

插值代码：

import numpy as np   
import h5py
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from scipy.interpolate import RegularGridInterpolator
from mpl_toolkits.mplot3d.art3d import Poly3DCollection, Line3DCollection
f = h5py.File('cube.h5', 'r')  #importing h5 file from gevolution output
list(f.keys())
dset = f[u'mesh_data']
dset.shape
dset.value.shape
dset[0:2,0:2,0:2]
x = np.linspace(-1, 1, 2)
y = np.linspace(-1, 1, 2)
z = np.linspace(-1, 1, 2)
def deriv_interp(x, y, z, dset):
    d_dset = np.gradient(dset)
    deriv_interpolators = [RegularGridInterpolator((x, y, z), d, method = 'linear') for d in d_dset]
    def interp(x,y,z):
        return np.dstack([d(x,y,z) for d in deriv_interpolators])
    return interp
pts = np.array([0.2, 0.9, 0.6]) 
deriv_interpolators(pts)

【问题讨论】：

标签： python numpy scipy interpolation h5py

【解决方案1】：

我认为RegularGridInterpolator 没有任何有用的属性来查找导数，因为它只能进行线性或最近邻插值，因此它并没有真正在引擎盖下创建任何类似导数的东西。你需要某种 n-d 多项式插值器，据我所知，它似乎没有在 scipy 中实现。当然，您总是可以在网格数据上使用np.gradient，然后为每个维度制作另一组迭代器。

def deriv_interp(x, y, z, dset):
    d_dset = np.gradient(dset)
    deriv_interpolators = [RegularGridInterpolator((x, y, z), d, method = 'linear') for d in d_dset]
    def interp(coord):
        return np.dstack([d(coord) for d in deriv_interpolators])
    return interp

还没有真正检查过，但这应该是一般的想法。

【讨论】：

根据您的建议，我已更改代码并编辑了我的问题（请参阅题为“插值代码”的问题部分中的第二个代码）。你能帮我找出我在这里的错误吗？ @丹尼尔

【解决方案2】：

不是 100% 确定您要达到的目标，也许您已经尝试过类似的方法，但是呢：

import jax
finterp=RegularGridInterpolator((x, y, z), meshdata, method = 'linear')
gradf=jax.grad(finterp)

或者你可以试试：

gradf=jax.grad(f)

【讨论】：