0

我想使用包括朗之万函数和对分布函数 [ 1 ]的磁模型来拟合具有超顺磁行为的磁滞曲线。为了适应这个方程，我必须求解一个定积分。我试图为此目的使用 scipy.integrate.quad 和 lmfit 的功能，但我没有得到 - 至少 - 来模拟合理的曲线（参见下面的代码）。可以用来模拟这个方程的一些实际物理参数是：Dm = 3.2E-9 m，w = 0.26，NT = 1.7E12，kB=1.38E-23 J/K 和 T=300K。使用此值进行模拟必须产生超顺磁曲线，如下面的链接中包含的那样。我将不胜感激任何使此代码工作并改进它的建议。

朗之万方程与 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()
4

0 回答 0