NumPy Polyfit: Fix Cubic Polynomial Coefficients Python

Sampling four points can give a different cubic because of Vandermonde ordering/conditioning and small numerical/modeling errors in your F(Dr) evaluations — tiny errors get amplified when you solve a 4×4 Vandermonde system. The robust fixes are: (A) be sure your coefficient ordering matches np.roots (highest-degree first), (B) use a stable fitter like np.polyfit or an orthogonal basis (Chebyshev) or (best if available) symbolically expand NT·Dr to get exact coefficients. Below I show why your particular code likely produced the wrong leading coefficient and give concise, practical Python/numpy snippets you can copy into Dynamo (CPython 3).

Contents

• Why sampling at Dr = 0…3 produced different cubic polynomial coefficients
• Root causes: Vandermonde ordering, conditioning, and model error
• Correct, robust ways to compute cubic polynomial coefficients (np.polyfit, np.vander, symbolic)
• Practical Python/numpy code snippets you can run
• Sources
• Conclusion

Why sampling at Dr = 0…3 produced different cubic polynomial coefficients

You solved A·coeffs = F_vals where A is a Vandermonde matrix built from Dr = [0,1,2,3]. In exact arithmetic, four exact samples of a true cubic uniquely recover the cubic. But in practice three things typically go wrong and produce a different leading coefficient (and thus a different root):

• Ordering mismatch: np.roots expects coefficients in descending-power order (a·x^3 + b·x^2 + c·x + d). If you accidentally produce coeffs in increasing order (d, c, b, a) and pass them straight to np.roots, the polynomial is wrong. The numpy polynomial APIs use different conventions, so watch out.
• Ill-conditioning of Vandermonde: Vandermonde matrices amplify tiny errors in the right-hand side into much larger errors in the solution for the leading coefficients; this is a classic numerical effect for polynomial interpolation (see Vandermonde references below).
• Model / evaluation error: if NT(Dr) isn’t an exact cubic (or your NT(Dr)Dr - 8M evaluation includes higher-order terms, approximations or numerical rounding), interpolating from only four samples recovers the unique cubic that fits those four sampled values — not necessarily your analytic cubic.

A quick sanity check: if you evaluate F using the analytic coefficients you quoted and then solve, you should recover the same coefficients exactly (within floating precision). If you don’t, ordering/ dtype / indexing is the likely culprit.

Root causes: Vandermonde ordering, conditioning, and model error

Short diagnostics to run before changing approach

• Check ordering: is your A built with columns [x3, x2, x, 1] (descending) or [1, x, x2, x3] (increasing)? np.vander(x, 4) by default uses descending order. np.polynomial.polynomial.polyfit returns coefficients in increasing order, so you must reverse before calling np.roots. See the numpy.polynomial docs for ordering details.
• Check dtype and scale: ensure x and F_vals are float64 (use dtype=float) so you don’t get integer truncation.
• Check matrix conditioning: np.linalg.cond(V) will tell you if the Vandermonde is poorly conditioned. If cond >> 1e12 you’ll likely see large coefficient errors.
• Verify the model: compute y_poly = np.polyval([5.694, -0.942, 2540.505, -73612.480], x) and compare to your evaluated F_vals with np.allclose(y_poly, F_vals, rtol=1e-10, atol=1e-8). If that fails, your NT(Dr) evaluations are not equal to the analytic cubic — fitting four points will therefore give a different cubic.

Relevant reading on Vandermonde ill-conditioning: Vandermonde matrix - Wikipedia and the empirical study How Bad Are Vandermonde Matrices?.

Correct, robust ways to compute cubic polynomial coefficients (np.polyfit, np.vander, symbolic)

Pick one of these depending on what you control:

1. Symbolic expansion (best / exact)
   If you have an analytic expression for NT(Dr) (not just a black-box evaluator), expand symbolically and collect coefficients. That yields exact coefficients and reproduces the hand-derived root. Use SymPy if available.

2. Numerically stable least-squares fit with many samples (practical, usually best when NT is numeric)
   Sample F(Dr) at many points across the domain that contains the root and use np.polyfit(x, y, 3) (least-squares + SVD). Because it’s overdetermined, tiny local evaluation noise has less influence than solving a 4×4 Vandermonde exactly.

3. Correct Vandermonde solve (works if F is exactly cubic and you manage numerical issues)
   Use V = np.vander(x, 4) (default descending), enforce float dtype, and solve coeffs = np.linalg.solve(V, y). If you must build A manually, make sure the column order matches np.roots.

4. Fit in an orthogonal basis (Chebyshev / Legendre) then convert to power basis
   This reduces ill-conditioning when nodes are poorly distributed. Example: Chebyshev.fit(...) then convert() to Polynomial and reverse coefficients for np.roots.

5. If you need extra numerical accuracy, use higher precision (mpmath / SymPy) or increase sampling + least-squares.

Which is best and why? If you can symbolically expand NT·Dr, do that — it’s exact. If not, np.polyfit on a dense sampling or a Chebyshev fit is the simplest robust numeric approach and will generally reproduce the analytic cubic (within floating precision) and the root ~17.360. Directly solving a small Vandermonde (4×4) is fragile unless you verify ordering, dtype and conditioning.

