[?] Curve fitting with several implicit functions?

Questions
, ,
1

Hello,
I am looking for a way to do a curve fitting by using several simultaneous implicit equations (and I couldn’t make this work with FindFit). Attached is an notebook that I prepared for generating a plotting the data of the following coupled equations and NSolve

NSolve[{10^pK1 == HD/(H*D), 10^pK2 == HG/(H*G), Ha == HD + HG + H, 
  Ga == HG + G, Da == HD + D}, {HD, HG, H, G, D}]

It followed making solution matrix for discrete values of Ga and picking out the correct parts of the solution of HD, HG, H, G and D that are real and fulfill some conditions given pK1, pK2, Ha and Da as fixed values, and Ga as the independent variable,

I have now the problem that I would like to do the opposite: Given a data list of the form {Ga, HD}, I would like fit the parameters pK1 and Ha, which best fulfill the equations 10^pK1 == HD/(HD), 10^pK2 == HG/(HG), Ha == HD + HG + H,
Ga == HG + G, Da == HD + D as above.

For instance, the 20-point data list (Ga, HD) below was created with Ha = 1, Da =1, pk1 = 2 and pK2 =3.

{{1.*10^-8, 0.904875}, {0.05, 0.878427}, {0.1, 0.846593}, {0.15, 
0.81062}, {0.2, 0.7719}, {0.25, 0.731591}, {0.3, 0.690551}, {0.35, 
0.649405}, {0.4, 0.608618}, {0.45, 0.568561}, {0.5, 
0.529541}, {0.55, 0.491828}, {0.6, 0.455658}, {0.65, 0.42124}, {0.7,
0.38875}, {0.75, 0.358325}, {0.8, 0.330057}, {0.85, 0.30399}, {0.9,
0.280116}, {0.95, 0.25838}}

If it was a set of differential equations, I would have used ParametricNDSolveValue to create a model and FindFit evaluating the model. But I have no good idea how to efficiently code for the set of "normal" equations shown above that give, when solved with Solve, massive long solutions. FindRoot appeared "dangerous" as it often found the "wrong" branch of the solution, depending on the initialization value.

Your help would be very much appreciated how I can fit the parameters for the set of equations shown above.
If you need, please feel free to create more list data with the interactive notebook that is attached. (The data list is called “HDpoints”).
best regards,
Frank

2

Your code does not evaluate, because you do not provide the starting values of Ha etc. I am not sure what you want, but perhaps you could try a Manipulate:

hdPts[Ha_, Da_, Gmax_, pK1_, pK2_, displaypoints_] := 
  Module[{HD, HG, H, G, D},
   HDpoints = {}; HGpoints = {}; Hpoints = {}; Gpoints = {}; 
   Dpoints = {}; 
   Quiet[Flatten[
     Flatten[Delete[
       Do[{solut = 
          N[NSolve[{10^pK1 == HD/(H*D), 10^pK2 == HG/(H*G), 
             Ha == HD + HG + H, Ga == HG + G, Da == HD + D}, {HD, HG, 
             H, G, D}], 30], HDfirst = HD /. Part[solut, 1, 4], 
         HDsecond = HD /. Part[solut, 2, 4], 
         HDthird = HD /. Part[solut, 3, 4], 
         HGfirst = HG /. Part[solut, 1, 5], 
         HGsecond = HG /. Part[solut, 2, 5], 
         HGthird = HG /. Part[solut, 3, 5], 
         Hfirst = H /. Part[solut, 1, 1], 
         Hsecond = H /. Part[solut, 2, 1], 
         Hthird = H /. Part[solut, 3, 1], 
         Gfirst = G /. Part[solut, 1, 2], 
         Gsecond = G /. Part[solut, 2, 2], 
         Gthird = G /. Part[solut, 3, 2], 
         Dfirst = D /. Part[solut, 1, 3], 
         Dsecond = D /. Part[solut, 2, 3], 
         Dthird = D /. Part[solut, 3, 3],
         If[
          0 <= HDfirst <= Da && 0 <= HGfirst <= Ga && 
           0 <= Dfirst <= Da && 0 <= Gfirst <= Ga && 
           0 <= Hfirst <= Ha && HGfirst <= Ha && HDfirst <= Ha && 
           Hfirst + HGfirst + HDfirst <= Ha && 
           Gfirst + HGfirst <= Ga && 
           Dfirst + HDfirst <= Da, {AppendTo[
            HDpoints, {Ga, HDfirst/Ha}], 
           AppendTo[HGpoints, {Ga, HGfirst/Ha}], 
           AppendTo[Hpoints, {Ga, Hfirst/Ha}], 
           AppendTo[Dpoints, {Ga, Dfirst/Da}], 
           AppendTo[Gpoints, {Ga, Gfirst/Ga}]}],
         If[
          0 <= HDsecond <= Da && 0 <= HGsecond <= Ga && 
           0 <= Dsecond <= Da && 0 <= Gsecond <= Ga && 
           0 <= Hsecond <= Ha && HGsecond <= Ha && HDsecond <= Ha && 
           Hsecond + HGsecond + HDsecond <= Ha && 
           Gsecond + HGsecond <= Ga && 
           Dsecond + HDsecond <= Da, {AppendTo[
            HDpoints, {Ga, HDsecond/Ha}], 
           AppendTo[HGpoints, {Ga, HGsecond/Ha}], 
           AppendTo[Hpoints, {Ga, Hsecond/Ha}], 
           AppendTo[Dpoints, {Ga, Dsecond/Da}], 
           AppendTo[Gpoints, {Ga, Gsecond/Ga}]}],
         If[
          0 <= HDthird <= Da && 0 <= HGthird <= Ga && 
           0 <= Dthird <= Da && 0 <= Gthird <= Ga && 
           0 <= Hthird <= Ha && HGthird <= Ha && HDthird <= Ha && 
           Hthird + HGthird + HDthird <= Ha && 
           Gthird + HGthird <= Ga && 
           Dthird + HDthird <= Da, {AppendTo[
            HDpoints, {Ga, HDthird/Ha}], 
           AppendTo[HGpoints, {Ga, HGthird/Ha}], 
           AppendTo[Hpoints, {Ga, Hthird/Ha}], 
           AppendTo[Dpoints, {Ga, Dthird/Da}], 
           AppendTo[Gpoints, {Ga, Gthird/Ga}]}]}, {Ga, 0.00000001, 
         Gmax, Gmax/displaypoints}], 1], 1], 1]]];
