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.
FindRoot
to solve the equation about $k_3$. $\endgroup$