2
$\begingroup$

The data to be fitted is for y vs x. We have y and x points. We do not have data for w. In fact, w variable is being varied during integration from 1e-9 to 100e-9. The fitting parameters are Ms,sigma,delta. For clarity I have uploaded the image showing whole equation in simplified form. Please see.

enter image description here

Part of the data is here.

data={{0.0049351, 887.55}, {0.014628, 2076.6}, {0.024377, 
  2684.6}, {0.034198, 3044.85}, {0.043943, 3281.3}, {0.053758, 
  3454.15}, {0.06349, 3585.85}, {0.073305, 3692.4}, {0.083048, 
  3778.2}, {0.092853, 3851.5}, {0.10261, 3914.95}, {0.11237, 
  3962.55}, {0.12216, 4006.85}, {0.13189, 4049.45}, {0.1417, 
  4089.2}, {0.15144, 4125.1}, {0.16124, 4157}, {0.17096, 
  4184.6}, {0.18077, 4211.9}, {0.1905, 4236.95}, {0.2003, 
  4260.15}, {0.21004, 4284.15}, {0.21982, 4302.95}, {0.22961, 
  4322.25}, {0.23933, 4340.35}, {0.24913, 4357.45}, {0.25886, 
  4375.4}, {0.26866, 4388.4}, {0.27837, 4405}, {0.28817, 
  4418.3}, {0.29791, 4434.5}, {0.30772, 4446.2}, {0.31745, 
  4457.3}, {0.32727, 4470.95}, {0.33704, 4483.4}, {0.34678, 
  4492.65}, {0.35659, 4499.4}, {0.36632, 4512.6}, {0.37614, 
  4523.2}, {0.38588, 4532.3}, {0.39568, 4541.95}, {0.40542, 
  4552.5}, {0.41522, 4557.2}, {0.42495, 4570.85}, {0.43473, 
  4574.45}, {0.44451, 4583.85}, {0.45423, 4590.65}, {0.46403, 
  4604.35}, {0.47377, 4609.8}, {0.48356, 4616.9}, {0.49329, 
  4624.6}, {0.49728, 4624.55}}

The following code with errors is after assistance by members and some of my effort after reading previous posts. It is better than earlier but not enough to bring results. I feel something is wrong in the definition of fitting model. Please help.

data = Import["C:/Users/aj/Desktop/LG once again/Mathematica/lgdata3.csv", "CSV", "HeaderLines" -> 7];

k = 1.38064852*10^-23;
T = 298;

model[Ms_?NumericQ, \[Sigma]_?NumericQ, \[Delta]_?NumericQ] := NIntegrate[Ms (Coth[(Ms x (1/6 Pi w^3))/(k *T)] - ((k*T)/(Ms x (1/6 Pi w^3)))) *(1/(Sqrt[2 Pi] w \[Sigma]) ) Exp[-{Log[w/\[Delta]]}^2/(2 \[Sigma]^2)], {w, 1*10^-10, 1*10^-7}]

