Fit an integral function with parametric limit to data with Python (Debye Model) -


i trying fit resistivity vs temperature data bloch-gruneisen formula resistivity in metals:

function

as can see there integral function parametric limit. don't know how implement algorithm run least squares fit. came with:

import matplotlib.pyplot plt import numpy np import pylab pl import scipy sp scipy.optimize import leastsq  #retrieve data file data = pl.loadtxt('salita.txt') temp = data[:, 0] res = data[:, 2]  def debye_func(p, t, r):     rho0, ad, td = p     coeff = ad*np.power(t, 5)/np.power(td, 4)     f = np.power(x^5)/np.power(np.sinh(x), 2) #function integrate     err_debye = r - rho0 - coeff * #integral???     return err_debye  p0 = sp.array([0.0001 , 0.00001, 50])  plsq = leastsq(debye_func, p0, args=(temp, res))  print plsq

ideas on how write it?

edit: code has become:

import matplotlib.pyplot plt import numpy np import pylab pl import scipy sp scipy.optimize import leastsq scipy.integrate import quad  #retrieve data file data = pl.loadtxt('salita.txt') temp = data[:, 0] res = data[:, 2]  def debye_integrand(x):     return np.power(x, 5)/np.power(np.sinh(x), 2)  def debye_func(p, t, r):     rho0, ad, td = p     coeff = ad*np.power(t, 5)/np.power(td, 4)     err_debye = r - rho0 - coeff * quad(debye_integrand, 0, td/(2*t))     return err_debye  p0 = sp.array([0.0001 , 0.00001, 50])  plsq = leastsq(debye_func, p0, args=(temp, res))  print plsq

now valueerror:

traceback (most recent call last):   file "debye.py", line 24, in <module>     plsq = leastsq(debye_func, p0, args=(temp, res))   file "/system/library/frameworks/python.framework/versions/2.7/extras/lib/python/scipy/optimize/minpack.py", line 348, in leastsq     m = _check_func('leastsq', 'func', func, x0, args, n)[0]   file "/system/library/frameworks/python.framework/versions/2.7/extras/lib/python/scipy/optimize/minpack.py", line 14, in _check_func     res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))   file "debye.py", line 19, in debye_func     err_debye = r - rho0 - coeff * quad(debye_integrand, 0, td/(2*t))   file "/system/library/frameworks/python.framework/versions/2.7/extras/lib/python/scipy/integrate/quadpack.py", line 247, in quad     retval = _quad(func,a,b,args,full_output,epsabs,epsrel,limit,points)   file "/system/library/frameworks/python.framework/versions/2.7/extras/lib/python/scipy/integrate/quadpack.py", line 296, in _quad     if (b != inf , != -inf): valueerror: truth value of array more 1 element ambiguous. use a.any() or a.all()

i think means i'm providing leastsq argument can't take, don't know how modify code.

edit2: solved equation analytically maxima , got

