Question: I am trying to estimate the Schott coefficients of a glass material given only its n_e(refraction index at e line) and V_e(Abbe number at e line). Schott is one way to represent the dispersion of a material, which is the different index of refraction (RI) at different wavelength.  In the figure above, the horizontal axis is the wavelength (in micrometer) and the vertical axis is index of refraction (This figure is based on the glass type named KZFH1). Because the glass dispersion have a common shape (higher at shorter wavelength and then tappers down), and the RI at key points (Fraunhofer lines) have a stable relationship, my thought is that I can use the definition of Abbe number and the general relation of different Fraunhofer line RI to create some data points, and use them to fit a curve: import numpy as np from scipy.optimize import curve_fit # Definition of the Schott function def _InvSchott(x, a0, a1, a2, a3, a4, a5): return np.sqrt(a0 + a1* x**2 + a2 * x**(-2) + a3 * x**(-4) + a4 * x**(-6) + a5 * x**(-8)) # Sample input, material parameter from a Leica Summilux patent n = 1.7899 V = 48 # 6 wavelengths, Fraunhofer symbols are not used due to there is another version that uses n_d and V_d shorter = 479.99 short = 486.13 neighbor = 546.07 middle = 587.56 longc = 643.85 longer = 656.27 # Refraction index of the corresponding wavelengths. # The linear functions are acquired from external regressions from 2000 glass materials n_long = 0.984 * n + 0.0246 # n_C' n_shorter = ( (n-1) / V) + n_long # n_F', from the definition of Abbe number n_short = 1.02 * n -0.0272 # n_F n_neighbor = n # n_e n_mid = 1.013 * n - 0.0264 # n_d n_longer = 0.982 * n + 0.0268 # n_C # The /1000 is to convert the wavelength from nanometer to micrometers x_data = np.array([longer, longc, middle, neighbor, short, shorter]) / 1000.0 y_data = np.array([n_longer, n_long, n_mid, n_neighbor, n_short, n_shorter]) # Provided estimate are average value from the 2000 Schott glasses popt, pcov = curve_fit(_InvSchott, x_data, y_data, [2.75118, -0.01055, 0.02357, 0.00084, -0.00003, 0.00001])  The x_data and y_data in this case are as follow: [0.65627 0.64385 0.58756 0.54607 0.48613 0.47999] [1.7844818 1.7858616 1.7867687 1.7899 1.798498 1.80231785]  And then I got the warning OptimizeWarning: Covariance of the parameters could not be estimated. The fit result were all but [inf inf inf inf inf inf]. I know this question has been asked a lot but I have not found a solution that works in this case yet. 6 data point is certainly a bit poor but this does satisfy the minimum, and Schott function is continuous, so I cannot figure out which part went wrong. TLDR: How do I find the coefficients for the function def _InvSchott(x, a0, a1, a2, a3, a4, a5): return np.sqrt(a0 + a1* x**2 + a2 * x**(-2) + a3 * x**(-4) + a4 * x**(-6) + a5 * x**(-8))  that fits the data below: x: [0.65627 0.64385 0.58756 0.54607 0.48613 0.47999] y: [1.7844818 1.7858616 1.7867687 1.7899 1.798498 1.80231785] 

Answer: Don't use sqrt during fitting, and don't fit this as a nonlinear model. Fit it as a linear model: import numpy as np from matplotlib import pyplot as plt schott_powers = np.arange(2, -9, -2) def inv_schott(lambd: np.ndarray, a: np.ndarray) -> np.ndarray: return np.sqrt(inv_schott_squared(lambd, a)) def inv_schott_squared(lambd: np.ndarray, a: np.ndarray) -> np.ndarray: terms = np.power.outer(lambd, schott_powers) return terms @ a def demo() -> None: # Sample input, material parameter from a Leica Summilux patent n = 1.7899 V = 48 # 6 wavelengths, Fraunhofer symbols are not used due to there is another version that uses n_d and V_d shorter = 479.99 short = 486.13 neighbor = 546.07 middle = 587.56 longc = 643.85 longer = 656.27 # Refraction index of the corresponding wavelengths. # The linear functions are acquired from external regressions from 2000 glass materials n_long = 0.984*n + 0.0246 # n_C' n_shorter = (n - 1)/V + n_long # n_F', from the definition of Abbe number n_short = 1.02*n - 0.0272 # n_F n_neighbor = n # n_e n_mid = 1.013*n - 0.0264 # n_d n_longer = 0.982*n + 0.0268 # n_C lambda_nm = np.array((longer, longc, middle, neighbor, short, shorter)) lambda_um = lambda_nm*1e-3 n_all = np.array((n_longer, n_long, n_mid, n_neighbor, n_short, n_shorter)) a, residuals, rank, singular = np.linalg.lstsq( a=np.power.outer(lambda_um, schott_powers), b=n_all**2, rcond=None, ) fig, ax = plt.subplots() lambda_hires = np.linspace(start=lambda_um.min(), stop=lambda_um.max(), num=501) ax.scatter(lambda_um, n_all, label='experiment') ax.plot(lambda_hires, inv_schott(lambda_hires, a), label='fit') ax.legend() plt.show() if __name__ == '__main__': demo()   You can observe the effect of truncating the polynomial order, though whether that's advisable is very difficult to say unless you produce more data: import numpy as np from matplotlib import pyplot as plt def inv_schott(lambd: np.ndarray, a: np.ndarray, powers: np.ndarray) -> np.ndarray: return np.sqrt(inv_schott_squared(lambd, a, powers)) def inv_schott_squared(lambd: np.ndarray, a: np.ndarray, powers: np.ndarray) -> np.ndarray: terms = np.power.outer(lambd, powers) return terms @ a def demo() -> None: # Sample input, material parameter from a Leica Summilux patent n = 1.7899 V = 48 # 6 wavelengths, Fraunhofer symbols are not used due to there is another version that uses n_d and V_d shorter = 479.99 short = 486.13 neighbor = 546.07 middle = 587.56 longc = 643.85 longer = 656.27 # Refraction index of the corresponding wavelengths. # The linear functions are acquired from external regressions from 2000 glass materials n_long = 0.984*n + 0.0246 # n_C' n_shorter = (n - 1)/V + n_long # n_F', from the definition of Abbe number n_short = 1.02*n - 0.0272 # n_F n_neighbor = n # n_e n_mid = 1.013*n - 0.0264 # n_d n_longer = 0.982*n + 0.0268 # n_C lambda_nm = np.array((longer, longc, middle, neighbor, short, shorter)) lambda_um = lambda_nm*1e-3 n_all = np.array((n_longer, n_long, n_mid, n_neighbor, n_short, n_shorter)) fig, ax = plt.subplots() lambda_hires = np.linspace(start=lambda_um.min(), stop=lambda_um.max(), num=501) ax.scatter(lambda_um, n_all, label='experiment') for lowest_power in range(0, -9, -2): powers = np.arange(2, lowest_power - 1, -2) a, residuals, rank, singular = np.linalg.lstsq( a=np.power.outer(lambda_um, powers), b=n_all**2, rcond=None, ) ax.plot(lambda_hires, inv_schott(lambda_hires, a, powers), label=f'powers to {lowest_power}') ax.legend() plt.show() if __name__ == '__main__': demo()