pythonscipycurve-fitting

SciPy curve_fit covariance of the parameters could not be estimated

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.

Glass dispersion plot

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]
Solution

  • 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()

    fit

    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()

    all orders