2
$\begingroup$

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.

$\endgroup$
10
  • $\begingroup$ The parameter is suspicious. $k_3$ seems to be a negative number, but I think it should not be, right? $\endgroup$
    – xzczd
    Oct 28, 2017 at 11:32
  • $\begingroup$ @xzczd, You are correct, k3 should not be negative. I have corrected a couple of the differential equation coefficients and updated the post. $\endgroup$ Oct 29, 2017 at 2:59
  • $\begingroup$ $k_3$ is still a negative number. You can use FindRoot to solve the equation about $k_3$. $\endgroup$
    – xzczd
    Oct 29, 2017 at 3:53
  • $\begingroup$ There appears to be something wrong with the system itself. At steady state four variables (a,b,c,d) are zero no matter how the initial conditions are set. $\endgroup$ Oct 29, 2017 at 15:12
  • $\begingroup$ @DanielLichtblau, this system is similar to that developed by Mondal et al. in this article. My system models the reduction of Fe2O3 to iron by CO through the intermediaries Fe3O4 and FeO. The creation of Fe3C is also part of the system. I would appreciate an answer to my question of how to expediently solve the system for the k-values. $\endgroup$ Oct 29, 2017 at 23:09

2 Answers 2

0
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ I can't get this to work. I entered in my eqns and initcond1 array as before, then ran the 'sys' solver and finally ran the two lines of code you provided. I receive several errors for FindRoot, ReplaceAll and General. $\endgroup$ Nov 4, 2017 at 5:42
  • $\begingroup$ @razzle_dazzle_84 Oh, I forgot I've modified the 2nd argument of DSolve. See my edit. $\endgroup$
    – xzczd
    Nov 4, 2017 at 6:36
  • $\begingroup$ Works great now. Thanks! It will take me a while to understand the syntax of the final line, but this is the answer I was looking for. $\endgroup$ Nov 6, 2017 at 15:39
2
$\begingroup$

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.

$\endgroup$
4
  • 1
    $\begingroup$ I got those same unphysical values for the reaction rates. Gave up when I saw it was not possible for me to get anything meaningful from the system. $\endgroup$ Nov 1, 2017 at 3:47
  • $\begingroup$ I have corrected a couple of coefficients. Now: db/dt=2*k1*a-k2*b and dd/dt=k3*c-k4*d. Running your solution with the rest of my code the same provides nice k values: {Subscript[k, 1] -> 0.0197907, Subscript[k, 2] -> 0.0350742, Subscript[k, 3] -> 0.00782303, Subscript[k, 4] -> 0.0282808}`. What is the purpose of the "1/10" and "50" values, though? $\endgroup$ Nov 4, 2017 at 5:48
  • $\begingroup$ They are initial guesses for each of the k's as required by FindRoot. If you know what the approximate values of the k's should be, you should use those. You can also use a range. Look up the help for FindRoot. $\endgroup$
    – Bill Watts
    Nov 4, 2017 at 6:03
  • $\begingroup$ This correction needs to be edited into the original post so that people do not waste time, like I just did, cutting and pasting from there under the misguided assumption that the correct values were now in place. $\endgroup$ Nov 6, 2017 at 17:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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