I am using SciPy’s optimize.curve_fit to get polynomial coefficients for an implicit function, and that was successful. I now have y(x). I need x(y). If I instead do a curve_fit for y, x (as opposed to x, y), the fit is not as good near the origin, and it needs to be as good as possible. I am in no danger of over-fitting.

I figured okay, the fit is bad, what if I instead get an inverse polynomial using the fit polynomial. I do not know how to accomplish this. There are some other posts on SO that deal with this, but the examples are kind of hard for me to work through as they are dealing with more complicated polynomials than I am. In the range that I am interested in, y(x) and x(y) are both one-to-one functions (for every x there is a unique y and vice versa) so there will be no weirdness with multiple outputs (x) for a given input (y).

Honestly, using curve_fit for x(y) should have worked, but it doesn’t and I don’t understand why. Here is my code:

import matplotlib.pyplot as plt
import numpy as np
import scipy as sp

# Get x, y from implicit equation

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 # LHS
G = Y*d / (a*(1-(1-d)**(X+1))) # RHS

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

for item in cs.collections: # I understand this will be deprecated soon
    for item in item.get_paths()
        v = i.vertices
        x = v[:, 0]
        y = v[:, 1]
        x_arr = np.append(x_arr, x)
        y_arr = np.append(y_arr, y)

c_plot = plt.plot(x_arr, y_arr, label="c_plot")

# Fit resulting curve 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, yx_pcov = sp.optimize.curve_fit(func_yx, x_arr, y_arr)
print(f'yx_popt: {yx_popt}')
yx_plot = plt.plot(x_arr, func_yx(x, *yx_popt), label="yx_plot")

# Fit resulting curve 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

xy_popt, xy_pcov = sp.optimize.curve_fit(func_xy, y_arr, x_arr)
print(f'xy_popt: {xy_popt}')
xy_plot = plt.plot(func_xy(y, *xy_popt), y_arr, label="xy_plot") # Note: AXES ARE FLIPPED HERE to compare to yx_plot

# Show everything

ax = plt.gca()
ax.set_aspect('equal', adjustable='box')
plt.legend()
plt.show()