1
$\begingroup$

exp - an experimental data, where expT - temperature, expChi - magnetic susceptibility

I would like to describe this data using a function chi[a, b, g, be, TC, c, o, T], where a, b, g, be, TC, c, o are parameters.

I try to find the best parameters using minimization of function chifun[a, b, g, be, TC, c, o]

The possible parameter ranges: a>0, b>0, c>0, d>0, C>0, TC>0, o<0 (TC ≈ 315, o ≈ -3, typical values from the literature: a=1.166; 1.228, b=0,0086; 0,0077, d=0.4)

Minimization occurs, but the theoretical curve does not resemble the experiment. What's wrong in the code?

ClearAll["Global`*"]
H = 1000;
(*c,o,a,b,g,be,TC are parameters*)
chilt[c_, o_, T_] := c/(T - o);
(*chilt[0.4,-3,10]*)
chiht[a_, b_, g_, be_, TC_, T_] := 
 Module[{}, 
  eq = a*{T - TC}*chih^{1/g} + b*chih^{1/g + 1/be}*H^{1/be} - 1 == 0;
  sol = FindRoot[eq, {chih, 1}];
  Return[sol]]
(*chiht[0.79,0.00893,1.38,0.428,320,10]*)
chi[a_, b_, g_, be_, TC_, c_, o_, T_] := 
 chilt[c, o, T] + chiht[a, b, g, be, TC, T][[1, 2]] // Chop