Practical Python/numpy code snippets you can run

Below are short, copy-pasteable snippets. Replace NT(…) and M with your implementation/values.

A — quick sanity test using the analytic polynomial (shows that a correct Vandermonde solve reproduces your coefficients)

import numpy as np

# analytic coefficients from your problem (descending order)
a, b, c, d = 5.694, -0.942, 2540.505, -73612.480

x = np.array([0., 1., 2., 3.], dtype=float)
y = a*x**3 + b*x**2 + c*x + d

V = np.vander(x, 4) # columns: x^3, x^2, x, 1 (descending powers)
coeffs = np.linalg.solve(V, y) # coeffs in descending order [a, b, c, d]
print('recovered coeffs:', coeffs)
print('roots:', np.roots(coeffs))
# you should see root ~ 17.360

B — recommended robust numeric method (np.polyfit over many samples)

import numpy as np

def NT(Dr):
 # your numeric NT(Dr) implementation here
 return ...

M = ... # your M value (numeric)
x = np.linspace(0.0, 40.0, 201) # sample across domain containing root
y = np.array([NT(xx)*xx - 8*M for xx in x], dtype=float)

coeffs = np.polyfit(x, y, 3) # returns [a,b,c,d] (descending)
print('coeffs (a,b,c,d):', coeffs)
roots = np.roots(coeffs)
print('all roots:', roots)
# pick the real root near expected value
real_roots = [r.real for r in roots if abs(r.imag) < 1e-8]
root_near = min(real_roots, key=lambda r: abs(r - 17.36))
print('root near 17.36:', root_near)
print('max residual (fit vs samples):', np.max(np.abs(y - np.polyval(coeffs, x))))

C — Chebyshev fit (more stable basis) and convert to power basis

import numpy as np
from numpy.polynomial import Chebyshev, Polynomial

def NT(Dr): ...
M = ...

x = np.linspace(0.0, 40.0, 201)
y = np.array([NT(xx)*xx - 8*M for xx in x], dtype=float)

cheb = Chebyshev.fit(x, y, 3) # Chebyshev fit (stable)
poly = cheb.convert(kind=Polynomial) # convert to power basis
# poly.coef is [c0, c1, c2, c3] (increasing); reverse for np.roots
coeffs_desc = poly.coef[::-1]
print('coeffs (descending):', coeffs_desc)
print('root near 17.36:', min([r.real for r in np.roots(coeffs_desc) if abs(r.imag) < 1e-8],
 key=lambda r: abs(r-17.36)))

D — symbolic exact extraction (if you can create NT symbolically)

from sympy import symbols, expand, Poly
Dr = symbols('Dr')

# NT_sym must be a SymPy expression built from your model
NT_sym = ... # e.g. some combination of symbols/operations
M_sym = symbols('M')
expr = expand(NT_sym*Dr - 8*M_sym)
P = Poly(expr, Dr)
coeffs_sym = [float(P.coeff_monomial(Dr**k)) for k in (3,2,1,0)] # descending
print('exact coefficients:', coeffs_sym)
# find numerical roots (either via SymPy or numpy)
import numpy as np
print('roots:', np.roots(coeffs_sym))

Practical tips while you try the snippets

• Always verify ordering before calling np.roots. If you used a function that returns increasing-order coefficients (like np.polynomial.polynomial.polyfit), reverse with [::-1].
• Check np.linalg.cond(V) if you solve a Vandermonde directly.
• Use float dtype everywhere: dtype=float or 1.0 literals.
• If np.polyfit raises a RankWarning, either increase sampling, reduce degree, or center/scale x (e.g. fit on (x - x.mean())/x.std() and convert back). See the NumPy docs for details.

Sources

• https://numpy.org/doc/stable/reference/generated/numpy.polyfit.html
• https://numpy.org/doc/stable/reference/generated/numpy.polynomial.polynomial.polyfit.html
• https://numpy.org/doc/stable/reference/generated/numpy.roots.html
• https://en.wikipedia.org/wiki/Vandermonde_matrix
• https://arxiv.org/abs/1504.02118
• https://people.sc.fsu.edu/~jburkardt/py_src/vandermonde_interp_1d/vandermonde_interp_1d.html
• https://stackoverflow.com/questions/18767523/fitting-data-with-numpy
• https://stackoverflow.com/questions/58651514/python-polynomial-roots-are-inaccurate
• https://github.com/numpy/numpy/issues/22104
• https://github.com/numpy/numpy/issues/9170

Conclusion

Your immediate fix is to verify coefficient ordering and dtype, then use a more stable method: either symbolically expand NT·Dr (exact) or fit numerically with np.polyfit (least-squares on many sample points) or fit in a Chebyshev basis and convert. Those approaches avoid Vandermonde fragility and will produce the cubic polynomial coefficients that reproduce the hand-derived root (Dr ≈ 17.360). If you want, paste your NT(Dr) implementation here and I’ll give the minimal code you can drop into Dynamo to produce the exact coefficients and the root.