Search code examples
wolfram-mathematicaimplicit-conversiondata-fittingminimizemodel-fitting

Fitting data for an implicit function with 2 parameters

I have an implicit function of independent variable sB and dependent variable H with two parameters kV and km, describing an chemical experiment. I tried to fit the experimentally obtained logarithmic points to the model provided by article using Wolfram Mathematica.
I turned the coordinates of the points into the values of the concentrations of sB(left) and H(right) for further use.

dotsA = {{0.5076923076923077, 7.751879699248121},{1.0153846153846153, 7.604511278195489},
{1.7769230769230773, 7.4},{2.576923076923077,7.209022556390978},
{3.2692307692307696, 7.070676691729323},{4.530769230769231, 7.003007518796992}};

dotsconcA = 10^-dotsA
{{0.310676, 1.7706*10^-8},{0.0965196, 2.48593*10^-8},
{0.0167139, 3.98107*10^-8},{0.00264897, 6.17984*10^-8},
{0.000537984, 8.49813*10^-8},{0.0000294599, 9.93099*10^-8}}

Constants and Model

kaA = 4.3*10^-7; kaB = 5.6*10^-10; kaW = 6.2*10^-8;
nA = 1; nB = 2; kH = 1; kW = 1; HB = 10^-7;
cBWA = 1*10^-3;
kV111 = 0.41*10^-3;
km111 = 13.1*10^-3;

kH*(HB - H) + kW*cBWA*(1/(1 + kaW/HB) - 1/(1 + kaW/H)) + 
(kV*1/2*(sB - kV - km + Sqrt[(sB - kV - km)^2 + 4*km*sB]))/(km + 1/2*(sB - kV - km + Sqrt[(sB - kV - km)^2 + 4*km*sB]))*
(nA/(1 + H/kaA) - nB/(1 + kaB/H)) == 0

Where kV111 and km111 are parameters that were obtained using the Least Squares Method in the article, but the article does not explain this in sufficient details, so I'm trying to find them on my own.
I found roots of function and made graph using parameters kv111 and km111

f1[H_, sB_] := kH*(HB - H) + kW*cBWA*(1/(1 + kaW/HB) - 1/(1 + kaW/H)) + 
(kv111*1/2*(sB - kv111 - km111 + Sqrt[(sB - kv111 - km111)^2 + 4*km111*sB]))/(km111 + 1/2*(sB - kv111 - km111 + Sqrt[(sB - kv111 - km111)^2 + 4*km111*sB]))*
(nA/(1 + H/kaA) - nB/(1 + kaB/H));

lst1 = Table[{-Log[10, sB], -Log[10, H /. FindRoot[f1[H, sB] == 0, {H, 10^-9}]]}, {sB, 10^-5, 1, 0.0001}]

and compared it to experimentally obtained points.

https://i.sstatic.net/3EHvg.png This is a link to my plot from my question on Mathematica stack exchange

This plot looks exactly the same as one in the article.
I tried to fit data to model using the FindFit but it gives errors and values that not even close to found ones. None of the parameters and variables can be negative.

FindFit[dotsconcA, f[sB, H, kV, km], {kV, km}, {sB, H}]
FindFit:Number of coordinates (1) is not equal to the number of variables (2).

fdot1 = f1[dotsconcA[[All, 1]], dotsconcA[[All, 2]], kV111, km111];
dataA1 = Transpose[Append[Transpose[dotsconcA], fdot1]];
FindFit[dataA1, f1[sB, H, kV, km], {kV, km}, {sB, H}]
FindFit:The function value {0.0000237276+0.I, -0.0000526427+0.I, -0.00132824+0.00202876I, 0.0122705+0.I, 0.0152136+0.I, 0.0161337+0.I} is not a list of real numbers with dimensions {6} at {kV,km} = {0.000352367,-0.0140961}.
{kV -> 0.000352367, km -> -0.0140961}

How can I find the kV and km parameters in Wolfram Math by myself?

UPD: So I tried the LSM that was given in the part 2 of article. This method implies that the sum of the squares of the differences between the values of H experimental and H calculated analytically should be minimal, i.e. as close as possible to 0.

f2[sB_,H_]=kH*(HB - H) + kW*cBWA*(1/(1 + kaW/HB) - 1/(1 + kaW/H)) + 
(kv111*1/2*(sB - kv111 - km111 + Sqrt[(sB - kv111 - km111)^2 +4*km111*sB]))
/(km111 + 1/2*(sB - kv111 - km111 + Sqrt[(sB - kv111 - km111)^2 +4*km111*sB]))*
(nA/(1 + H/kaA) - nB/(1 + kaB/H))

Sum[Log[10,
(H/.FindRoot[f2[dotsconcA[[i,1]],H]==0, {H,10^14}])/dotsconcA[[i,2]]]^2, {i,1,6}]
Clear[i]
0.0117239

Here I converted PH values to concentration and brought them under one logarithm. With the same sB value, we have two different Hs.
We must find this minimum by manipulating kV and km values. I tried to use Manipulate to manually find kV and km and it worked

Manipulate[Sum[Log[10,
(H/.FindRoot[f3[dotsconcA[[i,1]],H,kV1, km1]==0, {H,10^-14}])/dotsconcA[[i,2]]]^2,
{i,1,6}], {kV1, 0.41*10^-4, 0.41*10^-2.5, 0.41*10^-4},
{km1, 13.1*10^-4, 13.1*10^-2.5, 13.1*10^-4}]

But this

  1. Takes a lot of time,
  2. gives not accurate results,
  3. requires approximate ranges for manipulated parameters, which we don't know initially.

Is there a better way to find minimum of this Sum? I tried several "Minimization" functions but they're not working

Solution

  • So we found our fit and now I will describe how if someone will face the same problem. First of all, we simplified our model to a polynomial of the 4th degree.

    f1[sB_,H_,kV_,km_]:=(H+kaA) (H+kaB) (H+kaW) (HB+kaW) ((-H+HB) kH+(cBWA (-H+HB) kaW kW)/((H+kaW) (HB+kaW))+(kV (-kaA (H+kaB) nA+H (H+kaA) nB) (km+kV-sB-Sqrt[(km+kV-sB)^2+4 km sB]))/((H+kaA) (H+kaB) (km-kV+sB+Sqrt[(km+kV-sB)^2+4 km sB])))

    Secondly, we used Solve to find a solution for H and defined a function of sB and 2 parameters km and kV.

    poly4[sB_,kV_,km_]:=H/.Solve[f1[sB,H,kV,km]==0&&0<=sB<=1,H,PositiveReals,InverseFunctions->True][[1]]

    And finally the fitting

    FindFit[dotsconcA,poly4[sB,kV,km], {{km,0.00001}, {kV,0.00001}}, sB]
km->0.006484, kV->0.0003568

    Unfortunately, these parameters do not match those in the article, but when comparing the plots it looks great. Also, the fitting process takes about an hour, but that's a problem for the next question.