import matplotlib.pyplot plt import numpy np import pylab pl import scipy sp scipy.optimize import leastsq scipy.integrate import quad scipy.special import zetac mpmath import polylog  #retrieve data file data = pl.loadtxt('salita.txt') temp = data[:, 0] res = data[:, 2]  def debye_integrand(x):     return np.power(x, 5)/np.power(np.sinh(x), 2)  def debye_func(p, t, r, integral):     rho0, ad, td = p     coeff = ad*np.power(t, 5)/np.power(td, 4)     den = np.exp(td/t) -1     m1 = 5*((td/(2*t))**4)*np.log(np.exp(td/(2*t)+1)*(np.exp(td/t)-1)+120*polylog(5, np.exp(td/(t))*(1-np.exp(td/(2*t)))     m2 = 120*(td/(2*t))*polylog(4, np.exp(td/(2*t)))*(np.exp(np.exp(td/t))-1)+60*((td/(2*t))**2)*polylog(3, np.exp(td/(2*t))*(1-np.exp((td/(2*t)))     m3 = 20*((td/(2*t))**3)*polylog(2, np.exp(td/(2*t))*(np.exp(td/t)-1)+120**polylog(5, -np.exp(td/(2*t)))*(1-np.exp(td/t))     m4 = 120*(td/(2*t))*polylog(4, -np.exp(td/(2*t)))*(np.exp(td/t)-1)+60*((td/(2*t))**2)*polylog(3, -np.exp(td/(2*t)))*(1-np.exp(td/t))     m5 = 20*((td/(2*t))**3)*polylog(2, -np.exp(td/(2*t)))*(np.exp(td/t)-1) -2*((td/(2*t))**5)*np.exp(td/t)     m6 = 5*((td/(2*t))**4)*np.log(1-np.exp(td/(2*t))*(np.exp(td/t)-1)     zeta = 15.0*zetac(5)/2      integral = (m1+m2+m3+m4+m5+m6)/den +zeta      err_debye = r - rho0 - coeff * integral     return err_debye  #initizalied einstein model fit p0 = sp.array([0.00001 , 0.0000001, 70.0])    plsq = leastsq(debye_func, p0, args=(temp, res))  print plsq

it says syntaxerror: invalid syntax in m2. tried loops in numerical way, didn't succeed.

my .txt file here, if want try. first column temperature, third 1 resistivity.

you could, instance, separately define integrand function,

def debye_integrand(x,  n):         return x**n/((np.exp(x) - 1)*(1 - np.exp(-x)))

and use scipy.integrate.quad integration numerically,

from scipy.integrate import quad   def debye_func(p, t, r):     # [...] rest of code above here     err_debye = r - rho0 - coeff * quad(debye_integrand, 0, t/td, args=(n,))     return np.sum(err_debye**2)

that's general idea, , might need adapted further code. ideal solution find analytical solution integral, or rewrite classical integral functions scipy.special, but might not straightforward (see below).

also should use more general scipy.opitimize.minimize function instead of least-square fit, since provides algorithms more efficient , robust non-smooth optimizations. default optimization method bfgs start.

edit: actually, there analytical solution integral (for n=5), can obtain, instance, maxima,

>> integrate(x**5/((exp(x) - 1)*(1 - exp(-x))), x, 0, a)

solution part 1 enter image description here

where a integration limit, li_k polylogarithm function of order k (see mpmath.polylog) , ΞΆ riemann zeta function (see scipy.special.zetac).

although, depending on needs, might faster go numerical integration (or pre-calculated table lookup) rather puyting of together, , converting python.

edit 2: here final solution analytical calculation of integral,

import numpy np import mpmath mp scipy.optimize import minimize scipy.integrate import quad import matplotlib.pyplot plt  def debye_integral_sym_scalar(x):      """      calculate debye integral scalar using multi precision math,      otherwise overflows 64bit floats      """      exp_x = mp.exp(x)       m1 = -120*mp.polylog(5, exp_x)      m2 = 120*x*mp.polylog(4, exp_x)      m3 = -60*x**2*mp.polylog(3, exp_x)      m4 = 20*x**3*mp.polylog(2, exp_x)      m5 = 5*x**4*mp.log(1 - exp_x)      m6 = - x**5*exp_x       return m1 + m2 + m3 + m4 + m5 + m6/(exp_x - 1) + 120*mp.zeta(5)   # actual function can use  def debye_integral_sym(x):      f =  np.vectorize(debye_integral_sym_scalar, otypes=[np.complex])      return f(x).real  def debye_integrand(x,  n):      return x**n/((np.exp(x) - 1)*(1 - np.exp(-x)))  # test debye_integral_sym returns same result quad = 10.0  res0 =  quad(debye_integrand, 0, a, args=(5,))[0] res1 = debye_integral_sym(a) np.testing.assert_allclose(res0, res1)  def resistivity_fit(p, t):     rho0, ad, td = p     coeff = ad*np.power(t, 5)/np.power(td, 4)     return rho0 + coeff * debye_integral_sym(td/(2*t))   def debye_err_func(p, t, r):     return  np.sum((r - resistivity_fit(p, t))**2)  # wget "http://pastebin.com/raw.php?i=tvzcdxya" -o salita.txt data = np.loadtxt('salita.txt')  temp_exp = data[:, 0] res_exp = data[:, 2]  p0 = np.array([0.0001 , 0.00001, 50]) p_opt = minimize(debye_err_func, p0, args=(temp_exp, res_exp))  print p_opt  temp = np.linspace(temp_exp.min(), temp_exp.max(), 100)  plt.plot(temp_exp, res_exp, '.', label='experimental data') plt.plot(temp, resistivity_fit(p_opt.x, temp), 'r', label='bloch-gruneisen fit') plt.legend(loc='best') plt.xlabel('temperature [k]') plt.ylabel('resistivity') plt.show()

