1

I'm trying to fit Einstein approximation of resistivity in a solid in a set of experimental data. I have resistivity vs temperature (from 200 to 4 K)

import xlrd as xd
import matplotlib.pyplot as plt
import numpy as np
import pylab as pl
import scipy as sp
from scipy.optimize import curve_fit

#retrieve data from file
data = pl.loadtxt('salita.txt')
Temp = data[:, 1]
Res = data[:, 2]

#define fitting function
def einstein_func( T, ro0, AE, TE):
    nl = np.sinh(TE/(2*T))
    return ro0 + AE*nl*T

p0 = sp.array([1 , 1, 1])

coeffs, cov = curve_fit(einstein_func, Temp, Res, p0)

But I get these warnings

crio.py:14: RuntimeWarning: divide by zero encountered in divide
  nl = np.sinh(TE/(2*T))
crio.py:14: RuntimeWarning: overflow encountered in sinh
  nl = np.sinh(TE/(2*T))
crio.py:15: RuntimeWarning: divide by zero encountered in divide
  return ro0 + AE*np.sinh(TE/(2*T))*T
crio.py:15: RuntimeWarning: overflow encountered in sinh
  return ro0 + AE*np.sinh(TE/(2*T))*T
crio.py:15: RuntimeWarning: invalid value encountered in multiply
  return ro0 + AE*np.sinh(TE/(2*T))*T
Traceback (most recent call last):
  File "crio.py", line 19, in <module>
    coeffs, cov = curve_fit(einstein_func, Temp, Res, p0)
  File "/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/minpack.py", line 511, in curve_fit
    raise RuntimeError(msg)
RuntimeError: Optimal parameters not found: Number of calls to function has reached maxfev = 800.

I don't understand why it keeps saying that there is a divide by zero in sinh, since I have strictly positive values. Varying my starting guess has no effect on it.

EDIT: My dataset is organized like this:

4.39531E+0  1.16083E-7
4.39555E+0  -5.92258E-8
4.39554E+0  -3.79045E-8
4.39525E+0  -2.13213E-8
4.39619E+0  -4.02736E-8
4.43130E+0  -1.42142E-8
4.45900E+0  -2.60594E-8
4.46129E+0  -9.00232E-8
4.46181E+0  1.42142E-7
4.46195E+0  -2.13213E-8
4.46225E+0  4.26426E-8
4.46864E+0  -2.60594E-8
4.47628E+0  1.37404E-7
4.47747E+0  9.47612E-9
4.48008E+0  2.84284E-8
4.48795E+0  1.35035E-7
4.49804E+0  1.39773E-7
4.51151E+0  -1.75308E-7
4.54916E+0  -1.63463E-7
4.59176E+0  -2.36902E-9

where the first column is temperature and the second one is resistivity (negative values are due to noise in trial current since the sample is a PbIn alloy which becomes superconductive at temperature lower than 6.7-6.9K, here we are at 4.5K).

Argument I'm providing to sinh are Numpy arrays, with a linear function ro0 + AE*T my code works. I've tried with scipy.optimize.minimize but the result is the same. Now I see that I have almost nine hundred values in my file, may that be the problem?

I have edited my dataset removing some lines and now the only warning showing is 

RuntimeWarning: overflow encountered in sinh

How can I workaround it?

iacolippo
  • 4,133
  • 25
  • 37
  • It might help to post a sample of lines from `salita.txt` – xnx May 12 '15 at 20:32
  • It's not a divide by zero per say, but rather an overflow/underflow in the `sinh` function. Make sure the arguments you are providing to it are reasonable. If you are still having issues, try a different solver from `scipy.optimize.minimize` which is more flexible than `scipy.optimize.curve_fit`. – rth May 12 '15 at 21:07
  • I noticed now that minimize gives me a different error: `np.array has no attribute 'lower'`. I think that means I'm giving it an incorrect argument, but I can't figure how could I write it. I was able to do the fit in MATLAB, but now I'm just curious on how to do this in Python. – iacolippo May 14 '15 at 08:50

1 Answers1

1

Here are a couple of observations that could help:

  • You could try the least-squares fit directly with leastsq, providing the Jacobian, which might help tame it.

  • I'm guessing you don't want the superconducting temperatures in your data set at all if you're fitting to an Einstein model (do you have a source for this eqn, btw?)

  • Do make sure your initial guesses are as good as they could possibly be (ro0=AE=TE=1 probably won't cut it).

  • Plot your data and make sure there aren't any weird artefacts

  • You seem to be indexing your data array in the wrong way in your code example: if the data is structured as you say, you want:

    Temp = data[:, 0] Res = data[:, 1]

(Python indexes start at 0).

xnx
  • 24,509
  • 11
  • 70
  • 109
  • I have other columns in my file but I don't need them for this particular fit. My source for the equation is Solid State Physics, Mermin &Ashcroft. I surely need to remove points under transitions, I made it in Matlab, but forgot it here. In data there aren't artifacts, as soon as I get home, I will try correcting my dataset. I guess I could use Pb values as starting guess, since my sample is a PbIn alloy. – iacolippo May 14 '15 at 15:53
  • I've tried with `leastsq` and the Jacobian but it says `too many values to unpack`. Even with starting guess very similar to the ones I found with Matlab. I don't think I'm requiring a so much difficult fit, probably I'm making some mistake as I'm a beginner with Python, but I guess I will improve with experience. – iacolippo May 14 '15 at 19:28
  • The returned values from `leastsq` are a bit different and more complicated than `curve_fit`: [the docs](http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.leastsq.html} give more details. – xnx May 14 '15 at 19:36
  • Are you able to post the whole of `salita.txt` and the initial guess parameters somewhere like pastebin.com? – xnx May 14 '15 at 19:39
  • Warning I get is: `/System/Library/Frameworks/Python.framework/Versions/2.7/Extras/lib/python/scipy/optimize/minpack.py:393: RuntimeWarning: Number of calls to function has reached maxfev = 800. warnings.warn(errors[info][0], RuntimeWarning) (array([ 9.18063198e-04, -1.82655707e-04, 9.56012732e+00]), 5) ` (returns last iteration). Initial guess parameters are: 10*10^(-6) for residual resistivity, 10*10^(-8) for electron phonon coupling AE, 50 for Einstein temperature. – iacolippo May 14 '15 at 20:16
  • My txt file is here [link](http://pastebin.com/Ew4tqz3E). It is made of 7 columns, first one is temperature, third is resistivity in this version (I have multiple version of it and they have different order in columns, not my fault...). Thank you – iacolippo May 14 '15 at 20:23
  • Hmmm.. your objective function is monotonically decreasing but the data are increasing (pretty much linearly): are you sure about your `einstein_func`? – xnx May 14 '15 at 20:42
  • Oh sure, this is the ascent of sample from dewar with liquid He, maybe it is necessary to invert vectors of data. Linear increase is correct since at high temperature resistivity is linear in temperature, there is a so-called knee only at intermediate temperatures, but here isn't that evident since there is superconductive transition. I have surely made stupid mistakes but no one ever explained me how to do a data fit. – iacolippo May 14 '15 at 20:58
  • I have made a stupid mistake: einstein_func isn't correct as you said. Correct function is: ro0+AE*T*((TE/2T)/sinh(TE/2T))^2... I feel so stupid right now. I'm trying this – iacolippo May 14 '15 at 21:05
  • I MADE IT! Thank you so much – iacolippo May 14 '15 at 21:08
  • @IacopoPoli Glad to help a fellow scientist! – xnx May 14 '15 at 21:10