How do I find the best parameter to fit my data if the model is a interpolating function?

Question

Hi I have a question regarding to find the best parameters for my model to fit my data. I have 3 ordinary equation, and I now just picked some parameters (k1 = 7.32*10^-5; k2 = 8.09*10^-9) for my ODE.

Defined all the constant values:

k1 = 7.32*10^-5;
k2 = 8.09*10^-9;
Da = 1.33*10 - 4;
Nd0 = 0.001;
HR0 = 0.2;
H0908 = 0.004;
x0 = 0.007;

And define the 3 ODE equations:

Equation1908 = {Nd0*
     Caf908'[t] == -k1*(Nd0*HR0)/H0908*Caf908[t]*Cbf908[t]/Hf908[t] + 
     k2*(0.2 - 3*HR0*Cbf908[t]),
   HR0*Cbf908'[t] == -3*k1*(Nd0*HR0)/H0908*Caf908[t]*
      Cbf908[t]/Hf908[t] + 3*k2*(0.2 - 3*HR0*Cbf908[t]),
   H0908*Hf908'[t] == 
    3*k1*(Nd0*HR0)/H0908*Caf908[t]*Cbf908[t]/Hf908[t] - 
     3*k2*(0.2 - 3*HR0*Cbf908[t]),
   Caf908[0] == 1, Cbf908[0] == 1, Hf908[0] == 1};

Using NDSolve:

BC1908 = NDSolve[
  Equation1908, {Caf908[t], Cbf908[t], Hf908[t]}, {t, 0, 1000}]

At the end, I obtained 3 interpolating function, for Caf908[t], Cbf908[t], and Hf908[t]... I used the solution (Caf908[t]) to fit my data. X-axis (tfexp0908) versus y-axis (fexp0908):

fexp0908 = .001*{1.009790523`, 0.898335138`, 0.878948419`, 
    0.830114856`, 0.767123385`, 0.732170062`, 0.672106602`, 
    0.637589428`, 0.59141947`, 0.523944512`, 0.554584169`, 
    0.451444203`, 0.396545111`, 0.444125908`, 0.352355452`, 
    0.272913948`, 0.33877861`, 0.287900412`, 0.276936425`};
tfexp0908 = {0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 180, 210, 
   245, 270, 300, 330, 360, 390};

And I obtained a good fit. But I would like to know how I can use Mathematica to obtain the best parameter, k1 & k2, for my model (Caf908[t]) and get the best fit to my data. Thank you!

Answer 1

Let's clear the k1 and k2 variables and put in kk1 and kk2 the corresponding values

Clear[k1, k2];
kk1 = 7.32*10^-5;
kk2 = 8.09*10^-9;

Now define a set of equations to solve depending on k1,k2

Equation[k1_?NumberQ,k2__?NumberQ] = {Nd0*
     Caf908'[t] == -k1*(Nd0*HR0)/H0908*Caf908[t]*Cbf908[t]/Hf908[t] + 
     k2*(0.2 - 3*HR0*Cbf908[t]), 
     HR0*Cbf908'[t] == -3*k1*(Nd0*HR0)/H0908*Caf908[t]*
      Cbf908[t]/Hf908[t] + 3*k2*(0.2 - 3*HR0*Cbf908[t]), 
     H0908*Hf908'[t] == 
      3*k1*(Nd0*HR0)/H0908*Caf908[t]*Cbf908[t]/Hf908[t] - 
      3*k2*(0.2 - 3*HR0*Cbf908[t]), Caf908[0] == 1, Cbf908[0] == 1, 
   Hf908[0] == 1};

Now let us define a $\chi^2$ function

Clear[chi2];
chi2[k1_?NumberQ, k2_?NumberQ] :=
 (sol = 
   Caf908 /. NDSolve[
      Equation[k1, k2], {Caf908, Cbf908, Hf908}, {t, 0, 1000}] // 
    First;
  (sol /@ tfexp0908) - fexp0908 // #.# &
  )

and let's sample this function around what is supposedly the correct value

Table[chi2[kk1*i, kk2*j], {i, 0.1, 2.1, 0.1}, {j, 0.1, 5.1, 
   0.2}] // ListContourPlot

It seems to suggest i) the guess is not so good ii) somehow the $\chi^2$ does not depend on k2

Finding the right range of values for kk1 and kk2 and sorting out why the $\chi^2$ does not depend on k1 is beyond a Mathematica question?

Note that one could use in principle NMinimize to find the solution more efficiently than via sampling the likelihood function, but before these points are addressed there's no point.

EDIT

I believe you should remove the 0.01* in front of fexp0908

Clear[chi2];
chi2[k1_?NumberQ, k2_?NumberQ, 
  debug_: False] := (sol = 
   Caf908 /. 
     NDSolve[Equation[k1, k2], {Caf908, Cbf908, Hf908}, {t, 0, 
       1000}] // First;
  If[debug, {(sol /@ tfexp0908), fexp0908} // ListLinePlot // Print];
  (sol /@ tfexp0908) - fexp0908 // #.# &)

because without it I get

 Table[{i*kk1,chi2[kk1*i, kk2*j]},{i, 0.9, 2.3, 0.1},
    {j, 0.1, 5.1, 0.2}]//Transpose//ListLinePlot

and Indeed the fit seems okish:

 chi2[kk1*1.3, kk2, True]

Now the minimization

NMinimize[{chi2[k1, k2], k1 > 0, k2 > 0}, {k1, k2}]

(* {0.0339566,{k1->0.000110757,k2->0.} *)

yields a poorer fit but hasn't converged.

It seems your set of data poorly constrains k2