你可以试试 GNU 科学图书馆：https ://www.gnu.org/software/gsl/doc/html/interp.html#id5

我终于在 f90 中找到了几个函数，可以很好地为任何大小的网格执行双线性插值。这些例程可以在这里找到。

如果您有 1D 插值例程并且查询点位于常规网格上，则这是一个可行的想法。

您可以计算数据矩阵的 SVD，$Y = U\Sigma V^T$，插值的列$U$和$V$, 获得,$\hat{U}$和$\hat{V}$，然后将插值计算为$\hat{U} \Sigma \hat{V}^T$. 如果$Y$近似低秩，您可以丢弃对应于小的奇异值的奇异向量。

请注意，在插值的列时$U$和$V$，确定哪对节点最接近查询点只需要对所有列进行一次。

如果查询点的数量很大，通过这种矩阵-矩阵乘积形成输出可以比直接双线性插值更快。

下面是一些演示这个想法的 matlab 代码：

x = linspace(-1, 1, 100);
[X,Y] = meshgrid(x);
A = exp(-(X.^2+Y.^2)).*sin(2*pi*X.*Y);
xq = linspace(-1, 1, 1000);
[Xq,Yq] = meshgrid(xq);
tic
Aq = interp2(X, Y, A, Xq, Yq);
toc

tic
[U,S,V] = svd(A);
sv = diag(S);
Uq = interp1(x, U, xq);
Vq = interp1(x, V, xq);
AqSVD = Uq*(sv.*Vq');
toc
norm(AqSVD - Aq, inf)/norm(Aq, inf)

Elapsed time is 0.033942 seconds.
Elapsed time is 0.012933 seconds.

ans =

      1.29298230462322e-15