with output of optimization function,

   status: 0    success: true    njev: 5    nfev: 25    hess_inv: array([[  7.32764243e-01,  -4.89555962e-01,  -1.93879729e-08],    [ -4.89555962e-01,   3.27690582e-01,  -2.09510086e-08],    [ -1.93879729e-08,  -2.09510086e-08,   1.00000000e+00]])    fun: 1.784420370873494e-11    x: array([  9.96468440e-06,   7.40349389e-06,   5.00000000e+01])    message: 'optimization terminated successfully.'    jac: array([ -1.11880569e-06,   1.28115957e-06,   2.31303410e-12])

and resulting plot, the resulting plot


Comments

Post a Comment

Popular posts from this blog

c++ - Difference between pre and post decrement in recursive function argument -

i have following sample code used pre-decrement void function(int c) { if(c == 0) { return; } else { cout << "dp" << c << endl; function(--c); cout << c << endl; return; } } this function give output input 4 : dp3 dp2 dp1 dp0 0 1 2 3 but when use post decrement void function(int c) { if(c == 0) { return; } else { cout << "dp" << c << endl; function(c--); cout << c << endl; return; } } output goes in infinite loop dp4 can explain me in detail why happens ? this happens because function(c--) called same value of c , when finishes c decremented. recursively called, hence called same value no chance of return , hit stack overflow error. assume this, int x = 1, y = 0; y = x++ result y == 1 , x == 2. if y = ++x , bot x , y 2.
Read more

php - Nothing but 'run(); ' when browsing to my local project, how do I fix this? -

i added new project. made mistake in projectname.conf file in sites-available apache2, led me reinstalling apache2. after reinstalling , re-enabling projects on every page of every project .. <body> run(); </body> .. i've run composer update, restarted service, restarted pc, checked config files, tried disabling new project. there nothing in error log, , access.log: 127.0.0.1 - - [15/may/2015:01:40:47 +0200] "get / http/1.1" 304 181 "-" "mozilla/5.0 (x11; ubuntu; linux x86_64; rv:38.0) gecko/20100101 firefox/38.0" maybe i'm forgetting apache2 reinstall? don't think did different last time installed it. looks haven't setup php correctly. take @ this great guide , maybe you're forgetting link php config apache2.
Read more

php - How can I echo out this array? -

i don't know how go grabbing data in variable , put foreach loop. class tree{ private $info = array( array('name'=>'jim','state'=>'nu','company'=>'nu','phone'=>array('cell'=>'5615111111','office'=>'5611111111'),'email'=>array('primary'=>'exs@example.com','ex@there.com')), array('name'=>'joe smith','city'=>'phoenix','phone'=>'4805551111','email'=>'jsmith@some_email.com'), array('name'=>'john doe','city'=>'chandler','company'=>'doe co.','email'=>array('jdoe@gmail.com','personal'=>'email@email.com'),'phone'=>'6025550002') ); } if you're inside of class, can access private variable using $this->variablename . example, i...
Read more