Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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!

share|improve this question
its doable in principle but right now your set of equations does not seem to integrate properly. You forgot to give k1 a value – chris Nov 4 '12 at 20:23
Hi Chris, Thanks for your reply. I'm sorry that I forgot to include a value for k1. K1& k2 values are:k1 = 7.32*10^-5; k2 = 8.09*10^-9; – DumbleKo Nov 4 '12 at 21:46
Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2)Read the FAQs! 3) When you see good Q&A, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. ALSO, remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign` – chris Nov 5 '12 at 20:15
up vote 5 down vote accepted

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

chi2[k1_?NumberQ, k2_?NumberQ] :=
 (sol = 
   Caf908 /. NDSolve[
      Equation[k1, k2], {Caf908, Cbf908, Hf908}, {t, 0, 1000}] // 
  (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

Mathematica graphics

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.


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

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

Mathematica graphics

and Indeed the fit seems okish:

 chi2[kk1*1.3, kk2, True]

Mathematica graphics

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

share|improve this answer
Thanks for your help, Chris. Do you know if there is anyway I can find k1&k2 if the values are incorrect? I just started learning mathematica, and know there are findfit and fit. But I must have a function in order to use these two method, and I only have interpolating function from NDSolve. Could you show me how to use NMinimize? Thanks! – DumbleKo Nov 4 '12 at 21:50
Hi Chris, the problem is for solving chemical reaction, and k1&k2 are the reaction constants of forward and reverse reactions. I'm still having some problems though. When I use k1&k2 (k1->15.3735,k2->2.45184*10^-9) obtained from NMinimize to fit my data, it gave me a very bad fit, while the previous values (k1 = 7.32*10^-5; k2 = 8.09*10^-9) fits better still. Thanks. – DumbleKo Nov 5 '12 at 0:04
Thank you so much, Chris! I finally solved the problem after fixing some values and my equations!! I'll let you know if I have more questions! – DumbleKo Nov 5 '12 at 20:13
@chris What is the //#.#& at the end of the chi2 function? – Christopher Bowman Mar 21 '13 at 18:45
NVM. I just realized it's the square of the residuals between what's predicted by the interpolation and the actual data. – Christopher Bowman Mar 21 '13 at 19:03

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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