Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I am trying to solve following system of three differential equations

r = {
  -cab'[t] == k*cab[t] - kb*ca[t]*cb[t],
  ca'[t] == k*cab[t] - kb*ca[t]*cb[t],
  cb'[t] == k*cab[t] - kb*ca[t]*cb[t],
  cab[0] == 0.071176,
  ca[0] == 0.003697,
  cb[0] == 0.002845
};

and fit them to the following experimental data

cab = {0.071176, 0.05493, 0.04287, 0.03391, 0.02767, 0.02355, 
       0.02096,0.01825, 0.017029, 0.016582, 0.01421, 0.01432, 0.01431};
ca = {0.003697, 0.021051, 0.033007, 0.04154, 0.047609,0.052086, 0.05529, 
      0.057088, 0.05842, 0.05966, 0.06114, 0.06079,0.06153};
cb = {0.002845, 0.01949, 0.03078, 0.03714, 0.04476, 0.04939, 0.05214, 
      0.05327, 0.05526, 0.05573, 0.05859, 0.05774, 0.05868};
t = {0, 3600, 7200, 10800, 14400, 18000, 21600, 25200, 28800, 32400,82800, 86400, 90000};

I want to get the best k and kb values. Could you help me since I am novice in Mathematica. I would like a step-by-step tutorial on hoe to do it.

share|improve this question
    
welcome to SE. It is customary here to show code that has been written so far, could you share that? Otherwise, if you are looking for a staring point, have a look at ParametricNDSolve, FindFit and (List)Interpolation. –  user21 Aug 6 '13 at 15:41
add comment

2 Answers

Rather than looking for a solution of the differential equations including the initial conditions, I'd postpone imposing them to when one does the fit. The solution will contain 3 constants to be determined based on the initial conditions :

