2
$\begingroup$

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!

$\endgroup$
  • $\begingroup$ Welcome to Mathematica.SE! I suggest that: 1) You take the introductory Tour now! 2) When you see good questions and answers, 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, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! 3) As you receive help, try to give it too, by answering questions in your area of expertise. $\endgroup$ – bbgodfrey Apr 7 '15 at 2:05
  • $\begingroup$ I start with a simpler problem. eqn = {CA'[t] == -k1*CA[t]*CF[t], CB'[t] == k1 CA[t] CF[t] - k2 CB[t] CF[t], CF'[t] == -(k1*CA[t]*CF[t])}; vars = {CA[t], CB[t], CF[t]}; DSolve[eqn, vars, t] and that solves. Then I try adding initial conditions, one at a time, and find one or two that still solve, but much slower. Some parts of DSolve want exact rational coefficients. Other parts will convert decimal to rational, complain, but keep going. See if you can make any progress with this partial system. $\endgroup$ – Bill Apr 7 '15 at 3:47
2
$\begingroup$

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

enter image description here

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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