# Solve for several constants using a set of initial conditions

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.

• The parameter is suspicious. $k_3$ seems to be a negative number, but I think it should not be, right? – xzczd Oct 28 '17 at 11:32
• @xzczd, You are correct, k3 should not be negative. I have corrected a couple of the differential equation coefficients and updated the post. – razzle_dazzle_84 Oct 29 '17 at 2:59
• $k_3$ is still a negative number. You can use FindRoot to solve the equation about $k_3$. – xzczd Oct 29 '17 at 3:53
• 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. – Daniel Lichtblau Oct 29 '17 at 15:12
• @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. – razzle_dazzle_84 Oct 29 '17 at 23:09

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.

• 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. – razzle_dazzle_84 Nov 4 '17 at 5:42
• @razzle_dazzle_84 Oh, I forgot I've modified the 2nd argument of DSolve. See my edit. – xzczd Nov 4 '17 at 6:36
• 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. – razzle_dazzle_84 Nov 6 '17 at 15:39

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.

• 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. – Daniel Lichtblau Nov 1 '17 at 3:47
• 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? – razzle_dazzle_84 Nov 4 '17 at 5:48
• 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. – Bill Watts Nov 4 '17 at 6:03
• 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. – Daniel Lichtblau Nov 6 '17 at 17:09