0

I am trying to curve fit the following equation with parameters d, D, Ar, Tr each of them bounded in some range. The physical constants are: gamma = 26.76E7, n = 6.59E28, Ad = 2.099E-20 The equation is broken into several parts.

#Fitting function
def func(x, d, D, Ar, Tr):
    pi = math.pi
    #define the variables
    w1 = 2 * pi * x
    z1 = (2 * w1 * d**2 / D)**0.5
    w2 = 2 *w1
    z2 = ( 2 * w2 * d**2 / D)**0.5
    J1 = (
        (1 + 5 * z1 / 8 + z1**2 / 8) / 
        ( 1 + z1 + z1**2 / 2 + z1**3 / 6 + z1**5 / 81 + z1**6 / 648)
    )
    J2 = (
        (1 + 5 * z2 / 8 + z2**2 / 8) /
        ( 1 + z2 + z2**2 / 2 + z2**3 / 6 + z2**5 / 81 + z2**6 / 648)
    )

    gamma = 26.76E7      #unit = 1/(T*s)
    n = 6.59E28
    
    #define the normalization constant A_d
    Ad = (8 / 45) * pi * gamma**4* sc.hbar**2 * (
        sc.mu_0 / (4 * pi)
    )**2 * n
    
    R1_diff = Ad * (J1 + 4 * J2 ) / (d * D )
    
    R1_rot = Ar * (Tr / (1 + w1**2 * Tr**2) + 4 * Tr / (
        1 + w2**2 * Tr**2)
    )
    
    R1_IL = R1_diff + R1_rot
    
    return R1_IL

This is the x and y data sets -

#Experimental x and y data points    
xData = np.array([
2.00E+07
1.42E+07
1.01E+07
7.16E+06
7.16E+06
5.09E+06
3.61E+06
2.57E+06
1.82E+06
1.29E+06
9.20E+05
6.53E+05
4.63E+05
3.29E+05
2.34E+05
1.66E+05
1.18E+05
8.39E+04
5.96E+04
4.24E+04
3.00E+04])

yData = np.array([
7.89E+00
8.79E+00
9.58E+00
1.02E+01
1.03E+01
1.08E+01
1.13E+01
1.17E+01
1.20E+01
1.23E+01
1.25E+01
1.28E+01
1.29E+01
1.30E+01
1.31E+01
1.31E+01
1.35E+01
1.31E+01
1.33E+01
1.35E+01
1.34E+01])

The range for the parameters are: d = [3.00E-10, 4.00E-10], D = [12.5E-12, 14.00E-12], Ar = [3.00E9, 4.00E9], Tr = [1.00E-10, 2.5E-10]

This is the code that I wrote which is not giving me the correct fit and the right parameter values.

popt, pcov = curve_fit(func, xData, yData, bounds = ([3.00E-10,12.5E-12,2.5E9,1.00E-10],[4.0E-10,14.0E-12,4.00E9,2.00E-10]))
print(popt)

#x values for the fitted function
xFit = np.arange(3.00E+04, 2.00E+07, 10.00)
 
#Plot the fitted function
plt.loglog(xFit, func(xFit, *popt), 'r') #label='fit params: d=%5.3f, D=%5.3f, Ar=%5.3f, Tr=%5.3f' % tuple(popt))
 
#Plot experimental data points
plt.loglog(xData, yData, 'bo') #label='experimental-data')
plt.xlabel('f1')
plt.ylabel('R1')
plt.legend()
plt.show()

Anything would be really appreciated. The curve_fit is way off when I run this code.

mikuszefski
  • 3,943
  • 1
  • 25
  • 38
  • The bound keyword in `curve_fit` takes `( list_of_all_lower_bounds, list_of_all_upper_bounds)` Check the example given in the docs. – mikuszefski Jun 17 '21 at 09:39
  • Hint: It would be better if you put the scaling of the orders of magnitude inside your function and let the fitted parameters be of order 1. Much more stable in the fit algorithm. – mikuszefski Jun 17 '21 at 09:41

1 Answers1

0

For information only :

The power function below appears to be well fitted to the data.

enter image description here

JJacquelin
  • 1,529
  • 1
  • 9
  • 11