(*chi[0.79,0.00893,1.38,0.428,320,0.4,-3,10]*)
(*experimental data*)
exp = {{10.0052`, 0.0879716`}, {11.0433`, 
    0.08638660000000001`}, {12.359`, 0.0848015`}, {13.7127`, 
    0.0832164`}, {14.9939`, 0.0820276`}, {16.3367`, 
    0.0808388`}, {17.6282`, 0.0800463`}, {18.9464`, 
    0.0788575`}, {18.9934`, 0.0788575`}, {18.9969`, 
    0.0788575`}, {18.9962`, 0.0788575`}, {19.4221`, 
    0.07846120000000001`}, {20.8034`, 
    0.07766870000000001`}, {22.1306`, 0.0764799`}, {23.4352`, 
    0.0756873`}, {24.7503`, 0.075291`}, {26.0295`, 
    0.07449850000000001`}, {27.3692`, 
    0.07370600000000001`}, {30.0264`, 0.0729134`}, {31.3262`, 
    0.0725172`}, {32.6313`, 0.0721209`}, {33.9404`, 
    0.0717246`}, {35.2494`, 0.0713284`}, {36.551`, 
    0.07093210000000001`}, {37.8445`, 0.0705358`}, {40.5813`, 
    0.06974330000000001`}, {41.7696`, 
    0.06934699999999999`}, {43.0624`, 
    0.06934699999999999`}, {44.3609`, 0.0689507`}, {45.6536`, 
    0.0685545`}, {46.9616`, 0.0685545`}, {48.254`, 
    0.06815819999999999`}, {49.5428`, 0.0677619`}, {52.1298`, 
    0.0673657`}, {53.4147`, 0.0673657`}, {54.7091`, 
    0.0669694`}, {56.0001`, 0.0669694`}, {57.2883`, 
    0.0665731`}, {58.5865`, 0.0665731`}, {59.8911`, 
    0.06617690000000001`}, {62.4856`, 
    0.06617690000000001`}, {63.7665`, 
    0.06578060000000001`}, {65.0689`, 
    0.06578060000000001`}, {66.3588`, 
    0.06538429999999999`}, {67.6581`, 
    0.06538429999999999`}, {68.9546`, 
    0.06538429999999999`}, {70.2572`, 
    0.06498810000000001`}, {72.8482`, 
    0.06498810000000001`}, {74.1637`, 
    0.06498810000000001`}, {75.4546`, 0.0645918`}, {76.7495`, 
    0.0645918`}, {78.0347`, 0.0645918`}, {79.3376`, 
    0.0641955`}, {80.6437`, 0.0641955`}, {83.5861`, 
    0.0641955`}, {84.8957`, 0.0637993`}, {86.2123`, 
    0.0637993`}, {87.5317`, 0.0637993`}, {88.8468`, 
    0.0637993`}, {90.1674`, 0.0637993`}, {91.4775`, 
    0.063403`}, {94.2769`, 0.063403`}, {95.5863`, 
    0.063403`}, {96.9007`, 0.063403`}, {98.2013`, 
    0.063403`}, {99.523`, 0.0630067`}, {100.839`, 
    0.0630067`}, {102.147`, 0.0630067`}, {104.812`, 
    0.0630067`}, {106.131`, 0.0630067`}, {107.448`, 
    0.0626104`}, {108.789`, 0.0626104`}, {110.099`, 
    0.0626104`}, {111.414`, 0.0626104`}, {112.756`, 
    0.0626104`}, {115.404`, 0.0626104`}, {116.737`, 
    0.0626104`}, {118.059`, 0.0626104`}, {119.376`, 
    0.0622142`}, {120.7`, 0.0622142`}, {122.025`, 
    0.0622142`}, {123.367`, 0.0622142`}, {126.006`, 
    0.0622142`}, {127.352`, 0.0622142`}, {128.674`, 
    0.0622142`}, {129.999`, 0.0622142`}, {131.333`, 
    0.0622142`}, {132.652`, 0.0618179`}, {133.987`, 
    0.0618179`}, {136.652`, 0.0618179`}, {137.987`, 
    0.0618179`}, {139.324`, 0.0618179`}, {140.645`, 
    0.0618179`}, {141.989`, 0.0618179`}, {143.316`, 
    0.0618179`}, {144.646`, 0.0618179`}, {147.323`, 
    0.0618179`}, {148.651`, 0.0614216`}, {149.978`, 
    0.0614216`}, {151.303`, 0.0614216`}, {152.637`, 
    0.0614216`}, {153.949`, 0.0614216`}, {155.268`, 
    0.0614216`}, {157.929`, 0.0614216`}, {159.252`, 
    0.0614216`}, {160.581`, 0.0610254`}, {161.889`, 
    0.0610254`}, {163.219`, 0.0610254`}, {164.546`, 
    0.0610254`}, {165.873`, 0.0610254`}, {168.524`, 
    0.0610254`}, {169.852`, 0.060629100000000005`}, {171.173`, 
    0.060629100000000005`}, {172.494`, 
    0.060629100000000005`}, {173.83`, 
    0.060629100000000005`}, {175.16`, 
    0.060629100000000005`}, {176.473`, 
    0.060629100000000005`}, {179.122`, 
    0.060232799999999996`}, {180.459`, 
    0.060232799999999996`}, {181.78`, 
    0.060232799999999996`}, {183.106`, 0.0598366`}, {184.439`, 
    0.0598366`}, {185.757`, 0.0598366`}, {187.079`, 
    0.0598366`}, {189.73`, 0.0594403`}, {191.06`, 
    0.0594403`}, {192.386`, 0.059044`}, {193.717`, 
    0.059044`}, {195.037`, 0.059044`}, {196.372`, 
    0.0586478`}, {197.693`, 0.0586478`}, {200.361`, 
    0.058251500000000005`}, {201.689`, 0.0578552`}, {203.008`, 
    0.0578552`}, {204.335`, 0.057459`}, {205.651`, 
    0.057459`}, {206.996`, 0.0570627`}, {208.324`, 
    0.0570627`}, {210.977`, 0.056666400000000006`}, {212.292`, 
    0.056270200000000006`}, {213.602`, 0.0558739`}, {214.947`, 
    0.0558739`}, {216.259`, 0.0554776`}, {217.583`, 
    0.0550813`}, {218.895`, 0.0550813`}, {221.56`, 
    0.054288800000000005`}, {222.865`, 
    0.053892499999999996`}, {224.195`, 0.0534963`}, {225.505`, 
    0.0531`}, {226.824`, 0.0531`}, {228.192`, 0.0527037`}, {229.508`, 
    0.0523075`}, {232.191`, 0.051514899999999995`}, {233.507`, 
    0.051118699999999996`}, {234.836`, 0.0507224`}, {236.161`, 
    0.0503261`}, {237.479`, 0.049929900000000006`}, {238.807`, 
    0.049533600000000004`}, {240.13`, 0.048741`}, {242.774`, 
    0.047948500000000005`}, {244.096`, 0.0475522`}, {245.417`, 
    0.047155999999999997`}, {246.74`, 0.0467597`}, {248.057`, 
    0.0463634`}, {249.382`, 0.045570900000000004`}, {250.457`, 
    0.0462508`}, {251.117`, 0.0460365`}, {252.897`, 
    0.0454166`}, {253.112`, 0.0453033`}, {254.115`, 
    0.0449356`}, {255.112`, 0.044570700000000005`}, {256.107`, 
    0.044217900000000004`}, {257.936`, 0.043554`}, {259.077`, 
    0.0431208`}, {260.395`, 0.042633000000000004`}, {260.555`, 
    0.042568800000000004`}, {261.881`, 0.0420749`}, {263.515`, 
    0.0414803`}, {264.389`, 0.041108`}, {264.559`, 
    0.041047`}, {264.724`, 0.0409939`}, {264.891`, 
    0.0409305`}, {265.058`, 0.0408548`}, {265.219`, 
    0.040807`}, {265.384`, 0.0407401`}, {265.549`, 
    0.040679900000000005`}, {265.707`, 
    0.040611600000000005`}, {265.872`, 0.0405381`}, {266.041`, 
    0.040491900000000004`}, {266.218`, 0.0404131`}, {266.379`, 
    0.040365000000000005`}, {266.532`, 0.0402826`}, {266.705`, 
    0.0402386`}, {266.876`, 0.0401762`}, {267.038`, 
    0.040110900000000005`}, {267.201`, 0.0400532`}, {267.373`, 
    0.039997200000000004`}, {267.543`, 0.0399247`}, {267.71`, 
    0.0398555`}, {267.874`, 0.0397915`}, {268.035`, 
    0.0397336`}, {268.194`, 0.039681100000000004`}, {268.367`, 
    0.039605600000000005`}, {268.542`, 
    0.039538699999999996`}, {268.705`, 
    0.039476199999999996`}, {268.866`, 0.0393997`}, {269.029`, 
    0.039342300000000004`}, {269.201`, 0.0392861`}, {269.374`, 
    0.03922`}, {269.539`, 0.039153799999999996`}, {269.698`, 
    0.0390823`}, {269.87`, 0.0390231`}, {270.036`, 
    0.0389709`}, {270.194`, 0.038892499999999997`}, {270.371`, 
    0.038839399999999996`}, {270.536`, 0.0387638`}, {270.697`, 
    0.0387051`}, {270.868`, 0.038642499999999996`}, {271.038`, 
    0.038574`}, {271.208`, 0.0385115`}, {271.367`, 
    0.0384455`}, {271.54`, 0.0383782`}, {271.711`, 
    0.038315600000000005`}, {271.868`, 0.0382524`}, {272.03`, 
    0.038191699999999995`}, {272.2`, 0.0381262`}, {272.361`, 
    0.0380706`}, {272.525`, 0.0380021`}, {274.156`, 
    0.0373821`}, {274.229`, 0.037318800000000006`}, {274.362`, 
    0.037270000000000005`}, {274.524`, 
    0.037204799999999996`}, {274.705`, 0.0371411`}, {274.877`, 
    0.0370827`}, {275.029`, 0.0370104`}, {275.196`, 
    0.036960700000000006`}, {275.37`, 
    0.036885100000000004`}, {275.526`, 0.0368172`}, {275.683`, 
    0.0367629`}, {275.855`, 0.0366982`}, {276.025`, 
    0.036634599999999996`}, {276.185`, 0.0365695`}, {276.345`, 
    0.036506300000000005`}, {276.506`, 0.0364347`}, {276.68`, 
    0.036378900000000006`}, {276.845`, 0.036315`}, {277.012`, 
    0.036263199999999995`}, {277.184`, 
    0.036168900000000004`}, {277.35`, 0.0361295`}, {277.517`, 
    0.0360462`}, {277.676`, 0.0359975`}, {277.839`, 
    0.0359304`}, {278.009`, 0.0358579`}, {278.18`, 
    0.035811100000000005`}, {278.341`, 0.0357453`}, {278.5`, 
    0.0356707`}, {278.67`, 0.0356269`}, {278.836`, 
    0.0355542`}, {279.005`, 0.035486`}, {279.176`, 
    0.0354266`}, {279.333`, 0.0353622`}, {279.502`, 
    0.0352938`}, {279.671`, 0.0352304`}, {279.833`, 
    0.0351688`}, {279.997`, 0.0351014`}, {280.169`, 
    0.035043700000000004`}, {280.337`, 0.0349787`}, {280.498`, 
    0.0349181`}, {280.67`, 0.0348453`}, {280.837`, 
    0.0347763`}, {280.995`, 0.0347195`}, {281.164`, 
    0.0346492`}, {281.332`, 0.034588900000000006`}, {281.5`, 
    0.034523000000000005`}, {281.664`, 0.0344628`}, {281.83`, 
    0.0343953`}, {282.004`, 0.0343438`}, {282.164`, 
    0.034283499999999995`}, {282.324`, 0.034214`}, {282.493`, 
    0.0341503`}, {282.67`, 0.034087000000000006`}, {282.833`, 
    0.0340348`}, {282.996`, 0.0339549`}, {283.164`, 
    0.0339051`}, {284.809`, 0.033290999999999994`}, {284.864`, 
    0.0332491`}, {284.993`, 0.0331839`}, {285.16`, 
    0.033110999999999995`}, {285.323`, 
    0.033054099999999996`}, {285.487`, 0.0329876`}, {285.655`, 
    0.0329227`}, {285.821`, 0.0328639`}, {285.991`, 
    0.0328064`}, {286.164`, 0.032734200000000005`}, {286.329`, 
    0.0326806`}, {286.494`, 0.032606`}, {286.66`, 
    0.032562`}, {286.824`, 0.0324928`}, {286.985`, 
    0.0324283`}, {287.149`, 0.0323634`}, {287.321`, 
    0.032287300000000005`}, {287.487`, 0.0322322`}, {287.644`, 
    0.0321765`}, {287.817`, 0.032120800000000005`}, {287.988`, 
    0.0320552`}, {288.146`, 0.0319866`}, {288.306`, 
    0.0319345`}, {288.476`, 0.031862100000000004`}, {288.646`, 
    0.0318045`}, {288.815`, 0.0317536`}, {288.985`, 
    0.0316971`}, {289.14`, 0.0316201`}, {289.307`, 
    0.031550999999999996`}, {289.481`, 0.0314893`}, {289.645`, 
    0.031444`}, {289.81`, 0.0313815`}, {289.975`, 
    0.031311`}, {290.141`, 0.031252`}, {290.306`, 
    0.0311829`}, {290.466`, 0.0311326`}, {290.636`, 
    0.0310702`}, {290.807`, 0.031011`}, {290.967`, 
    0.030944600000000003`}, {291.139`, 0.0308823`}, {291.304`, 
    0.0308283`}, {291.47`, 0.0307563`}, {291.646`, 
    0.030706`}, {291.809`, 0.030651599999999998`}, {291.974`, 
    0.030578`}, {292.141`, 0.0305155`}, {292.306`, 
    0.0304571`}, {292.47`, 0.030398599999999998`}, {292.634`, 
    0.0303283`}, {292.805`, 0.0302715`}, {292.971`, 
    0.030212`}, {293.127`, 0.0301523`}, {293.299`, 
    0.0300891`}, {293.475`, 0.0300263`}, {293.634`, 
    0.0299634`}, {293.794`, 0.0299016`}, {295.41`, 
    0.029346900000000002`}, {295.467`, 
    0.029285900000000004`}, {295.612`, 0.0292319`}, {295.787`, 
    0.029188600000000002`}, {295.954`, 0.0291198`}, {296.116`, 
    0.0290857`}, {296.277`, 0.029017800000000003`}, {296.445`, 
    0.0289598`}, {296.617`, 0.0289062`}, {296.78`, 
    0.0288415`}, {296.946`, 0.0287719`}, {297.111`, 
    0.028703700000000002`}, {297.277`, 0.0286442`}, {297.444`, 
    0.0285661`}, {297.603`, 0.0285065`}, {297.772`, 
    0.028437999999999998`}, {297.945`, 0.0283781`}, {298.11`, 
    0.0283183`}, {298.27`, 0.0282553`}, {298.436`, 
    0.0281991`}, {298.608`, 0.0281362`}, {298.775`, 
    0.028089`}, {298.934`, 0.0280215`}, {299.089`, 
    0.027953600000000002`}, {299.26`, 
    0.027889800000000003`}, {299.431`, 
    0.027837900000000002`}, {299.596`, 0.0277736`}, {299.77`, 
    0.0277224`}, {299.939`, 0.0276707`}, {300.1`, 
    0.027607500000000004`}, {300.268`, 0.0275584`}, {300.436`, 
    0.0274834`}, {300.595`, 0.027445`}, {300.761`, 
    0.027380500000000002`}, {300.928`, 0.0273191`}, {301.085`, 
    0.0272599`}, {301.248`, 0.0272049`}, {301.416`, 
    0.0271372`}, {301.582`, 0.0270904`}, {301.742`, 
    0.0270381`}, {301.902`, 0.026959499999999997`}, {302.077`, 
    0.0269182`}, {302.254`, 0.026863900000000003`}, {302.409`, 
    0.0267965`}, {302.573`, 0.026751`}, {302.744`, 
    0.0266902`}, {302.904`, 0.0266289`}, {303.075`, 
    0.0265635`}, {303.24`, 0.026511600000000003`}, {303.398`, 
    0.0264599`}, {303.571`, 0.0263952`}, {303.75`, 
    0.0263447`}, {303.915`, 0.0262879`}, {304.07`, 
    0.0262416`}, {304.236`, 0.0261701`}, {304.396`, 
    0.026108799999999998`}};
expT = exp /. {f_, l_} :> f;
expChi = exp /. {f_, l_} :> l;
chifun[a_, b_, g_, be_, TC_, c_, o_] := 
 Sum[Abs[chi[a, b, g, be, TC, c, o, expT[[i]]] - expChi[[i]]], {i, 
   Length[exp]}]
ip = {0.79, 0.00893, 1.38, 0.428, 320, 0.4, -3};
param = NMinimize[{chifun[a, b, g, be, TC, c, 
    o], {(0 < a < 3) && (0 < b < 3) && (0 < g < 3) && (0 < be < 
       3) && (280 < TC < 350) && (0 < c < 5) && (-100 < o < 0)}}, {a, 
   b, g, be, TC, c, o}, 
  Method -> {"Automatic", "InitialPoints" -> ip}, 
  EvaluationMonitor :> {Print["a=", a, "   b=", b, "   g=", g, 
     "   be=", be]}]
{2.7194945088238707`, {a -> 2.2630679452300386`, 
  b -> 2.9906490183683014`, g -> 1.378825817551753`, 
  be -> 2.3093381293268127`, TC -> 302.85700145284545`, c -> 5.`, 
  o -> -67.17312860068925`}}

gexp = ListPlot[exp]

g2 = Show[
  Plot[chi[param[[2, 1, 2]], param[[2, 2, 2]], param[[2, 3, 2]], 
    param[[2, 4, 2]], param[[2, 5, 2]], param[[2, 6, 2]], 
    param[[2, 7, 2]], T], {T, 0, 300}], gexp, 
  PlotRange -> {{0, 300}, {0, 0.08}}]
$\endgroup$
7
  • $\begingroup$ what is chih in the definition of chiht? $\endgroup$
    – ydd
    Commented Oct 24, 2023 at 18:08
  • $\begingroup$ @ydd, thanks for the comment. chih is a magnetic susceptibility at high temperatures. a*{T - TC}*chih^{1/g} + b*chih^{1/g + 1/be}*H^{1/be} - 1=0 - the equation is implicit, so I use Module to find chih. $\endgroup$
    – Mam Mam
    Commented Oct 24, 2023 at 18:12
  • $\begingroup$ Did you try to minimize the sum of squares (a*{T - TC}*chih^{1/g} + b*chih^{1/g + 1/be}*H^{1/be} - 1)^2 over the exp? $\endgroup$
    – user64494
    Commented Oct 24, 2023 at 19:45
  • $\begingroup$ @user64494, thanks for the comment. Here I minimize the sum of the absolute values of the difference between the theoretical and experimental values chifun[a_, b_, g_, be_, TC_, c_, o_] := Sum[Abs[chi[a, b, g, be, TC, c, o, expT[[i]]] - expChi[[i]]], {i, Length[exp]}]. Isn't this the same approach as their square? $\endgroup$
    – Mam Mam
    Commented Oct 24, 2023 at 21:31
  • $\begingroup$ No. Square is differentiable, 'RealAbs' isn't. $\endgroup$
    – user64494
    Commented Oct 24, 2023 at 21:35