Manipulate[
 ListLinePlot[hdPts[Ha, Da, Gmax, pK1, pK2, displaypoints], 
  PlotRange -> {0, 1}],
 {Ha, 0, 1},
 {Da, 0, 1},
 {Gmax, 0, 1},
 {pK1, 0, 1},
 {pK2, 0, 1},
 {{displaypoints, 100}, 1, 200}]

As it is, it does not work, perhaps because of bad initial values. It may not be a good idea to use capital letters as variables, specially one like D, which has its own meaning.

3

Well, your problem is somewhat mysterious because you don’t say which are experimental inputs.

I assume Ha and Da are values which can be “fixed” experimentally. With this I do

eqs = {10^pK1 == HD/(H*D), 10^pK2 == HG/(H*G), Ha == HD + HG + H, Ga == HG + G, Da == HD + D};

Now I get rid of some parameters taking the experimentally fixed values into account:

In[6]:= t1 = Eliminate[eqs /. {Ha -> 1, Da -> 1}, {H, G, D, HG, Ha}]

Out[6]= -10^(2 pK1) + 10^pK1 HD + 3 10^(2 pK1) HD - 
  10^(pK1 + pK2) HD + 2^(1 + pK1 + pK2) 5^(pK1 + pK2) HD^2 - 
  10^pK1 HD^2 - 3 10^(2 pK1) HD^2 + 10^pK2 HD^2 + 10^(2 pK1) HD^3 - 
  10^(pK1 + pK2) HD^3 == 10^(pK1 + pK2) Ga (-1 + HD) HD

Then with t1 I define a function to get HD for a given Ga “knowing” that 0 < HD < 1 (because Ha == HD + HG + H and Ha == 1 )

