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}]

ParametricNDSolve
,FindFit
and(List)Interpolation
. $\endgroup$ – user21 Aug 6 '13 at 15:41