# Fitting a numerical integral via NonLinearModelFit to magnetic data

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 )

-
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 '14 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 '14 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 '14 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 '14 at 19:41
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 '14 at 6:39

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}]


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/(#)&

-
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 '14 at 9:24