Curve.fit optimization error: 'Covariance of the parameters could not be estimated

Question

I'm trying to fit some data into a function, but scipy curve_fit gives me the following message: "OptimizeWarning: Covariance of the parameters could not be estimated" and won't plot or print the parameters of the fit. I believe the power inside the exp function might be the problem, but I don't know how to work around this error. I tried changing the initial parameters but it won't help. Here's the data and the function I want to fit.

Data below as comment format to paste in excel

Thanks in advance.

import pandas as pd #Module for reading data
import numpy as np #Math module
import matplotlib.pyplot as plt #Plot module
from scipy.optimize import curve_fit #Curve fit module 

var = pd.read_excel("Data.xlsx") # Read in same directory as file
print(var) #Print File

x1 = list(var['H(T)']) # x values 
y1 = list(var['J (A/cm^2)']) # y values

######### FUNCTION ########

def func2(b, Jc1, Jc2, Bl, Bm, y):
   return Jc1*np.exp(-b/Bl)+Jc2*(b/Bm)*np.exp((1/y)*(1-(b/Bm)**y))

#p0 = [23987.20423,12345.20423,.34,1,1.5] 

popt1, pcov1 = curve_fit(func2, x1, y1, #p0) # Call the module for fit 
print(popt1) # Print values in terminal 

plt.plot(x1, y1,'b*', label='Data 2-Step') # Plot data 
plt.plot(x1, func2(x1, *popt1) , 'k--') # plot fit

plt.xlabel('Field(T)') # name of x label 
plt.ylabel('Jc (A/cm^2)') # name of y label 
plt.title('Jc at 77K') # name of graph 

plt.legend() # legends 
plt.grid(True) # grid ? 

plt.tight_layout() # to look good 

plt.savefig('fitJc77K.png') # saves image 
plt.show() # plot 

#            H(T)    J (A/cm^2)
#    0  -0.000625  53831.204230
#    1   0.099721  48796.160850
#    2   0.200126  40844.694180
#    3   0.300065  32798.048680
#    4   0.400528  27049.519580
#    5   0.500758  23531.580950
#    6   0.600639  21316.537570
#    7   0.701246  19782.205290
#    8   0.801331  18735.695240
#    9   0.901387  17920.634920
#    10  1.001922  16973.087830
#    11  1.101876  16800.300530
#    12  1.202004  16390.637040
#    13  1.302481  15986.217990
#    14  1.402464  15600.469840
#    15  1.502709  15291.792590
#    16  1.603156  14915.839150
#    17  1.703227  14638.497350
#    18  1.803399  14316.516400
#    19  1.903890  13946.912170
#    20  2.003945  13532.901590
#    21  2.104074  13162.857140
#    22  2.204638  12655.183070
#    23  2.304548  12171.305820
#    24  2.404604  11015.720630
#    25  2.505168  11324.922750
#    26  2.605005  10736.685710
#    27  2.705279  10095.360850
#    28  2.805771   9444.359790
#    29  2.905536   8200.046140
#    30  3.005882   8263.225400
#    31  3.106156   7593.037460
#    32  3.206066   6908.264970
#    33  3.306484   6322.300950
#    34  3.406685   5758.770370
#    35  3.506523   4535.793440
#    36  3.607015   3777.405290
#    37  3.707288   4103.506880
#    38  3.807199   3611.456510
#    39  3.907544   3186.223490
#    40  4.007819   2081.366770
#    41  4.107584   1669.132320
#    42  4.208147   1527.562750
#    43  4.308421   1374.063490
#    44  4.408404   1500.655660
#    45  4.509114   1260.510480
#    46  4.609096   1081.811560
#    47  4.708935    726.452490
#    48  4.809935    846.458580
#    49  4.910063    410.075915
#    50  5.010046     93.699340
#    51  5.110609     65.085630
#    52  5.210593     33.309590
#    53  5.310866     10.023870
#    54  5.411431      8.439660
#    55  5.511414      0.161990
#    56  5.611542      1.470010
#    57  5.711960     -0.574940
#    58  5.812089      1.173210
#    59  5.912217     -1.641310

Answer 1

Part of your problem is that you weren't able to find good initial parameters. More of your problem is that you didn't provide sane bounds, and bounds are crucial to prevent divide-by-zero and negative powers. You must disallow negative magnetic field values for the same reason.

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


field, density = np.array((
    [-0.000625,   53831.204230],
    [ 0.099721,   48796.160850],
    [ 0.200126,   40844.694180],
    [ 0.300065,   32798.048680],
    [ 0.400528,   27049.519580],
    [ 0.500758,   23531.580950],
    [ 0.600639,   21316.537570],
    [ 0.701246,   19782.205290],
    [ 0.801331,   18735.695240],
    [ 0.901387,   17920.634920],
    [ 1.001922,   16973.087830],
    [ 1.101876,   16800.300530],
    [ 1.202004,   16390.637040],
    [ 1.302481,   15986.217990],
    [ 1.402464,   15600.469840],
    [ 1.502709,   15291.792590],
    [ 1.603156,   14915.839150],
    [ 1.703227,   14638.497350],
    [ 1.803399,   14316.516400],
    [ 1.903890,   13946.912170],
    [ 2.003945,   13532.901590],
    [ 2.104074,   13162.857140],
    [ 2.204638,   12655.183070],
    [ 2.304548,   12171.305820],
    [ 2.404604,   11015.720630],
    [ 2.505168,   11324.922750],
    [ 2.605005,   10736.685710],
    [ 2.705279,   10095.360850],
    [ 2.805771,    9444.359790],
    [ 2.905536,    8200.046140],
    [ 3.005882,    8263.225400],
    [ 3.106156,    7593.037460],
    [ 3.206066,    6908.264970],
    [ 3.306484,    6322.300950],
    [ 3.406685,    5758.770370],
    [ 3.506523,    4535.793440],
    [ 3.607015,    3777.405290],
    [ 3.707288,    4103.506880],
    [ 3.807199,    3611.456510],
    [ 3.907544,    3186.223490],
    [ 4.007819,    2081.366770],
    [ 4.107584,    1669.132320],
    [ 4.208147,    1527.562750],
    [ 4.308421,    1374.063490],
    [ 4.408404,    1500.655660],
    [ 4.509114,    1260.510480],
    [ 4.609096,    1081.811560],
    [ 4.708935,     726.452490],
    [ 4.809935,     846.458580],
    [ 4.910063,     410.075915],
    [ 5.010046,      93.699340],
    [ 5.110609,      65.085630],
    [ 5.210593,      33.309590],
    [ 5.310866,      10.023870],
    [ 5.411431,       8.439660],
    [ 5.511414,       0.161990],
    [ 5.611542,       1.470010],
    [ 5.711960,      -0.574940],
    [ 5.812089,       1.173210],
    [ 5.912217,      -1.641310],
)).T
field = np.clip(field, a_min=0, a_max=None)


def func2(b: np.ndarray,
          Jc1: float, Bl: float,
          Jc2: float, Bm: float, y: float) -> np.ndarray:
    return (
        Jc1*np.exp(-b/Bl) +
        Jc2*b/Bm*np.exp(
            (1 - (b/Bm)**y)/y
        )
    )


p0 = (5e4, 2, -1e4, 1, 2)
popt1, _ = curve_fit(
    f=func2,
    xdata=field, ydata=density,
    p0=p0,
    bounds=(
        (-np.inf,  0, -np.inf,  0,  0.1),
        ( np.inf, 10,  np.inf, 10, 10.0),
    ),
)
print(popt1)

fig, ax = plt.subplots()
ax.scatter(field, density, marker='+', label='Data 2-Step')
ax.plot(field, func2(field, *p0), label='p0')
ax.plot(field, func2(field, *popt1), label='k--')

ax.set_xlabel('Field (T)')
ax.set_ylabel('Jc (A/cm²)')
ax.set_title('Jc at 77°K')

ax.legend()
ax.grid(True)
plt.show()

[ 5.55993112e+04  1.44996434e+00 -1.54454798e+04  5.87121864e-01
  1.71151312e+00]