Clear[ff]    
    ff[pK1_?NumericQ, pK2_?NumericQ, Ga_?NumericQ] := Module[{xx},
      xx = NSolve[
        -10^(2 pK1) + 10^pK1 HD + 3 10^(2 pK1) HD - 10^(pK1 + pK2) HD + 
          2^(1 + pK1 + pK2) 5^(pK1 + pK2) HD^2 - 10^pK1 HD^2 - 
          3 10^(2 pK1) HD^2 + 10^pK2 HD^2 + 10^(2 pK1) HD^3 - 
          10^(pK1 + pK2) HD^3 == 10^(pK1 + pK2) Ga (-1 + HD) HD, HD];
      Part[Select[Flatten[HD /. xx], 0 < # < 1 &], 1]
      ]

Test :
In[24]:= ff[2, 3, .6]

Out[24]= 0.455658

Then with the values of yours given above

data = {{1.*10^-8, 0.904875}, {0.05, 0.878427}, {0.1, 
    0.846593}, {0.15, 0.81062}, {0.2, 0.7719}, {0.25, 0.731591}, {0.3,
     0.690551}, {0.35, 0.649405}, {0.4, 0.608618}, {0.45, 
    0.568561}, {0.5, 0.529541}, {0.55, 0.491828}, {0.6, 
    0.455658}, {0.65, 0.42124}, {0.7, 0.38875}, {0.75, 
    0.358325}, {0.8, 0.330057}, {0.85, 0.30399}, {0.9, 
    0.280116}, {0.95, 0.25838}};


In[26]:= FindFit[data, ff[a, b, x], {a, b}, x]

Out[26]= {a -> 2., b -> 3.}

This should be checked with other combinations of parameters.

If possible (and you know him) give my greetings to Prof. Horst Hahn

4

Dear Hans Dolhaine,
Thank you very much for your help. This is what I was looking for and will help me with my work.

Can I ask you another short question about your code for which I couldn’t figure out the answer by browsing through the forum and documentation?

When I run your code as it is, it works nicely but it requires first evaluating

t1 = Eliminate[eqs /. {Ha -> 1, Da -> 1}, {H, G, D, HG, Ha}],

which gives

-10^(2 pK1) + 10^pK1 HD + 3 10^(2 pK1) HD - 10^(pK1 + pK2) HD + 
  2^(1 + pK1 + pK2) 5^(pK1 + pK2) HD^2 - 10^pK1 HD^2 - 
  3 10^(2 pK1) HD^2 + 10^pK2 HD^2 + 10^(2 pK1) HD^3 - 
  10^(pK1 + pK2) HD^3 == 10^(pK1 + pK2) Ga (-1 + HD) HD

and then to manually copy the t1 output into the xx equation that you have defined in the ff function.

 ff[pK1_?NumericQ, pK2_?NumericQ, Ga_?NumericQ] := Module[{xx},
      xx = NSolve[
        -10^(2 pK1) + 10^pK1 HD + 3 10^(2 pK1) HD - 10^(pK1 + pK2) HD + 
          2^(1 + pK1 + pK2) 5^(pK1 + pK2) HD^2 - 10^pK1 HD^2 - 
          3 10^(2 pK1) HD^2 + 10^pK2 HD^2 + 10^(2 pK1) HD^3 - 
          10^(pK1 + pK2) HD^3 == 10^(pK1 + pK2) Ga (-1 + HD) HD, HD];
      Part[Select[Flatten[HD /. xx], 0 < # < 1 &], 1]
      ]

Why can I not put “t1” directly into the ff function ?

ff[pK1_?NumericQ, pK2_?NumericQ, Ga_?NumericQ] := 
 Module[{xx}, xx = NSolve[t1, HD];
  Part[Select[Flatten[HD /. xx], 0 < # < 1 &], 1]]

Clearly, this doesn’t work because Select[Flatten[HD /. xx], 0 < # < 1 &]
gives an empty list. I suspect because Mathematica is evaluating 0 < # < 1 & before substituting the values for Ga, pK1 and pK2 inside.
Would you / anyone know how I can improve this part of the sample code that you gave? Thanks a lot.

Thanks again
Frank

(P.S. Yes, I know Horst Hahn as he is the Institute director of the INT, where my Emmy Noether junior research group is based. I will pass on your greetings when I see him. Did you have some joined projects in the past, and/or do you know him privately?)

5

OK, I “generalized” the Code a bit, so that it should work with other values for the experimentally fixed Ha and Da. See here and last part of the attachement.

Clear[t2]
t2[pK1_, pK2_, Ga_, Haa_, Daa_] := Evaluate[
  Eliminate[eqs /. {Ha -> Haa, Da -> Daa}, {H, G, D, HG, Ha}]]


Clear[f1]
f1[pK1_?NumericQ, pK2_?NumericQ, Ga_?NumericQ, Ha_?NumericQ, 
  Da_?NumericQ] := Module[{xx},
  xx = HD /. NSolve[t2[pK1, pK2, Ga, Ha, Da], HD];
  Part[Select[Flatten[xx], 0 < # < Ha &], 1]
  ]

In[14]:= FindFit[data, f1[a, b, x, 1, 1], {a, b}, x]

Out[14]= {a -> 2., b -> 3.}

But you should check it out with a couple of trials, meaning you should generate “experimental data” like your data above and check which answers are given by FindFit.

(And yes, years ago I worked together with Prof Hahn in a technical project and had some private contacts as well)

6

Thank you! Your generalized code works very well and evaluates fast. I included this into some larger routine for fitting of our signal - conc.(titrant) data and the fitting results are as they should have been. You saved us a lot of time as our previously preferred fitting programme (Origin) couldn’t handle some of the complex numbers resulting from the root-equations. With Mathematica, everything works like charm.