Matrix Exponentiation
Algorithms for the exponentiation of a matrix have become available in recent years. It is not generally known that matrix exponentiation may be used to solve sets of kinetic differential equations where all rates are first order. Most kinetic processes of interest consist of a collection of first and pseudofirst order processes. Consider the simple first order decay:
dA/dt = rate(1)*A = kA.
The solution to this is A = Ao.e(kt), i.e., the solution is of exponential form.
Consider now the case of A → B, for which the differential equations are:
dA/dt = kA
dB/dt = kA,
which can be written in matrix form as:

A 
B 
dA/dt 
k 
0 
db/dt 
+k 
0 
Here, the matrix columns represent the respective initial concentrations and the rows their derivatives with respect to time. Note that in this case the initial concentration of B is zero, and that the B column entries are zero because there is no back reaction.
In general such a matrix contains only rate constants and zeros. The available algorithms for calculation of the exponential of a matrix produce a solution which, when applied to the above matrix, corresponds to:
This matrix resulting from the exponentiation process requires further simplification. Recalling that exponentials are normalized to unity the initial concentrations of A and B are the 1 and 0. Accordingly the resultant matrix is multiplied by a vector with contents of 1 and 0 yielding the concentrations for A and B at time t as exp(kt) and 1exp(kt). This result is the same is obtained by the usual algebraic methods and for a case of N time points yields an Nx2 concentration matrix. This concentration matrix is then inverted or used to calculate the function as required.
Thus the differential equations have been solved by applying a special type of exponentiation to a matrix of rate constants representing the differential equations. Matrix Exponentiation is thus a more general kind of exponentiation: the beauty of the Exponential Matrix method is that it may be easily extended to much more complicated cases with many first order processes including reversible steps. Application to complex cases in practice is not computationally limited but is strongly limited in terms of signal to noise ratio of the data. In other word, increasing the number of significantly determined rates requires an increasingly higher signal to noise ratio.
The Exponential Matrix method is markedly superior to the numerical integration method of RungeKutta for the computation of complex arrays of first order processes. It may be readily shown for the simple and common case of A → B → C the exponential matrix method is more than 20 times faster than RungeKutta. Thus not only does RungeKutta have well known stiffness problems it is so much slower that it makes the program no longer real time.
Formation of Differential Matrix
Mechanism
Rate equation in differential form
dA/dt = k_{1}(A) 
+0 
+0 
dB/dt = +k_{1}(A) 
k_{2}(B) 
+0 
dC/dt = 0 
+k_{2}(B) 
+0 
Matrix derived from above by removing the concentration terms
k1 
0 
0 
+k1 
k2 
0 
0 
+k2 
0 
It is this matrix itself which is solved by Matrix Exponentiation!
In this way, Matrix Exponentiation provides an exact solution of the rate equation.
We could say "instantaneous" and be fully accurate!
