This Answer has been revised in two respects: First, boundary condition values are Rationalize
d 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.