How to solve a cubic function with more than hundred pairs of known dependent and independent variable in python

Question

I try to find out how to solve a cubic function with independent variable (x) and dependent variable f(x) known but coefficient a, b, c and constant d unknown. I tried sympy but realized it only works for 4 pairs. Now I try to explore some possibilities to solve this (by finding the actual value for coefficient a/b/c and constant d). Any advice are greatly appreciated.

Below is the code which obviously does not work and returned an empty list because there are more than hundred of pairs.

from sympy import Eq, solve
from sympy.abc import a,b,c,d, x

formula = a*x**3 + b*x**2 + c*x + d  # general cubic formula

xs = [28.0, 29.0, 12.0, 12.0, 42.0, 35.0, 28.0, 30.0, 32.0, 46.0, 18.0, 28.0, 28.0, 64.0, 
38.0, 18.0, 49.0, 37.0, 25.0, 24.0, 42.0, 50.0, 12.0, 64.0, 23.0, 35.0, 22.0, 16.0, 44.0, 
77.0, 26.0, 44.0, 38.0, 37.0, 45.0, 42.0, 24.0, 42.0, 12.0, 46.0, 12.0, 26.0, 37.0, 15.0, 
67.0, 36.0, 43.0, 36.0, 45.0, 82.0,
44.0, 30.0, 33.0, 51.0, 50.0]

fxs = [59.5833333333333, 59.5833333333333, 10.0, 10.0, 47.0833333333333, 51.2499999999999, 
34.5833333333333, 88.75, 63.7499999999999, 34.5833333333333, 51.2499999999999, 10.0, 
63.7499999999999, 51.0, 59.5833333333333,47.0833333333333, 49.5625, 43.5624999999999, 
63.7499999999999, 10.0, 76.25, 47.0833333333333,10.0, 51.2499999999999,47.0833333333333,10.0, 
35.0, 51.2499999999999, 76.25, 100.0, 51.2499999999999, 59.5833333333333, 63.7499999999999, 
76.25, 100.0, 51.2499999999999, 10.0, 22.5, 10.0, 88.75, 10.0, 59.5833333333333, 
47.0833333333333, 34.5833333333333, 51.2499999999999, 63.7499999999999,63.7499999999999, 10.0, 
76.25, 62.1249999999999, 47.0833333333333, 10.0, 76.25, 47.0833333333333, 88.75]

sol = solve([Eq(formula.subs(x, xi), fx) for xi, fx in zip(xs, fxs)])
print(sol)  
[]

Answer 1

For curve-fitting, I recommend using scipy.optimize.curve_fit:

import numpy as np
import plotly.graph_objects as go  # only used to show output; not needed in answer
from scipy.optimize import curve_fit


def cubic(x, a, b, c, d):
    return a * x**3 + b * x**2 + c * x + d


(a, b, c, d), _ = curve_fit(cubic, xs, fxs)  # `xs` and `fxs` copied from the OP


x = np.linspace(min(xs), max(xs), 1000)


fig = go.Figure(go.Scatter(x=xs, y=fxs, mode='markers', name='data'))
fig.add_trace(go.Scatter(x=x, y=cubic(x, a, b, c, d), name='fit'))
fig.show()

Output: