Solving ODES in order to find unknown constants

Question

I am trying to solve a set of ODES which represent a second order consecutive chemical reaction.

eqn = {
   CA'[t] == -k1*CA[t]*CF[t],
   CB'[t] == k1*CA[t]*CF[t] - k2*CB[t]*CF[t],
   CC'[t] == k2*CB[t]*CF[t] - k3*CC[t]*CF[t],
   CD'[t] == k3*CC[t]*CF[t] - k4*CD[t]*CF[t],
   CE'[t] == k4*CD[t]*CF[t],
   CF'[t] == -(k1*CA[t]*CF[t] + k2*CB[t]*CF[t] + k3*CC[t]*CF[t] + 
       k4*CD[t]*CF[t])};

I want to solve these equations either for a generic expression for k1, k2, k3, or k4, or for a particular value with these particular boundary conditions, n.b. not to sure if both t= 1 and t = 4 needs to be used...

bcs = {
   CA[0] == 0.052952,
   CB[0] == 0,
   CC[0] == 0,
   CD[0] == 0,
   CE[0] == 0,
   CF[0] == 0.0226939,
   CA[1] == 0.047442,
   CB[1] == 4.08*10^-3,
   CC[1] == 1.24*10^-3,
   CD[1] == 0.19*10^-3,
   CE[1] == 0,
   CF[1] == 0.0191289,
   CA[4] == 0.043412,
   CB[4] == 5.03*10^-3,
   CC[4] == 3.54*10^-3,
   CD[4] == 0.90*10^-3,
   CE[4] == 0.07*10^-3,
   CF[4] == 0.0151489};

I have tried to use both NDsolve and Dsolve, NDSolve requires that I have initial values for k1, k2, k3, k4 and Dsolve doesn't seem to want to help... the rest of my code for DSolve...

vars = {CA[t], CB[t], CC[t], CD[t], CE[t], CF[t]};

solution[t_] = DSolve[{eqn, bcs}, vars, t]

Any help would be appreciated!! Chemistry is my forte not Mathematica!!

Cheers!

Answer 1

This Answer has been revised in two respects: First, boundary condition values are Rationalized to postpone roundoff issues until the very end of the calculation. Second, equations ans3 are solved for k1 τ1 instead of Exp[-k1 τ1], etc.

Define

var = {CA[t], CB[t], CC[t], CD[t], CE[t], CF[t]};
bc0 = {CA[0] == 0.052952, CB[0] == 0, CC[0] == 0, 
  CD[0] == 0, CE[0] == 0, CF[0] == 0.0226939}/. z_Real :> Rationalize[z, 0];
bc1 = {CA[1] == 0.047442, CB[1] == 4.08*10^-3, CC[1] == 1.24*10^-3, 
  CD[1] == 0.19*10^-3, CE[1] == 0, CF[1] == 0.0191289}/. z_Real :> Rationalize[z, 0];

and eqn as in the Question.

These equations become linear, if CF[t] can be eliminated from the right side. This is accomplished by the change of independent variable t to τ such that D[τ[t] ,t] == CF[t]. The equations then become

eqn /. CF[t] -> 1 /. t -> τ;

which are solved by

ans1 = First@DSolve[{eqn /. CF[t] -> 1, bc0} /. t -> τ, var /. t -> τ, τ]
(* {CA[τ] -> 6619/(125000*E^(k1*τ)), CB[τ] -> 
  (-6619*E^(-(k1*τ) - k2*τ)*(-E^(k1*τ) + E^(k2*τ))*k1)/(125000*(k1 - k2)), ...} *)

To apply bc1, first isolate the four exponentials

ans2 = Collect[Expand[ans1], {Exp[-k1 τ], Exp[-k2 τ], Exp[-k3 τ], Exp[-k4 τ]},
  Simplify] /. Rule -> Equal;

and evaluate ans2 at τ1, corresponding to t = 1.

ans3 = Simplify[ans2 /. τ -> τ1 /. (bc1 /. Equal -> Rule /. (1 -> τ1)) 
         /. {k1 -> k1/τ1, k2 -> k2/τ1, k3 -> k3/τ1, k4 -> k4/τ1}]
(* {26476/E^k1 == 23721, 
     ((510 + 6619/E^k1 - 6619/E^k2)*k1 - 510*k2)/(k1 - k2) == 0, ...} *)

Note the substitution {k1 -> k1/τ1, k2 -> k2/τ1, k3 -> k3/τ1, k4 -> k4/τ1}, which has the effect of eliminating the unknown quantity, τ1, from ans3. The first four equations now can be solved for the four k. Because NSolve can handle only the first two transcendental equations, we proceed as follows.

a1 = First@NSolve[ans3[[1]], k1, Reals]
(* {k1 -> 0.109878} *)
a2 = First@NSolve[ans3[[2]] /. a1, k2]
(* {k2 -> 0.621713} *)
a3 = FindRoot[ans3[[3]] /. a1 /. a2, {k3, .4}]
(* {k3 -> 0.417505} *)
a4 = FindRoot[ans3[[4]] /. a1 /. a2 /. a3, {k4, 0}]
(* {k4 -> 2.28314*10^-15} *)

(Probably k4 should be zero, but this is the best that FindRoot can do here.) Inserting these k values into ans3[5] leaves a residual of order 10^-17, essentially zero. Why this should be expected is explained later. On the other hand, inserting these k values into ans3[6] leaves a residual of 0.003565. Thus, the system of equations defined in the Question has no solution.

Minor modifications to the boundary conditions can, of course, make the system consistent. For instance, subtracting 0.003565 from bc1[6] does so.

Finally, we remark that summing the first five eqn and their corresponding boundary conditions yields

{Simplify[Total[Cases[Take[eqn, 5], (z_ == _) -> z, {0, Infinity}]]] ==
   Simplify[Total[Cases[Take[eqn, 5], (_ == z_) -> z, {0, Infinity}]]], 
   Total[Cases[Take[bc0, 5], (_ == z_) -> z, {0, Infinity}]], 
   Total[Cases[Take[bc1, 5], (_ == z_) -> z, {0, Infinity}]]}
(* {Derivative[1][CA][t] + Derivative[1][CB][t] + Derivative[1][CC][t] + 
    Derivative[1][CD][t] + Derivative[1][CE][t] == 0, 6619/125000, 6619/125000}

Hence,

CA[t]  + CB[t]  + CC[t] + CD[t] + CE[t] == 6619/125000

which indicates that one of the first five eqn is redundant and could have been deleted. In fact, CE is the logical variable to drop.

Numerical Solution

Once the ratios of the K are known, it is straightforward to integrate the original eqn and from there determine τ1 = 52.3361. Rescaling the k by τ1 then gives

Solutions at t = 1 differ from bc1 by of order 10^-9, except for CF[1], which differs by 0.003565, as expected.