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
    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
    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
    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
    Oct 24, 2023 at 21:31
  • $\begingroup$ No. Square is differentiable, 'RealAbs' isn't. $\endgroup$
    – user64494
    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
    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$ Oct 27, 2023 at 13:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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