Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm new to Mathematica and I'm currently trying to fit
$$m_T (H,T) = N_T \int\limits_0^{\infty} \frac{x k_\text{B} T}{\mu_0 H} \mathcal{L}(x) \text{pdf}(D_\text{mag}) \text{d}D_\text{mag}$$
with $x=\mu_0 M_S \pi D_\text{mag}^3 H/(6k_\text{B}T)$ and $\mathcal{L}(x)=\coth(x)-1/x$ (Langevin function)
to experimental data (this is a fit for the total magnetic moment of superparamagnetic nanoparticles).

So I tried the following:

kB = 1.38*10^-23; (*Boltzmann constant*)

pdf[d_?NumericQ, d0_?NumericQ, w_?NumericQ] = 
  1/(w*Sqrt[2 Pi])*1/d*E^(-1/(2*w^2)*(Log[d/d0])^2);

mt[b_?NumericQ, Nt_?NumericQ, Ms_?NumericQ, d0_?NumericQ, w_?NumericQ] := 
  Nt*NIntegrate[(Ms*Pi*d^3/6)*(Coth[b*Ms*Pi*d^3/(6*kB*300)] - 
    1/(b*Ms*Pi*d^3/(6*kB*300)))*pdf[d, d0, w], {d, 0, 1}];

As I think starting values are quite important here, I tried to find reasonable values and ended with

The blue dots are the experimental data and the red line is 1000*mt[0.0001*b, 0.35*10^14, 20150, 2.0246171576454905*10^-8, 0.3].

Note 1: The 1000*... and 0.0001*b are due to conversion of cgs units <-> SI units.

Note 2: I know that some of the nanoparticles exhibit ferromagnetic behaviour, thus the hysteresis, but for now I would be happy to just fit the superparamagnetic total magnetic moment to these points.

So I tried to fit mt to the data:

nlmB1 = 
  NonlinearModelFit[data, 
    {mt[b, Nt, Ms, d0, w], 
      {10^(15) > Nt > 10^(13) && 20200 > Ms > 20000 && 5*10^(-8) > d0 > 1*10^(-8) && 0.5 > w > 0.29}}, 
    {{Nt, 0.35*10^14}, {Ms, 20150}, {d0, 2.0246171576454905*10^-8}, {w, 0.3}}, 
    b]

But I get an error saying that There are no points that satisfy the constraints {}. So I removed most of the constraints and tried

nlmB1 = 
  NonlinearModelFit[data, 
    {mt[b, Nt, Ms, d0, w], {10^(15) > Nt && Nt > 10^(13)}}, 
    {{Nt, 0.35*10^14}, {Ms, 20150}, {d0, 2.0246171576454905*10^-8}, {w, 0.3}}, 
    b]

This time the fitting procedure takes forever and mathematica tells me

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in d near {d} = {2.07884*10^-8}. NIntegrate obtained 8.417808128264057*^-20 and 1.8814309233927094*^-16 for the integral and error estimates.
General::stop: "Further output of NIntegrate will be suppressed during this calculation."

What do I need to do to get Mathematica really fit this? I have also tried several Methods and I've changed the AccuracyGoal without getting a successful fit.


The data is included via

dataInput = Import["Measurement.dat", "Table", "FieldSeparators" -> ","];
data = dataInput[[All, {2, 4}]];

(If you want to try, here is the data: https://www.dropbox.com/s/z55iktrt0poki5y/Measurement.dat )

share|improve this question
    
Can you recheck your post, you write 1000*mt[0.0001*b, 300, 0.35*10^14, 20150, 2.0246171576454905*10^-8, 0.3] however it looks like mt only takes 5 arguments. –  bobthechemist Apr 2 at 16:59
    
@bobthechemist You're right, thanks, the 300 was from the temperature which I had as a variable in the beginning. Post updated. –  John Apr 2 at 17:09
    
the "no points satisfy" error comes from starting with w=.3 which conflicts with w>.3. By the way you can do one statement eg 20200 > Ms > 20000 instead of 20200 > Ms && Ms > 20000. Its not wrong as you have it but its a bit cumbersome to read. –  george2079 Apr 2 at 18:28
    
@george2079 I'll edit it zomorrow. I thought that, due to the numerical nature of the fitting process all strictly greater or lesser than are converted to greater or equal / lesser or equal, I may be wrong though. If I recall correctly it didn't work with 0.29, but I can't check atm, as I have no pc at hands. And thanks so far! –  John Apr 2 at 19:41
2  
Usually physicists do 2 steps to make things much easier. (1) rescale to dimensionless units - kB for example shouldn't be in the picture at all. (2) rescale to characteristic system scales so numbers like {d0, 2.0246171576454905*10^-8} become something like {d0, 2.02} where d0 is measured in some sort of a new scale-unit. Not only you avoid things like precision loss and messy representation - you also gain a better understanding of the physics of the problem. Characteristic scales are different between astronomy and solid state - they bring clarity to your fitting parameters. –  Vitaliy Kaurov Apr 3 at 6:39

1 Answer 1

up vote 3 down vote accepted

Following @Vitaliy's comment try this formulation:

kB = QuantityMagnitude@
    UnitConvert[Quantity["BoltzmannConstant"] , "Joule/Kelvin"]

pdf[dp_?NumericQ, d0p_?NumericQ, w_?NumericQ] = 
    kB^(-1/3)/(w Sqrt[2 Pi]) 1/dp *E^(-1/(2  w^2) (Log[dp /d0p])^2);
mt[b_?NumericQ, Nt_?NumericQ, Ms_?NumericQ, d0p_?NumericQ, w_?NumericQ] := 
    Nt kB^(4/3) Ms Pi /6  NIntegrate[  
        dp^3 (Coth[#] - 1/(#) &@(b Ms Pi dp^3/(6 300))) pdf[dp  , d0p, w],
          {dp, 0, Infinity}]

Note I've replaced d,d0 with dp,d0p which are normalized by kB^(1/3) (You should verify!!) Here is your fit plot, with the new d0p argument .844 ~ 2E-8/kB^(1/3)

Plot[ 1000*mt[0.0001*b, 0.35*10^14, 20150, .844 , 0.3] , {b, -10000, 10000}]

enter image description here

this could use some more work , but the essential issue here is your integrand was changing very rapidly for very small values of d. (The fact that setting the Infinite integration limit to 1 worked is a clue to that issue..). I didn't try the model fit but it should work much better. (note the integration limit now has to be infinity)

another thing..

that Coth[x]-1/x may be causing errors for small x (even though the limit behavior is well behaved and is x/3 + O(x^3). Try this:

Piecewise[{{#/3, Abs[#] < 10^-5}, {Coth[#] - 1/(#), True}}]&

in place of:

Coth[#] - 1/(#)&
share|improve this answer
    
Thanks a lot, this actually worked, even though fitting takes like ~24 hours and doesn't reach the AccuracyGoal, the result is quite sensible. –  John Apr 10 at 9:24

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.