在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

适合2维矢量场的多项式将是两个二元多项式 - 一个用于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) 

而且，我想我发现了一个bug在​​- 这个版本提供了更精确的结果：

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 

它的工作原理是