使用四边形 - 磁性磁滞循环拟合数据

Question

我很难获得融合，因为它不能收敛或给出NAN错误，具体取决于我的开始参数。我正在使用Quad使用LMFIT集成和拟合。任何帮助都将受到赞赏。

我将数据拟合到langevin功能，并由对数正态分布加权。由于我的声誉得分，Stackoverflow不允许我发布该功能的图像，但它在下面的代码中。

我正在插入MS，DM和Sigma的H（字段）和拟合，而MU_0，MSB，KB和T都是常数。

这是我使用的内容，使用一些示例数据：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy
from numpy import vectorize, sqrt, log, inf, exp, pi, tanh
from scipy.constants import k, mu_0
from lmfit import Parameters
from scipy.integrate import quad

x_data = [-7.0, -6.5, -6.0, -5.5, -5.0, -4.5, -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, 
          -0.95, -0.9, -0.85, -0.8, -0.75, -0.7, -0.65, -0.6, -0.55, -0.5, -0.45, -0.4, 
          -0.35, -0.3, -0.25, -0.2, -0.1,-0.05, 3e-6, 0.05, 0.1, 0.15, 0.2, 0.25, 0.3, 
          0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 
          1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 5.5, 6.0, 6.5, 7.0]
y_data = [-61.6, -61.6, -61.6, -61.5, -61.5, -61.4, -61.3, -61.2, -61.1, -61.0, -60.8, 
          -60.4, -59.8, -59.8, -59.7, -59.5, -59.4, -59.3, -59.1, -58.9, -58.7, -58.4, 
          -58.1, -57.7, -57.2, -56.5, -55.6, -54.3, -52.2, -48.7, -41.8, -27.3, 2.6, 
          30.1, 43.1, 49.3, 52.6, 54.5, 55.8, 56.6, 57.3, 57.8, 58.2, 58.5, 58.7, 59.0, 
          59.1, 59.3, 59.5, 59.6, 59.7, 59.8, 59.9, 60.5, 60.8, 61.0, 61.2, 61.3, 61.4, 
          61.4, 61.5, 61.6, 61.6, 61.7, 61.7]

params = Parameters()
params.add('Dm' , value = 8e-9   , vary = True, min = 0, max = 1) # magnetic diameter (m)
params.add('s'  , value = 0.4    , vary = True, min = 0.0, max = 10.0) # sigma, unitless
params.add('Ms' , value = 61.0    , vary = True) #, min = 30.0 , max = 100.0) # saturation magnetization (emu/g)

params.add('Msb', value = 446000 * 1e-16, vary = False) # Bulk magnetite saturation magnetization (A/m)
params.add('T'  , value = 300    , vary = False) # Temperature (K)

def Mag(x_data, params):
    
    v = params.valuesdict() # put parameters into a dictionary

    def numerator(D, x_data, params):
        
        # langevin 
        a_numerator = pi * v['Msb'] * x_data * D**3
        a_denominator = 6*k*v['T']
        a = a_numerator / a_denominator
        langevin = (1/tanh(a)) - (1/a)
        
        # PDF
        exp_num = (log(D/v['Dm']))**2
        exp_denom = 2 * v['s']
        exponential = exp(-exp_num/exp_denom)
        pdf = exponential/(sqrt(2*pi) * v['s'] * D)
        
        return D**3 * langevin * pdf
    
    def denominator(D, params):
        
        # PDF
        exp_num = (log(D/v['Dm']))**2
        exp_denom = 2 * v['s']
        exponential = exp(-exp_num/exp_denom)
        pdf = exponential/(sqrt(2*pi) * v['s'] * D)
        return D**3 * pdf
    
    # return integrals
    return v['Ms'] * quad(numerator, 0, inf, args=(x_data, params))[0]  / quad(denominator, 0, inf,args=(params))[0]

# vectorize
vcurve = np.vectorize(Mag, excluded=set([1]))

plt.plot(x_data, vcurve(x_data, params))
plt.scatter(x_data, y_data)

这绘制了带有启动参数的数据和拟合方程。我在Langevin中的某个地方有一个

from lmfit import minimize, Minimizer, Parameters, Parameter, report_fit

def fit_function(params, x_data, y_data):
    model1 = vcurve(x_data, params)
    resid1 = y_data - model1
    return resid1

minner = Minimizer(fit_function, params, fcn_args=(x_data, y_data))
result = minner.minimize()

report_fit(result)
result.params.pretty_print()

问题，积分不会收敛，给出以下错误：

/var/folders/pz/tbd_dths0_512bm6l43vpg680000gp/T/ipykernel_68003/1413445460.py:39: IntegrationWarning: The algorithm does not converge.  Roundoff error is detected
  in the extrapolation table.  It is assumed that the requested tolerance
  cannot be achieved, and that the returned result (if full_output = 1) is 
  the best which can be obtained.
  return v['Ms'] * quad(numerator, 0, inf, args=(x_data, params))[0]  / quad(denominator, 0, inf,args=(params))[0]

我坚持为什么拟合不收敛。这是一个问题，因为我使用的数字很小，还是Quad/LMFit的问题？谢谢你！

Answer 1

具有更接近订单1的参数（例如，1E-7和1E7之间）是一个好主意。如果您期望参数在1.e-9（或1.E-16！）范围内，则肯定可以将其缩放（在拟合功能中），以便拟合算法来回传递的值更接近订单1。但是，我有点怀疑这是您遇到的主要问题。

在我看来，您的mag函数对变量参数的值dm和s都不是很敏感的。我不是100％确定为什么这样。您是否使用“ mag”或“ vcurve”验证了计算，您期望他们做的事情？