1 Answer 1

1
$\begingroup$

We can turn an implicit function into an explicit one using Nest as follows

eq = a*(T - TC)*x^(1/g) + b*x^(1/g + 1/be)*H^(1/be) - 1; ini0 = {0.79,
   0.00893, 1.38, 0.428, 320, 0.4, -3}; data = {304.396`, 
  0.026108799999999998`}; ini = {a -> 0.846702954609469`, 
  b -> 0.006132039075224503`, g -> 1.4706605549510186`, 
  be -> 0.40313022491307543`, TC -> 319.9976258605692`, 
  c -> 0.07647955132052105`, o -> -3.006399549360516`, H -> 10^3, 
  T -> data[[1]]};

(*x=((1-a*(T-TC)*x^(1/g))/(b H^(1/be)))^(g be/(g+be))*)
f = Nest[((1 - 
         a*(T - TC)*#^(1/g))/(b H^(1/be)))^(g be/(g + be)) &, .087, 
   7];

Let check that f works same as FindRoot

FindRoot[eq == 0 /. ini, {x, data[[2]]}]

Out[]= {x -> 0.0283069}

f /. ini

Out[]= 0.0283069

Using f with FindMinimum we can optimize model. Note that experimental data exp is too long for optimization, therefore we use part of exp in a form (please, take exp from the main post)

data = Join[Take[Drop[exp, -200], {1, -1, 10}], 
   Take[Take[exp, -200], {1, -1, 50}]];
expT = data[[All, 1]];
expChi = data[[All, 2]];
fun = ((f + c/(T - o)) /. T -> expT) - expChi;
var = {a, b, g, be, 
  TC, c, o}; ip = Take[ini, Length[var]][[All, 2]] 

Solution

con = {(10^-8 < a < 3) && (10^-8 < b < 3) && (0.1 < g < 
     3) && (0.1 < be < 3) && (280 < TC < 350) && (10^-8 < c < 
     5) && (-100 < o < -10^-8)}; H = 1000; param = 
 FindMinimum[{fun . fun, con}, 
  Table[{var[[i]], ip[[i]]}, {i, Length[var]}], MaxIterations -> 1000]

(*Out[]= {0.0000991514, {a -> 2.67021, b -> 7.59139*10^-6, 
  g -> 1.19684, be -> 0.229313, TC -> 293.396, c -> 0.190081, 
  o -> -0.0411701}}*)

Visualization

Show[ListLinePlot[
  Transpose[{exp[[All, 
     1]], ((f + c/(T - o)) /. T -> exp[[All, 1]]) /. param[[2]]}], 
  PlotRange -> All], ListPlot[exp, PlotStyle -> Red]]

Figure 1

Note that the coincidence is not the best. The model can be improved.

$\endgroup$
2
  • $\begingroup$ Thank you very much for the detailed explanation! $\endgroup$
    – Mam Mam
    Commented Oct 27, 2023 at 8:44
  • 1
    $\begingroup$ You are welcome. Please note, this model is not perfect fit for experimental data exp. The data was likely combined from several sources. $\endgroup$ Commented Oct 27, 2023 at 13:35

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