Python中逆多项式拟合的系数求解问题

一、直接拟合x(y)原点效果差的原因 #

1. 误差权重失衡：curve_fit默认采用均方误差，所有数据点权重均等。若数据集中大y（对应大x）的点数量更多或波动更大，拟合会优先降低这些区域的误差，牺牲原点附近（小y、小x）的精度。
2. 多项式基函数局限性：你使用的8次多项式a + b*y + c*y² + ...在原点的拟合能力依赖低次项系数，但全局误差平衡会导致低次项偏离最优值。
3. 数据分布不均：从等高线提取的点在原点附近密度较低，进一步拉低了该区域的拟合优先级。

二、两种解决方案 #

方案1：给x(y)拟合添加区域权重 #

通过curve_fit的sigma参数给原点附近的点分配更高权重，强制拟合优先保证该区域精度：

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize

# （省略原代码中提取x_arr、y_arr的部分，保持不变）

# 定义x(y)拟合函数
def func_xy(y, a, b, c, d, e, f, g, h, i):
    return a + b*y + c*y**2 + d*y**3 + e*y**4 + f*y**5 + g*y**6 + h*y**7 + i*y**8

# 设置权重：小y对应小sigma，权重更高（加epsilon避免除以0）
epsilon = 1e-6
sigma = y_arr + epsilon
xy_popt_weighted, _ = optimize.curve_fit(func_xy, y_arr, x_arr, sigma=sigma)

# 绘制带权重的拟合结果
plt.plot(func_xy(y_arr, *xy_popt_weighted), y_arr, label="带权重x(y)拟合")

方案2：从已拟合的y(x)逆推x(y) #

利用y(x)的单调性（一一映射），对目标y值求解方程func_yx(x) - y = 0得到对应x，再用这些高精度点拟合x(y)多项式：

# （省略原代码中拟合y(x)的部分，yx_popt已获取）

# 1. 加密原点附近的y采样点
y_sample = np.concatenate([
    np.linspace(0, 10, 200),  # 原点区域密集采样
    np.linspace(10, 800, 100)
])

# 2. 对每个y值求解对应的x
x_sample = []
for y_val in y_sample:
    def target_func(x):
        return func_yx(x, *yx_popt) - y_val
    # 牛顿迭代法求解，初始值设为y_val的粗略估计
    x_sol = optimize.newton(target_func, x0=y_val)
    x_sample.append(x_sol)
x_sample = np.array(x_sample)

# 3. 用高精度采样点拟合x(y)
xy_popt_from_yx, _ = optimize.curve_fit(func_xy, y_sample, x_sample)

# 绘制逆推拟合结果
plt.plot(func_xy(y_sample, *xy_popt_from_yx), y_sample, label="从y(x)逆推x(y)拟合")

三、完整整合代码 #

import matplotlib.pyplot as plt
import numpy as np
from scipy import optimize

# 从隐函数提取x、y数据
delta = 0.5
xrange = np.arange(0.0, 800.0, delta)
yrange = np.arange(0.0, 800.0, delta)
X, Y = np.meshgrid(xrange, yrange)

a = 0.04
d = 0.02
F = X
G = Y*d / (a*(1-(1-d)**(X+1)))

x_arr = np.zeros(0)
y_arr = np.zeros(0)
cs = plt.contour(X, Y, (F-G), [0])

for item in cs.collections:
    for i in item.get_paths():
        v = i.vertices
        x_arr = np.append(x_arr, v[:, 0])
        y_arr = np.append(y_arr, v[:, 1])

plt.plot(x_arr, y_arr, label="原始曲线")

# 拟合y(x)
def func_yx(x, a, b, c, d, e, f, g, h, i):
    return a + b*x + c*x**2 + d*x**3 + e*x**4 + f*x**5 + g*x**6 + h*x**7 + i*x**8

yx_popt, _ = optimize.curve_fit(func_yx, x_arr, y_arr)
plt.plot(x_arr, func_yx(x_arr, *yx_popt), label="y(x)拟合曲线")

# 方案1：带权重拟合x(y)
def func_xy(y, a, b, c, d, e, f, g, h, i):
    return a + b*y + c*y**2 + d*y**3 + e*y**4 + f*y**5 + g*y**6 + h*y**7 + i*y**8

epsilon = 1e-6
sigma = y_arr + epsilon
xy_popt_weighted, _ = optimize.curve_fit(func_xy, y_arr, x_arr, sigma=sigma)
plt.plot(func_xy(y_arr, *xy_popt_weighted), y_arr, label="带权重x(y)拟合")

# 方案2：从y(x)逆推x(y)
y_sample = np.concatenate([np.linspace(0, 10, 200), np.linspace(10, 800, 100)])
x_sample = []
for y_val in y_sample:
    def target_func(x):
        return func_yx(x, *yx_popt) - y_val
    x_sol = optimize.newton(target_func, x0=y_val)
    x_sample.append(x_sol)
x_sample = np.array(x_sample)

xy_popt_from_yx, _ = optimize.curve_fit(func_xy, y_sample, x_sample)
plt.plot(func_xy(y_sample, *xy_popt_from_yx), y_sample, label="从y(x)逆推x(y)拟合")

# 显示设置
ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.legend()
plt.show()

内容的提问来源于stack exchange，提问作者kriegersan