Solve for several constants using a set of initial conditions

Question

I have a set of differential equations that I am trying to solve for chemical reaction rate constants. So far I can get DSolve to solve the ODEs with initial conditions which provides the chemical species' concentration equations in terms of the rate constants k1, k2, k3, k4 and time. I can then solve for the first rate constant (k1) using Solve and a known value of species A at t=60. However, this result seems to be more symbolic as k1 is given as a number divided by 't'.

What I would like to have happen is solve the ODEs, then solve for all the rate constants so that I get numeric values for k1, k2, k3 and k4. Here is what I have so far:

eqns = {a'[t] == -3*Subscript[k, 1]*a[t], 
b'[t] == 2/3*Subscript[k, 1]*a[t] - Subscript[k, 2]*b[t], 
c'[t] == 3*Subscript[k, 2]*b[t] - Subscript[k, 3]*c[t], 
d'[t] == Subscript[k, 3]*c[t] - 3*Subscript[k, 4]*d[t], 
e'[t] == 1/3*Subscript[k, 4]*d[t]}

initcond1 = {a[0] == 3.8523, b[0] == 0, c[0] == a[0]*(2*10^-9)*700^2.6857, 
d[0] == 0, e[0] == 0}

endcond1 = {0.1093, 0.587, 4.763, 0.7754, 0.3484}

sys = First@DSolve[{eqns,initcond1},{a[t],b[t],c[t],d[t],e[t]},t]

a[t_] := sys[[1,2]]
Solve[a[60] == 0.1093, Subscript[k, 1]]

Which gives me the output:

{{Subscript[k, 1] -> 1.18744/t}}

Which I can turn into the correct answer by dividing 1.18744 by 60, but I would rather get the result:

{Subscript[k, 1] -> <a number>, Subscript[k, 2] -> <another number>, Subscript[k, 3] -> <yet another number>, Subscript[k, 4] - <one last number>}

I have spent a great deal of time tonight looking at the docs for Solve, DSolve and NDSolve and many forum posts, but I am still missing something. Any help is appreciated.

UPDATE I can now get the solution for k1 & k2 without the delayed set command:

round1 = Solve[{sys[[1, 2]] == endcond1[[1]], sys[[2, 2]] == endcond1[[2]]} /. t -> 60, {Subscript[k, 1], Subscript[k, 2]}]

Then I globally assign the solutions to the k1 and k2 values and try to solve for k3:

{Subscript[k, 1], Subscript[k, 2]} = {round1[[1, 1, 2]], round1[[1, 2, 2]]}

round2 = Solve[{FullSimplify[sys[[3,2]] /. t->60] == endcond1[[3]], Subscript[k, 3]>0}, Subscript[k, 3]]

However, the two solutions returned for k3 are 3*k1 and k2 which break the Solve for k3 since the denominator of c[t] is (-3*k1+k3)*(-k2+k3). I tried specifying that k3 could not equal 3*k1 or k2:

round2 = Solve[{FullSimplify[sys[[3,2]] /. t->60] == endcond1[[3]], Subscript[k, 3]>0, Subscript[k, 3] != 3*Subscript[k, 1], Subscript[k, 3] != Subscript[k, 2]}, Subscript[k, 3]]

But I get the same results with the warning:

Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. 

So...I would still like to find a more elegant way to solve the entire 'sys' system using the 'endcond1' array. If the solution for k3 keeps puking, then either my differential equation coefficients are incorrect or my endcond1 values are bad.

Answer 1

Since you just want to know "if there is a way to solve the system with a single command or 2 instead of piecemeal", here's the answer:

sys = First@DSolve[{eqns, initcond1}, {a, b, c, d, e}, t];

endcond1 = {a[60] == 0.1093, b[60] == 0.587, c[60] == 4.763, d[60] == 0.7754, 
   e[60] == 0.3484};

Reap[FoldList[#[[2 ;;]] /. (FindRoot[#[[1]], {#2, 0}][[1]] // Sow) &, endcond1 /. sys, 
   Subscript[k, #] & /@ Range@5];]

Other part of the code is the same as yours.

BTW, I believe if the parameters are correctly given, it won't be necessary to solve the equation system one by one. FindRoot should be able to deal with them all at once.

Answer 2

I believe either the data or the equations are inconsistent. I did the following using the same data as given.

sys = First@
DSolve[{eqns, initcond1}, {a[t], b[t], c[t], d[t], e[t]}, t] // Simplify;

a[t_] = a[t] /. sys;
b[t_] = b[t] /. sys;
c[t_] = c[t] /. sys;
d[t_] = d[t] /. sys;
e[t_] = e[t] /. sys;

Solving for the k's involve transcendental equations, so FindRoot or some other numerical method is generally required.

values = {};
values = Join[values, 
   FindRoot[(a[60] /. values) - endcond1[[1]], {Subscript[k, 1], 1/10}] // Flatten];
values = Join[values, 
  FindRoot[(b[60] /. values) - endcond1[[2]], {Subscript[k, 2],1/10}] // Flatten];
values = Join[values, 
  FindRoot[(c[60] /. values) - endcond1[[3]], {Subscript[k, 3], 1/10}] // Flatten];
values = Join[values, 
  FindRoot[(d[60] /. values) - endcond1[[4]], {Subscript[k, 4], 50}] // Flatten]

{Subscript[k, 1]->0.0197907,Subscript[k, 2]->0.00788753,Subscript[k, 3]->-0.0346387,Subscript[k, 4]->7093.36}

The value for k4 is bogus. A plot reveals the line for the equations never crosses 0, probably due to a negative k3.

A negative k no doubt indicates violation of conservation of mass. The end value of c3 is larger than the beginning values of all the products combined which probably accounts for the negative k3.

I could have used the final condition of e to find k4. Consistent data should give the same answer as the d final condition, which is not the case here.

Using a smaller value for the final c3 gave me reasonable k's, but of course, the final value of e did not match.

In the equations, it also seems strange that a is being reduced by 3 k1 a, while b is growing at only 2/3 k1 a, although not knowing the reactions involved, I suppose it is possible for it to take 9 a's to make 2 b's.