在 Python 中使用 Vandermonde 矩阵进行多元 (3D) 插值

Question

我似乎无法弄清楚如何将 Vandermonde 矩阵实现为多元插值。我能够得到实际的矩阵，但我不明白如何获取值（数组） c00,c01,c02... 。我知道 c = V/z，但我觉得我错过了一些东西（也许，不是除法？）。我也知道我需要以某种方式建立一个方程组（V 的列是每个 cij）。

你如何在 python 中做到这一点？

这是我目前所拥有的：

import numpy as np
x = [1, 1, 1, 2, 2, 2, 3, 3, 3]
y = [1, 2, 3, 1, 2, 3, 1, 2, 3]
z = [3.2, 4.4, 6.5, 2.5, 4.7, 5.8, 5.1, 3.6, 2.9]
numpy.polynomial.polynomial.polyvander2d(x, y, [2,2])
>>>array([[  1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.],
   [  1.,   2.,   4.,   1.,   2.,   4.,   1.,   2.,   4.],
   [  1.,   3.,   9.,   1.,   3.,   9.,   1.,   3.,   9.],
   [  1.,   1.,   1.,   2.,   2.,   2.,   4.,   4.,   4.],
   [  1.,   2.,   4.,   2.,   4.,   8.,   4.,   8.,  16.],
   [  1.,   3.,   9.,   2.,   6.,  18.,   4.,  12.,  36.],
   [  1.,   1.,   1.,   3.,   3.,   3.,   9.,   9.,   9.],
   [  1.,   2.,   4.,   3.,   6.,  12.,   9.,  18.,  36.],
   [  1.,   3.,   9.,   3.,   9.,  27.,   9.,  27.,  81.]])

np.dot(V, c.flat)和numpy.polynomial.polynomial.polyval2d(x, y, c)我认为必须以某种方式融入其中，但我不知道该怎么做。请帮忙！ 我应该输出： c = V \ z c = 0.97500 -5.27500 5.95000 -3.92500 19.82500 -21.55000 3.40000 -14.70000 18.50000

这是我得到这个例子的网站（他们使用了 MatLab）： https://ece.uwaterloo.ca/~dwharder/NumericalAnalysis/05Interpolation/multi/ 希望这会有所帮助！

Answer 1

给定位置 x 和 y：

x = [1, 1, 1, 2, 2, 2, 3, 3, 3]
y = [1, 2, 3, 1, 2, 3, 1, 2, 3]

和函数值z：

z = [3.2, 4.4, 6.5, 2.5, 4.7, 5.8, 5.1, 3.6, 2.9]

V 将是 Vandermonde 矩阵：

V = numpy.polynomial.polynomial.polyvander2d(x, y, [2,2])

V = array([[  1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.,   1.],
           [  1.,   2.,   4.,   1.,   2.,   4.,   1.,   2.,   4.],
           [  1.,   3.,   9.,   1.,   3.,   9.,   1.,   3.,   9.],
           [  1.,   1.,   1.,   2.,   2.,   2.,   4.,   4.,   4.],
           [  1.,   2.,   4.,   2.,   4.,   8.,   4.,   8.,  16.],
           [  1.,   3.,   9.,   2.,   6.,  18.,   4.,  12.,  36.],
           [  1.,   1.,   1.,   3.,   3.,   3.,   9.,   9.,   9.],
           [  1.,   2.,   4.,   3.,   6.,  12.,   9.,  18.,  36.],
           [  1.,   3.,   9.,   3.,   9.,  27.,   9.,  27.,  81.]])

每一行：

a = x[i]
b = y[i]
V[i,:] = [ 1, b, b², a, a*b, a*b², a², a²b, a²b²]

线性插值旨在解决：

$$z = V \cdot c$$

因此我们需要解决：

$$c = V^{-1} z$$

c = np.linalg.solve(V, z)

c 现在以与 Vandermonde 矩阵相同的顺序保存系数。

您可以手动评估它：

def poly_eval(x, y, z):
    return z[0] + z[1]*y + z[2]*np.power(y,2) + z[3]*x + z[4]*x*y + z[5]*x*np.power(y,2) + z[6]*np.power(x,2) + z[7]*np.power(x,2)*y + z[8]*np.power(x,2)*np.power(y,2)

或使用

np.polynomial.polynomial.polyval2d([1,2], [3,4], c)
Out[22]: array([-65.5, -88.4])