{solca[t_], solcab[t_], solcb[t_]} = 
  DSolve[{-cab'[t] == k*cab[t] - kb*ca[t]*cb[t], 
            ca'[t] == k*cab[t] - kb*ca[t]*cb[t], 
            cb'[t] == k*cab[t] - kb*ca[t]*cb[t]}, 
    {cab[t], ca[t], cb[t]}, t][[1, All, 2]]

You might notice that the three solutions only differ by a constant (and a sign), which is obvious when looking at your system. Our model is a function which reproduces each solution to the system for index=1,2 or 3; we also need to rename the integration constants so we can manipulate them later.

f[index_, t_] = KroneckerDelta[index, 1] solca[t] + 
                KroneckerDelta[index, 2] solcab[t] + 
                KroneckerDelta[index, 3] solcb[t] /. {C[1] -> c1, C[2] -> c2, C[3] -> c3};

The next step would be to use FindFit to get the parameters k, kb; however they appear in 3 different functions and we want to have a common fit. I'll use a trick I learned here.

datacab = {0.071176, 0.05493, 0.04287, 0.03391, 0.02767, 0.02355, 0.02096, 0.01825, 0.017029, 0.016582, 0.01421, 0.01432, 0.01431};
dataca = {0.003697, 0.021051, 0.033007, 0.04154, 0.047609, 0.052086, 0.05529, 0.057088, 0.05842, 0.05966, 0.06114, 0.06079, 0.06153};
datacb = {0.002845, 0.01949, 0.03078, 0.03714, 0.04476, 0.04939, 0.05214, 0.05327, 0.05526, 0.05573, 0.05859, 0.05774, 0.05868};
datat = {0, 3600, 7200, 10800, 14400, 18000, 21600, 25200, 28800, 32400, 82800, 86400, 90000};

data = Join[Thread[List[1, datat, dataca]], 
            Thread[List[2, datat, datacab]], 
            Thread[List[3, datat, datacb]]]

Now data is in a suitable form to be taken as our experimental data for our model. We can also impose the initial conditions as part of the fitting.

solfit = 
FindFit[data, {f[index, t], f[2, 0] == 0.071176, f[1, 0] == 0.003697, f[3, 0] == 0.002845}, 
 {k, kb, c1, c2, c3}, {index, t}, Method -> NMinimize]
(* {k -> 0.298832, kb -> 3.07733, c1 -> 0.074873, c2 -> -0.00085201, c3 -> -1.95533} *)

Check :

f[1, 0] /. solfit // Chop
f[2, 0] /. solfit // Chop
f[3, 0] /. solfit // Chop
(* 0.003697
   0.071176
   0.00284499 *)

Finally we can put model and best fit together, do some manipulations to remove apparent imaginary parts and define our final solution ;

solution[index_, t_] = 
 Simplify[ComplexExpand[f[index, t] /. solfit, TargetFunctions -> {Re, Im}], 
  Assumptions -> {KroneckerDelta[__] \[Element] Reals}]

Plot[{solution[1, t], solution[2, t], solution[3, t]}, {t, 0, 10}]

plot

share|improve this answer
    
Thanks a lot for your reply. I have tried to implement your solution (simply buy copy and paste method) but I had some errors. –  Lolek Aug 7 '13 at 10:05
    
Could you try with a fresh kernel and be more specific about the errors ? –  b.gatessucks Aug 7 '13 at 10:12
    
I tried with a fresh kernel and I got the same results as you. But values of k = 0.29 and kb = 3.077 are not what I expected. This plot which can be seen does not fit to my experimental data except for the first point at x=0. There is some thing wrong with x-axis values. Please have a look at the datat values (time). This values are from 0 to 90000 while on the plotted graph these values are between 0 and 10. Is it possible to overlay my experimental data with the graph being obtained from calculations to see what is the correlation? –  Lolek Aug 7 '13 at 11:38
    
Nice to see another implementation of indexing multiple functions for solving a system of equations. @oleksandr-r used a slightly different method here just for reference. –  bobthechemist Aug 7 '13 at 12:09
    
@Lolek You can use solfit2 = NonlinearModelFit[ data, {f[index, t], f[2, 0] == 0.071176, f[1, 0] == 0.003697, f[3, 0] == 0.002845}, {k, kb, c1, c2, c3}, {index, t}, Method -> NMinimize] and look at all the details of the fit, for instance solfit2["ParameterConfidenceIntervalTable"] and solfit2["FitResiduals"]. –  b.gatessucks Aug 7 '13 at 14:11
add comment

I have modyfied my script. But now, the fitting is done taking into account only one set of date - datacab. Do you know how to change it so that minimization was done based on 3 sets of data (datacab, dataca, datacb)

datacab = {0.071176, 0.05493, 0.04287, 0.03391, 0.02767, 0.02355, 0.02096, 0.01825,
0.017029, 0.016582, 0.01421, 0.01432, 0.01431};
dataca = {0.003697, 0.021051, 0.033007, 0.04154, 0.047609, 0.052086, 0.05529,
0.057088, 0.05842, 0.05966, 0.06114, 0.06079, 0.06153};
datacb = {0.002845, 0.01949, 0.03078, 0.03714, 0.04476, 0.04939, 0.05214, 0.05327,
0.05526, 0.05573, 0.05859, 0.05774, 0.05868};
datat = {0, 3600, 7200, 10800, 14400, 18000, 21600, 25200, 28800, 32400, 82800,
86400, 90000};
ll = Length[datat];
dane = ListPlot[Table[{{datat[[i]], datacab[[i]]}, {datat[[i]], 
 dataca[[i]]}, {datat[[i]], datacb[[i]]}},{i, 1, ll}], PlotStyle -> PointSize[0.02]];
Clear[k1, k2];
kk1 = 0.00007;
kk2 = 0.000274;
Equation[k1_?NumberQ, k2__?NumberQ] = {-cab'[t] == k1*cab[t] - k2*ca[t]*cb[t], ca'[t] == k1*cab[t] - k2*ca[t]*cb[t], cb'[t] == k1*cab[t] - k2*ca[t]*cb[t], cab[0] == 0.071176,ca[0] == 0.003697, cb[0] == 0.002845};
Clear[chi2];
chi2[k1_?NumberQ, k2_?NumberQ] := (sol = 
cab /. NDSolve[Equation[k1, k2], {cab, ca, cb}, {t, 0, 90000}] // 
 First; (sol /@ datat) - datacab // #.# &);
Clear[k1, k2]
st = NMinimize[{chi2[k1, k2], k1 > 0, k2 > 0}, {k1, k2}]
k1 = st[[2, 1, 2]];
k2 = st[[2, 2, 2]];
eq = {-cab'[t] == k1*cab[t] - k2*ca[t]*cb[t], ca'[t] == k1*cab[t] - k2*ca[t]*cb[t], cb'[t] == k1*cab[t] - k2*ca[t]*cb[t], cab[0] == 0.071176, ca[0] == 0.003697, cb[0] == 0.002845};
sol = NDSolve[eq, {cab, ca, cb}, {t, 0, 90000}];
ff = Plot[Evaluate[{cab[t], ca[t], cb[t]} /. sol], {t, 0, 90000}, PlotStyle -> {Red, Thick}];
Show[ff, dane]
share|improve this answer
    
This shouldn't be posted as an answer. It's a follow-up question. Maybe you want to edit your question instead, and post a comment to @b.gatessucks to let him know of the change. –  Jens Aug 8 '13 at 18:11
add comment

Your Answer

 
discard

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.