fit = FindFit[data,model, {{Ms, 5500}, {\[Sigma], 2/10}, {\[Delta], 10*10^(-9)}}, x,MaxIterations -> 1000, Method -> NMinimize,NormFunction -> (Norm[#, Infinity] &), Gradient -> "FiniteDifference"]

FindFit::nrnum: The function value Max[Abs[-4624.6+model],Abs[-4624.55+model],Abs[-4616.9+model],Abs[-4609.8+model],Abs[-4604.35+model],Abs[-4590.65+model],Abs[-4583.85+model],Abs[-4574.45+model],<<35>>,Abs[-3778.2+model],Abs[-3692.4+model],Abs[-3585.85+model],Abs[-3454.15+model],Abs[-3281.3+model],Abs[-3044.85+model],Abs[-2684.6+model],<<2>>] is not a real number at {Ms,\[Sigma],\[Delta]} = {5500.,0.2,1.*10^-8}. >>

NMinimize::nnum: The function value Experimental`NumericalFunction[{Hold[(Norm[#1,\[Infinity]]&)[{-887.55+model,-2076.6+model,-2684.6+model,-3044.85+model,-3281.3+model,-3454.15+model,-3585.85+model,-3692.4+model,-3778.2+model,<<34>>,-4570.85+model,-4574.45+model,-4583.85+model,-4590.65+model,-4604.35+model,-4609.8+model,-4616.9+model,<<2>>}]],Block},<<4>>,{<<1>>}][{<<1>>}] is not a number at {Ms,\[Delta],\[Sigma]} = {5500.,1.*10^-8,0.2}. >>

FindFit::nrnum: The function value Max[Abs[-4624.6+model],Abs[-4624.55+model],Abs[-4616.9+model],Abs[-4609.8+model],Abs[-4604.35+model],Abs[-4590.65+model],Abs[-4583.85+model],Abs[-4574.45+model],<<35>>,Abs[-3778.2+model],Abs[-3692.4+model],Abs[-3585.85+model],Abs[-3454.15+model],Abs[-3281.3+model],Abs[-3044.85+model],Abs[-2684.6+model],<<2>>] is not a real number at {Ms,\[Sigma],\[Delta]} = {5500.,0.2,1.*10^-8}. >>

NMinimize::nnum: The function value Experimental`NumericalFunction[{Hold[(Norm[#1,\[Infinity]]&)[{-887.55+model,-2076.6+model,-2684.6+model,-3044.85+model,-3281.3+model,-3454.15+model,-3585.85+model,-3692.4+model,-3778.2+model,<<34>>,-4570.85+model,-4574.45+model,-4583.85+model,-4590.65+model,-4604.35+model,-4609.8+model,-4616.9+model,<<2>>}]],Block},<<4>>,{<<1>>}][{<<1>>}] is not a number at {Ms,\[Delta],\[Sigma]} = {5500.,1.*10^-8,0.2}. >>

FindFit::nrnum: The function value Max[Abs[-4624.6+model],Abs[-4624.55+model],Abs[-4616.9+model],Abs[-4609.8+model],Abs[-4604.35+model],Abs[-4590.65+model],Abs[-4583.85+model],Abs[-4574.45+model],<<35>>,Abs[-3778.2+model],Abs[-3692.4+model],Abs[-3585.85+model],Abs[-3454.15+model],Abs[-3281.3+model],Abs[-3044.85+model],Abs[-2684.6+model],<<2>>] is not a real number at {Ms,\[Sigma],\[Delta]} = {5500.,0.2,1.*10^-8}. >>

General::stop: Further output of FindFit::nrnum will be suppressed during this calculation. >>

NMinimize::nnum: The function value Experimental`NumericalFunction[{Hold[(Norm[#1,\[Infinity]]&)[{-887.55+model,-2076.6+model,-2684.6+model,-3044.85+model,-3281.3+model,-3454.15+model,-3585.85+model,-3692.4+model,-3778.2+model,<<34>>,-4570.85+model,-4574.45+model,-4583.85+model,-4590.65+model,-4604.35+model,-4609.8+model,-4616.9+model,<<2>>}]],Block},<<4>>,{<<1>>}][{<<1>>}] is not a number at {Ms,\[Delta],\[Sigma]} = {5500.,1.*10^-8,0.2}. >>

General::stop: Further output of NMinimize::nnum will be suppressed during this calculation. >>


Plot[Evaluate[model[Ms, \[Sigma], \[Delta], x] /. fit], {x, 5*^-3, 
0.5}, Epilog -> Point[data], PlotRange -> All]

And then further list of errors is long. That's why I am putting here just its initial part.

During evaluation of In[6]:= FindFit::nrnum: The function value Max[Abs[-4624.6+model],Abs[-4624.55+model],Abs[-4616.9+model],Abs[-4609.8+model],Abs[-4604.35+model],Abs[-4590.65+model],Abs[-4583.85+model],Abs[-4574.45+model],<<35>>,Abs[-3778.2+model],Abs[-3692.4+model],Abs[-3585.85+model],Abs[-3454.15+model],Abs[-3281.3+model],Abs[-3044.85+model],Abs[-2684.6+model],<<2>>] is not a real number at {Ms,[Sigma],[Delta]} = {5500.,0.2,1.*10^-8}.

>

During evaluation of In[6]:= NMinimize::nnum: The function value Experimental`NumericalFunction[{Hold[(Norm[#1,[Infinity]]&)[{-887.55+model,-2076.6+model,-2684.6+model,-3044.85+model,-3281.3+model,-3454.15+model,-3585.85+model,-3692.4+model,-3778.2+model,<<34>>,-4570.85+model,-4574.45+model,-4583.85+model,-4590.65+model,-4604.35+model,-4609.8+model,-4616.9+model,<<2>>}]],Block},<<4>>,{<<1>>}][{<<1>>}] is not a number at {Ms,[Delta],[Sigma]} = {5500.,1.*10^-8,0.2}. >>

$\endgroup$
  • 2
    $\begingroup$ I can see three immediate problems. 1) D is a reserved symbol in Mathematica and can't be used as a variable. 2) α, lang, lognorm, and z are not defined as functions of their variables. 3) sdata = data; model = z; are useless. However, I think the error message are being caused by the formatting of your CSV file. $\endgroup$ – m_goldberg Jul 19 '16 at 7:45
  • $\begingroup$ I have modified program with all three corrections, as you suggested. I have also enabled output of data import. But still the same error appears. However, the program works fine and gives fitted parameters (Ms and w) when we remove integration and lognormal part from the model. When we include integration and lognormal and make w as 'dw' for integration then outcome is infeasible. $\endgroup$ – ajay Jul 19 '16 at 10:09
  • $\begingroup$ thank you m_goldberg.Please look once again into the problem.If possible, please give me your email-id , I will send my calculated notebook.my email is ajayshankar0(AT)gmail(DOT)com $\endgroup$ – ajay Jul 19 '16 at 10:10
  • $\begingroup$ Here, is a glimpse of the mathematica output when I import data in the notebook. {{0.0049351, 887.55}, {0.014628, 2076.6}, {0.024377, 2684.6}, {0.034198, 3044.85}, {0.043943, 3281.3}, {0.053758, 3454.15}, {0.06349, 3585.85}, {0.073305, 3692.4}, {0.083048, 3778.2}, {0.092853, 3851.5}, {0.10261, 3914.95}, {0.11237, 3962.55}, {0.12216, 4006.85}, {0.13189, 4049.45}, {0.1417, 4089.2}, {0.15144, 4125.1}, {0.16124, 4157}, {0.17096, 4184.6}, {0.18077, 4211.9}, {0.1905, 4236.95}, {0.2003, 4260.15}, {0.21004, 4284.15}, {0.21982, 4302.95}, $\endgroup$ – ajay Jul 19 '16 at 10:20
  • $\begingroup$ If you have made the corrections I suggested, you should also edit your question and make them here. $\endgroup$ – m_goldberg Jul 19 '16 at 14:33
5
$\begingroup$

Use of NIntegrate and ?NumericQ were helpful in reducing errors. Take a look at the referenced questions at the bottom for additional information.

ClearAll["Global`*"]

data = {{0.0049351, 887.55}, {0.014628, 2076.6}, {0.024377, 
    2684.6}, {0.034198, 3044.85}, {0.043943, 3281.3}, {0.053758, 
    3454.15}, {0.06349, 3585.85}, {0.073305, 3692.4}, {0.083048, 
    3778.2}, {0.092853, 3851.5}, {0.10261, 3914.95}, {0.11237, 
    3962.55}, {0.12216, 4006.85}, {0.13189, 4049.45}, {0.1417, 
    4089.2}, {0.15144, 4125.1}, {0.16124, 4157}, {0.17096, 
    4184.6}, {0.18077, 4211.9}, {0.1905, 4236.95}, {0.2003, 
    4260.15}, {0.21004, 4284.15}, {0.21982, 4302.95}, {0.22961, 
    4322.25}, {0.23933, 4340.35}, {0.24913, 4357.45}, {0.25886, 
    4375.4}, {0.26866, 4388.4}, {0.27837, 4405}, {0.28817, 
    4418.3}, {0.29791, 4434.5}, {0.30772, 4446.2}, {0.31745, 
    4457.3}, {0.32727, 4470.95}, {0.33704, 4483.4}, {0.34678, 
    4492.65}, {0.35659, 4499.4}, {0.36632, 4512.6}, {0.37614, 
    4523.2}, {0.38588, 4532.3}, {0.39568, 4541.95}, {0.40542, 
    4552.5}, {0.41522, 4557.2}, {0.42495, 4570.85}, {0.43473, 
    4574.45}, {0.44451, 4583.85}, {0.45423, 4590.65}, {0.46403, 
    4604.35}, {0.47377, 4609.8}, {0.48356, 4616.9}, {0.49329, 
    4624.6}, {0.49728, 4624.55}};

k = 1.38064852*10^-23;
T = 298;

δ = 5.2*^-8;

model[Ms_?NumericQ, σ_?NumericQ, x_?NumericQ] := 
  NIntegrate[ Ms (Coth[(Ms x (1/6 Pi w^3))/(k T)] - 
  ((k T)/(Ms x (1/ 6 Pi w^3)))) (1/(Sqrt[2 Pi] w σ)) 
  Exp[-(Log[ w/δ])^2/(2 σ^2)], {w, 10*^-10, 10*^-7}];

fit = FindFit[data, model[Ms, σ, x], {{Ms, 5000}, {σ, 0.5}}, x]

Plot[Evaluate[model[Ms, σ, x] /. fit], {x, 0.005, 0.5}, 
 Epilog -> Point[data], PlotRange -> All, AxesOrigin -> {0, 0}]

enter image description here

{Ms -> 4811.78, σ -> 0.397626}

As requested by the OP, here is a table of fitted values:

Table[{x, First@Evaluate[model[Ms, σ, x] /. fit]}, {x, 0.005, 0.5, 0.01}]

{{0.005, 1040.15}, {0.015, 2064.21}, {0.025, 2614.87}, {0.035, 2970.88}, {0.045, 3223.37}, {0.055, 3413.07}, {0.065, 3561.38}, {0.075, 3680.8}, {0.085, 3779.16}, {0.095, 3861.68}, {0.105, 3931.94}, {0.115, 3992.52}, {0.125, 4045.3}, {0.135, 4091.71}, {0.145, 4132.84}, {0.155, 4169.56}, {0.165, 4202.53}, {0.175, 4232.31}, {0.185, 4259.34}, {0.195, 4283.98}, {0.205, 4306.53}, {0.215, 4327.26}, {0.225, 4346.37}, {0.235, 4364.04}, {0.245, 4380.44}, {0.255, 4395.69}, {0.265, 4409.91}, {0.275, 4423.2}, {0.285, 4435.65}, {0.295, 4447.33}, {0.305, 4458.32}, {0.315, 4468.67}, {0.325, 4478.45}, {0.335, 4487.68}, {0.345, 4496.43}, {0.355, 4504.72}, {0.365, 4512.59}, {0.375, 4520.07}, {0.385, 4527.19}, {0.395, 4533.97}, {0.405, 4540.44}, {0.415, 4546.63}, {0.425, 4552.53}, {0.435, 4558.19}, {0.445, 4563.6}, {0.455, 4568.79}, {0.465, 4573.77}, {0.475, 4578.56}, {0.485, 4583.15}, {0.495, 4587.57}}

Reference:

FindFit with a sophisticated function (integral)

FindFit working with Integrate but not NIntegrate - I am using "?NumericQ"

Where can I find examples of good Mathematica programming practice?

$\endgroup$
3
$\begingroup$

I found an error with the syntax in the model where a curly brace was used rather than parenthesis

Exp[-{Log[w/δ]}^2/(2 σ^2)] -> Exp[-(Log[w/δ])^2/(2 σ^2)]

With this change model is defined as:

model[Ms_?NumericQ, δ_?NumericQ, σ_?NumericQ, x_?NumericQ] := 
 NIntegrate[
  Ms (Coth[(Ms x (1/6 Pi w^3))/(k T)] - ((k T)/(Ms x (1/6 Pi w^3))))
   1/(Sqrt[2 Pi] w σ) Exp[-(Log[w/δ])^2/(2 σ^2)],
  {w, 10*^-10, 10^-7}]

The numerical optimization functions need a reasonable starting point in order to work properly. I recommend using Manipulate to try to get starting values that are in the ball park.

Manipulate[
 Module[
  {
   dataRecon = Transpose@Join[{data[[All, 1]],
       Map[model[Ms, δ, σ, #] &, data[[All, 1]]]}]
   },

  Show[
   ListPlot[data, PlotStyle -> Black],
   ListPlot[dataRecon, PlotStyle -> Red],
   Graphics[{
     Red,
     Line[dataRecon]
     }],
   PlotRange -> All
   ]
  ],

 {{Ms, 8000.}, 5000, 15000},
 {{δ, 6*10^-8}, 5*10^-9, 5*10^-7},
 {{σ, 1.}, 0.01, 10}
 ]

Mathematica graphics

I tried FindFit but after 30 minutes aborted. Sometimes FindMinimum works better in fussy cases.

One has to define a function as input to FindMinimum. For fitting problems use the sum of the difference between the measured and reconstructed data squared.

I call this the objectiveFuncion. Also when attempting to fit parameters that are orders of magnitude different from one another, I normally scale the parameters so that they are in the same range.

objectiveFunction[data_, MsScaled_?NumericQ,
  δScaled_?NumericQ, σ_?NumericQ] := Module[
  {
   Ms = 1000*MsScaled,
   δ = δScaled*10^-8,
   dataRecon,
   residual
   },

  dataRecon = Map[model[Ms, δ, σ, #] &, data[[All, 1]]];
  residual = dataRecon - data[[All, 2]];

  (* Sum of the residuals squared *)

  Total@Map[#^2 &, residual]
  ]

Now run FindMinimum to estimate the parameters as follows:

FindMinimum[
 {
  objectiveFunction[data, MsScaled, δScaled, σ],
  5 < MsScaled < 15,
  0.5 < δScaled < 50,
  0.1 < σ < 2
  },
 {
  {MsScaled, 8},
  {δScaled, 6},
  {σ, 1.0}
  }
 ]

{53784.1, {MsScaled -> 5.34602, δScaled -> 5.45066, σ -> 0.465785}}

In order to validate it plot the measured and reconstructed data

Module[
 {
  Ms = 5346.0,
  δ = 5.45*10^-8,
  σ = 0.4658,
  dataRecon
  },

 dataRecon = Transpose@Join[{data[[All, 1]],
     Map[model[Ms, δ, σ, #] &, data[[All, 1]]]}];

 Show[
  ListPlot[data, PlotStyle -> Black],
  ListPlot[dataRecon, PlotStyle -> Red],
  Graphics[{
    Red,
    Line[dataRecon]
    }],
  PlotRange -> All
  ]
 ]

Mathematica graphics

Compared to the Manipulate guess this fits the data better.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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