0

I would like to fit an hysteresis curve, with a superparamagnetic behavior, using a magnetic model which includes a Langevin function and a pair distribution function [1]. To fit with this equation I must solve a definite integral. I was trying to use scipy.integrate.quad for this purpose and the features of lmfit, but I do not get -at least- to simulate a reasonable curve (see the code below). Some real physical parameters that can be used to simulate this equation are: Dm = 3.2E-9 m, w = 0.26, NT = 1.7E12, kB=1.38E-23 J/K and T=300K. Simulating with this values it must result in a superparamagnetic curve as the contained in the link below. I would appreciate any suggestion to make this code work and to improve it.

Langevin equation with PDF

www.dropbox.com/pri/get/superpara.dat?_subject_uid=197016565&w=AADkHqW1w-gE9pQkG0oLoE7tNG1J-rWxN0lcIM9ioXWiLA

from lmfit import minimize, Minimizer, Parameters, Parameter, report_fit
from numpy import loadtxt, vectorize, sqrt, log, log10, inf, exp, pi, tanh
from scipy.integrate import quad
import matplotlib.pyplot as plt

dat = loadtxt('superpara.dat')

u0_H   = dat[:, 0]
dat1 = dat[:, 1]

def PDF(Dmag, u0_H, params):
    v = params.valuesdict()
    pdf = (1/(v['w']*sqrt(2*pi)*Dmag))*exp(-(log(Dmag/v['Dm']))**2 / (2*v['w']**2))
    x = (pi/(6*v['kB']*v['T']))*v['Ms']*(u0_H*1e-4)*Dmag**3
    return (pi/6)*v['Ms']*(Dmag**3)*( (1/tanh(x))-(1/x) )*pdf

def curve(u0_H, params):
    return params['NT']*quad(PDF, 0.0, inf, args=(u0_H, params))[0]

vcurve = vectorize(curve, excluded=set([1]))

def fit_function(params, u0_H, dat1):
    model1 = vcurve(u0_H, params)
    resid1 = dat1 - model1
    return resid1.flatten()

params = Parameters()
params.add('Dm' , value= 3.2e-9  , vary= True)
params.add('w'  , value= 0.26    , vary= True)
params.add('Ms' , value= 0.1     , vary= True)
params.add('T'  , value= 300     , vary= False)
params.add('kB' , value= 1.38e-23, vary= False)
params.add('NT' , value= 1.7e12  , vary= True)

minner = Minimizer(fit_function, params, fcn_args=(u0_H, dat1))
result = minner.minimize()

report_fit(result)

y1_fit = vcurve(u0_H, result.params)
y1_init = vcurve(u0_H, params)

plt.plot(u0_H, dat1, 'k+', u0_H, y1_init, 'b-', u0_H, y1_fit, 'r-')
plt.show()
Ernesto
  • 21
  • 4
  • 1
    it would be helpful to know what isn't working, including output and any error messages, and what you expected. If you call your `vcurve()` with reasonable parameters, do you get something resembling your data? – M Newville Jun 11 '17 at 03:18
  • It does not show any error messages, and the initial values do not change after the fit. Here the output: Dm: 3.2000e-09 +/- 0 (0.00%) (init= 3.2e-09) w: 0.26000000 +/- 0 (0.00%) (init= 0.26) Ms: 0.10000000 +/- 0 (0.00%) (init= 0.1) T: 300 (fixed) kB: 1.38e-23 (fixed) NT: 1.7000e+12 +/- 0 (0.00%) (init= 1.7e+12) The initial parameters have physical sense and I do not get something resembling my data! I have perceived that when I give a value close to unit for Dm, kB and NT, It shows a fit result but with a very high error for some parameters. – Ernesto Jun 11 '17 at 15:22
  • 1
    I would suggest that you first verify that calling your `curve()` works as you expect (giving reasonable output resembling the data) and then verify that calling your `vcurve()` works as you expect. Only after `curve()` turns your input into the correct output should you expect the fitting methods to be able to refine your parameter values. – M Newville Jun 11 '17 at 19:09
  • @MNewville I tried to call the vcurve() (I cannot call the curve() due to the need to vectorize after use the quad function). It does not work! Could it be a problem related to the inability of scipy.integrate.quad to solve this function? – Ernesto Jun 12 '17 at 17:02
  • 1
    Well, if `vcurve()` does not work, then curve-fitting using this as the model definitely won't work. But you know where to look next! And, yes the probably could definitely be related to `integrate.quad`. I do not see anything obviously wrong, but your function looks sort of complicated.... – M Newville Jun 15 '17 at 01:53
  • It works after rescaling Dmag by a factor of 1e-9. It seems that is necessary to be careful with very high potentials such as those introduced by kB into the equation. – Ernesto Jun 25 '17 at 22:11

0 Answers0