1

I need to fit a graph using a complicated function including numerical integration. If I do this using a simpler function (which is also an integral), this method works, but for the function I'm trying, it's not working. Although the program is running, it is making a very bad fit. You can see that the adjustment the software makes is a constant line y = 0 ... this is not expected.

The code is:

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()

My data are: dados= https://github.com/Henrique0501/Dadosbalanca/blob/main/dadosbalanca.xlsx

Could you help me to fix this?

HENRIQUE MARTINS FERREIRA
  • 11
  • 4
  • The error message is quite clear. Can you plot the fit function with the initial guess? For readability I would clean up a bit more. `(...)**-1` would be a division. `cosh/sinh = 1/tanh` Moreover, the terms `6 * Kb * T` and `h_ / Kef` `Msv * H` and `D**3` appear several times. defining a variable for those would probably be a good idea. – mikuszefski Oct 22 '20 at 09:48
  • Is the `6**3` actually a `D**3` that would make more sense and avoid divergence as this the approaches something like `1/x - 1/x`....and the units of u0 in the comment are wrong – mikuszefski Oct 22 '20 at 11:01
  • also check your closing parenthesis in `integrando3` as for now there is an exponential of an exponential....is this supposed to be like this? – mikuszefski Oct 22 '20 at 11:05
  • I am also buffled by the interation limits. First of all: if `Dm` as units of meters it does not make sense to integrate from `1/Dm` to `10 * DM`. Next: if `Dm` is of order `1e-9` you are integrating from `10**9` to `10**(-8)` Is the reversed integration order intended? – mikuszefski Oct 22 '20 at 11:10
  • @mikuszefski Thank you very much for your comments. I adhered to some of your suggestions and I think it has improved a lot. The ideal is not to use initial guess, but I decided to use p0 = 1, because if I don't specify p0, the error `ValueError: Unable to determine number of fit parameters.` happens. But I don't want to use p0 = 1. Do you know how I could fix this? In fact, the integration limits are Dm / 10 and Dm * 10. You are right. – HENRIQUE MARTINS FERREIRA Oct 22 '20 at 16:09
  • And, really, in integrando3 the function has exponentials within an exponential. I put a picture of what the function is like. The meanings of h, K and x are `h = h_ / Kef`, `K=Kef` and `x = D` – HENRIQUE MARTINS FERREIRA Oct 22 '20 at 16:09
  • Of course, Kef is not 1. The console returned Kef=1 because of p0=1. And I'm so upset by this RuntimeWarning – HENRIQUE MARTINS FERREIRA Oct 22 '20 at 16:18
  • Ok...sounds good. Can you update your code in the post accordingly so everyone is on the same page? Do not put is as an answer! Also the `6**3` is a `D**3`, right (otherwise units would go out the window). That's for a vibrating sample magnetometer? – mikuszefski Oct 23 '20 at 06:42
  • Can you provide a reference for the formula?....I still don't get the reason for the `(...)**(-1)` in the denominator of the last exponential. That just puts it into the numerator. Why not writing it like this? – mikuszefski Oct 23 '20 at 07:05
  • What would the unit of the pre-factor be? What is the physical meaning? – mikuszefski Oct 23 '20 at 08:53
  • `Dm` is of order nanometer? Is this correct? – mikuszefski Oct 23 '20 at 13:41
  • @mikuszefski I did the changes you sugested. Thank you. The program is better now, but you can see that the fit is not good yet. Could you help me? – HENRIQUE MARTINS FERREIRA Nov 05 '20 at 00:28
  • This function describes the apparent mass registered by a scale whose plate rises over time due to the magnet on it being attracted by a sample of magnetic nanoparticles. – HENRIQUE MARTINS FERREIRA Nov 05 '20 at 00:30
  • Hi, I can have a look, but have limited time. First question. Can you check the units of your factor. I get Kg/m. Ids this correct?...the rest is unit free – mikuszefski Nov 05 '20 at 10:41
  • @mikuszefski Which factor are you talking about? The factor t0 is a frequency factor (proportional to the period of attempts to jump the barrier by the magnetic moment ... the unit is s). The sigmad factor is dimensionless (a dispersion factor for Dm measurements) – HENRIQUE MARTINS FERREIRA Nov 06 '20 at 00:35
  • Are you talking about Multiplicador? I checked the unit. In fact, it is kg/m – HENRIQUE MARTINS FERREIRA Nov 06 '20 at 00:45
  • ok...so it is a mass-linedensity? why not. The value is huge though. Can you check if this makes sense? – mikuszefski Nov 06 '20 at 06:42
  • Ok...had some time to check the functions in the integrator. First one is a peak that basically is between 0 and 3 `Dm`. The last function only makes some sense if `Kef` is of order `h_`. with `u0 Msv H = edZ` beeing the Zeeman energy density, we have `h_ = edZ / 2`. Anyhow to see something we have to put approx 15000 for `Kef`. Then to see a behavior in time we have to look for `0 < t < 1e-8` if your data is mass vs time there is something wrong by orders of magnitude. (If you are still interested in sorting this out, tell me. Otherwise I leave this as is.) – mikuszefski Nov 12 '20 at 09:28

0 Answers0