1

我需要使用包括数值积分在内的复杂函数来拟合图形。如果我使用更简单的函数（这也是一个积分）来执行此操作，则此方法有效，但对于我正在尝试的函数，它不起作用。尽管该程序正在运行，但它非常不适合。您可以看到软件所做的调整是一条恒定线 y = 0 ......这不是预期的。

代码是：

import numpy as np
from scipy.integrate import quad
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt
import pandas as pd

dados=pd.read_excel("C:/Users/Samaung/Documents/Henrique/DadosExcel/dadosbalanca.xlsx")
t=dados["t"].values
m=dados["m"].values

mt=3.75*10**-5 # kg massa total empregada da amostra
rho=5961 #kg/m**3 densidade da amostra
v0=mt/rho #m**3 volume da amostra
uim=0.453 #A*m**2 momento magnético ímã
g=9.82 #m/s**2 aceleração da gravidade na Terra
z0=0.0063 #m distância centro ímã parte inferior da amostra
z1=0.0093 #m distância centro ímã partesuperior da amostra
R=0.0049 #m raio do porta amostra cilíndrico
dis=1/z0**4-1/z1**4-(z0**2+R**2)/(3*(z0**2+R**2)**3)+(z1**2+R**2)/(3*(z0**2+R**2)**3)
Msv=470644.784 #A/m magnetização volumétrica
Dm=17.6*10**-9 #m diâmetro médio
sigmad=0.13 #fator de dispersão adimensional
u0=1.256*10**-6 #N/A**2 permeabilidade magnética
H=70823.95 #A/m campo aplicado
Kb=1.3807*10**-23 #mm*2*kg/(s*K) constante de Boltzman
T=300 #K temperatura
Hk_=2/(u0*Msv) #campo de anisotropia dividido por Kef
t0=10**-9 #s fator de frequência
h_=H/Hk_
v1=6*Kb*T

def f(t,Kef):
    Multiplicador=3*u0*v0*uim/(32*np.pi*g)*dis*Msv*np.exp(2*sigmad**2)/(Dm*np.sqrt(2*np.pi)*sigmad)
    Integrando1=lambda D: np.exp(-(np.log(D/Dm))**2/(2*sigmad**2))
    Integrando2= lambda D: 1/np.tanh(u0*np.pi*Msv*H*D**3/v1)-v1/(u0*np.pi*Msv*H*D**3)
    Integrando3=lambda D: 1-np.exp(-t*((1-(h_/Kef)**2)*((1-h_/Kef)*np.exp(-Kef*np.pi*D**3/v1*(1-h_/Kef)**2)+(1+h_/Kef)*np.exp(-Kef*np.pi*D**3/v1*(1+h_/Kef)**2)))/(t0))
    fn=lambda D: Integrando1(D)*Integrando2(D)*Integrando3(D)
    ma=np.asarray(quad(fn,Dm/10,10*Dm)[0])
    return Multiplicador*ma
func=np.vectorize(f)

parametro, erro=curve_fit(func,t,m, p0=(139178.61676287447,), bounds=(0.1, np.inf))
print("Kef =", parametro[0])
xFit=np.arange(0.1,27105.06,1)
y = func(xFit, *parametro)
plt.plot(t, m, 'o')
plt.plot(xFit, y, '-')
plt.show()

我的数据是：dados= https://github.com/Henrique0501/Dadosbalanca/blob/main/dadosbalanca.xlsx

你能帮我解决这个问题吗？

4

0 回答 0