具有多个多项式拟合的曲面图答案

【问题标题】：Surface plot with multiple polynomial fits具有多个多项式拟合的曲面图
【发布时间】：2016-10-12 12:58:48
【问题描述】：

我所问的可能是不可能的，但我希望你们能帮助我。

所以我有两个二维数组，f1(x1) = y1 和 f2(x2) = y2。我想绘制这些比率的曲面图，因此 z 维度为 (f2(x2)/f1(x1))。不幸的是，无论我从哪个方向解决问题，我都会碰壁。

我的主要问题是每个数组的范围不同，x1 从 300 到 50,000，x2 从 300 到 1200。现在我很高兴假设所有的 f2(x2) = f2(1200) x2 > 1200。但是这个界限意味着我不可能以任何现实的方式拟合这个数据的多项式（我的第一组数据很好地由五阶多项式再现，而我的第二组数据最适合一阶多项式）。是否有另一种方法可以将函数拟合到 (x2,y2) 以便它获取边界外所有点的边界值？

我极其错误的尝试是这样的，

# x1 is an array from 300 to 50,000 in steps of 50
# x2 is an array from 300 to 1,150 in steps of 50

f1_fit = np.poly1d(np.polyfit(x1, y1, 5))
f2_fit = np.poly1d(np.polyfit(x2, y2, 1))

X, Y = np.meshgrid(x1, x2)
Z = (f2_fit(x2) / f1_fit(x1))

有趣的是，看似无害的问题实际上是一种痛苦。 :D

编辑：这是玩具数据的数量，

x1 = x2 = [  300.   350.   400.   449.   499.   548.   598.   648.   698.   748.
             798.   848.   897.   947.   997.  1047.  1097.  1147.  1196.  1246.
             1296.  1346.  1396.  1446.  1496.  1546.  1595.  1645.  1695.  1745.]

y1 = [  351.   413.   476.   561.   620.   678.   734.   789.   841.   891.
        938.   982.  1023.  1062.  1099.  1133.  1165.  1195.  1223.  1250.
        1274.  1298.  1320.  1340.  1360.  1378.  1395.  1411.  1426.  1441.]

y2 = [ 80.  75.  70.  65.  62.  58.  58.  52.  48.  46.  44.  41.  38.  35.  32.
       32.  29.  30.  30.  30.  30.  30.  30.  30.  30.  30.  30.  30.  30.  30.]

【问题讨论】：

我看到了解决这个问题的两种方法。您可以进行手动插值以获取“中间”数据点的值，也可以使用 scipy 的 interp2d 函数插值 docs.scipy.org/doc/scipy/reference/generated/…
如果你给我一些输入数据来玩，那么我会告诉你怎么做。
我已经添加了数据。
呃，您的玩具数据与您的描述不符 - x1 和 x2 相同，而您在问题中对它们的描述不同。
它们是相同的，因为我假设所有 x2 > 1150 的边界值为 f2(1150)。

标签： python matplotlib curve-fitting surface

【解决方案1】：

所以我设法解决了我的问题。如上所述，我通过设置 x1 = x2 对我的数据进行了一些预处理，从而推断出 f(x2) 的边缘值。

import numpy as np
import scipy.interpolate as interp
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

data1 = np.loadtxt('data1.dat')
data2 = np.loadtxt('data2.dat')

X = []
Y = []
Z = []

for i in data1:
    for j in data2:
        X.append(i[0])
        Y.append(j[0])
        Z.append(i[1]/j[1])

x_mesh,y_mesh, = np.meshgrid(np.linspace(300,50000,200), np.linspace(300,50000,200))

z_mesh = interp.griddata((X,Y),Z,(x_mesh,y_mesh),method='linear')

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
surf = ax.plot_surface(x_mesh,y_mesh,z_mesh,cstride=10,rstride=10,cmap='OrRd_r')
plt.show()

【讨论】：