在python中为向量场拟合多项式函数

Question

首先，感谢大家为 stackoverflow 所做的出色工作……你们太棒了，已经帮助我很多次了。关于我的问题：我有一系列格式为（VectorX、VectorY、StartingpointX、StartingpointY）的向量

data = [(-0.15304757819399128, -0.034405679205349315, -5.42877197265625, 53.412933349609375), (-0.30532995491023485, -0.21523935094046465, -63.36669921875, 91.832427978515625), (-0.15872430479453215, -0.077999419482978283, -67.805389404296875, 81.001983642578125), (-0.36415549211687903, -0.33757147194808113, -59.015228271484375, 82.976226806640625), (0.0, 0.0, 0.0, 0.0), (-0.052973530805275004, 0.098212384392411423, 19.02667236328125, -13.72125244140625), (-0.34318724086483599, 0.17123742336019632, 80.0394287109375, 108.58499145507812), (0.19410169197834648, -0.17635303976555861, -55.603790283203125, -76.298828125), (-0.38774018337716143, -0.0824692384322816, -44.59942626953125, 68.402496337890625), (0.062202543524108478, -0.37219011831012949, -79.828826904296875, -10.764404296875), (-0.56582988168383963, 0.14872365390732512, 39.67657470703125, 97.303192138671875), (0.12496832467900276, -0.12216653754859408, 24.65948486328125, -30.92584228515625)]

当我绘制矢量场时，它看起来像这样：

import numpy as np
import matplotlib.pyplot as plt

def main():
    # Format Data...
    numdata = len(data)
    x = np.zeros(numdata)
    y = np.zeros(numdata)
    u = np.zeros(numdata)
    v = np.zeros(numdata)
    for i,el in enumerate(data):
        x[i] = el[2]
        y[i] = el[3]
        # length of vector 
        z[i] = math.sqrt(el[0]**2+el[1]**2)
        u[i] = el[0]
        v[i] = el[1]

    # Plot
    plt.quiver(x,y,u,v )
    # showing the length with color 
    plt.scatter(x, y, c=z)
    plt.show()


main()

我想创建一个多项式函数来拟合整个区域的连续向量场。经过一些研究，我发现了以下用于在二维中拟合多项式的函数。问题是，它只接受一个值作为拟合值。

def polyfit2d(x, y, z, order=3):
    ncols = (order + 1)**2
    G = np.zeros((x.size, ncols))
    ij = itertools.product(range(order+1), range(order+1))
    for k, (i,j) in enumerate(ij):
        G[:,k] = x**i * y**j
    m, _, _, _ = np.linalg.lstsq(G, z)
    return m

def polyval2d(x, y, m):
    order = int(np.sqrt(len(m))) - 1
    ij = itertools.product(range(order+1), range(order+1))
    z = np.zeros_like(x)
    for a, (i,j) in zip(m, ij):
        z += a * x**i * y**j
    return z

此外，当我尝试拟合向量的一维长度时，从 polyval2d 返回的值完全关闭。有谁知道一种方法来获得一个拟合函数，该函数将为网格中的任何点返回一个向量 (x,y)？

谢谢！

Answer 1

拟合二维向量场的多项式将是两个二元多项式 - 一个用于 x 分量，一个用于 y 分量。换句话说，您的最终多项式拟合将类似于：

P(x,y)  = ( x + x*y, 1 + x + y ) 

所以您必须拨打polyfit2d 两次。这是一个例子：

import numpy as np
import itertools

def polyfit2d(x, y, z, order=3):
    ncols = (order + 1)**2
    G = np.zeros((x.size, ncols))
    ij = itertools.product(range(order+1), range(order+1))
    for k, (i,j) in enumerate(ij):
        G[:,k] = x**i * y**j
    m, _, _, _ = np.linalg.lstsq(G, z)
    return m

def fmt1(x,i):
  if i == 0:
    return ""
  elif i == 1:
    return x
  else:
    return x + '^' + str(i)

def fmt2(i,j):
  if i == 0:
    return fmt1('y',j)
  elif j == 0:
    return fmt1('x',i)
  else:
    return fmt1('x',i) + fmt1('y',j)

def fmtpoly2(m, order):
  for (i,j), c in zip(itertools.product(range(order+1), range(order+1)), m):
    yield ("%f %s" % (c, fmt2(i,j)))

xs = np.array([ 0, 1, 2, 3] )
ys = np.array([ 0, 1, 2, 3] )

zx = np.array([ 0, 2, 6, 12])
zy = np.array([ 1, 3, 5, 7])

mx = polyfit2d(xs, ys, zx, 2)
print "x-component(x,y) = ", ' + '.join(fmtpoly2(mx,2)) 

my = polyfit2d(xs, ys, zy, 2)
print "y-component(x,y) = ", ' + '.join(fmtpoly2(my,2))

在这个例子中，我们的向量场是：

at (0,0):    (0,1)
at (1,1):    (2,3)
at (2,2):    (6,5)
at (3,3):    (12,7)

另外，我想我在polyval2d 中发现了一个错误——这个版本给出了更准确的结果：

def polyval2d(x, y, m):
    order = int(np.sqrt(len(m))) - 1
    ij = itertools.product(range(order+1), range(order+1))
    z = np.zeros_like(x)
    for a, (i,j) in zip(m, ij):
        z = z + a * x**i * y**j
    return z