Gauss-Jordan elimination gives inconsistent matrix for a consistent system?

Question

I am trying to get an analytical expression for a steady state of an ODE system governing a chemical reaction network via symbolic computer algebra systems. As an example for this question I'll take a fully connected chemical reaction network with three species $N_0, N_1, N_2$ and three bidirection reactions $N_0 \rightleftarrows^{k_1}_{k_2} N_1 \rightleftarrows^{k_3}_{k_4} N_2 \rightleftarrows^{k_5}_{k_6} N_0$ governed by ODEs $$ \frac{dN_0}{dt} = -k_1N_0 -k_6N_0 + k_2N_1 + k_5N_2 \\ \frac{dN_1}{dt} = k_1N_0 -k_2N_1 - k_3N_1 + k_4N_2 \\ \frac{dN_2}{dt} = k_6N_0 +k_3N_1 - k_4N_2 - k_5N_2 \\ $$ with mass constraint $N_0 + N_1 + N_2 = N_T$. From chemical reaction network theory I know that it can only have one equilibrium given $N_T > 0$. Alternatively, one can confirm that the rank of the matrix below is 3. However, deriving it analytically is a huge pain in the ass, therefore I was hoping to use symbolic algebra systems to solve this system of equations at equilibrium, namely: $$ \begin{equation} \begin{cases} 0 = -k_1N_0 -k_6N_0 + k_2N_1 + k_5N_2 \\ 0 = k_1N_0 -k_2N_1 - k_3N_1 + k_4N_2 \\ 0 = k_6N_0 +k_3N_1 - k_4N_2 - k_5N_2 \\ N_T = N_0 + N_1 + N_2 \end{cases} \end{equation} $$ which can be written as $Ax=b$ as follows: $$\left[\begin{array}{c, c, c} -k_1 - k_6 & k_2 & k_5 \\ k_1 & -k_2-k3 & k_4 \\ k_6 & k_3 & -k_4-k_5 \\ 1 & 1 & 1\\ \end{array}\right] \left[\begin{array}{c} N_0\\ N_1\\ N_2\\ \end{array}\right] = \left[\begin{array}{c} 0\\ 0\\ 0\\ N_T\\ \end{array}\right] $$ Now this can be solved algebraically by transforming it into a reduced row echelon form. However, a standard Gauss-Jordan elimination results in a matrix that says that the system is inconsistent, the extended matrix is: $$\left[\begin{array}{c, c, c, |c} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & N_T\\ \end{array}\right] $$ The working steps are omitted because there are a lot of them and they are lengthy. But the basic idea is that the last row simply gets eliminated and $N_T$ stays because all other elements are 0.

So, what is actually going on here? Is the basic Gauss-Jordan elimination (e.g. as given here) known to fail in such cases? Are there simple modifications to the algorithm that would take care of this? Any tips appreciated, if there are any questions/inaccuracies in my question please ask and I will try to clarify/expand.

Answer 1

No way of knowing what the computer is doing, especially as symbolic quantities are involved. Your first three rows sum to a zero row, meaning redundant (dependent). Drop the third row

$$\left[\begin{array}{c, c, c} -k_1 - k_6 & k_2 & k_5 \\ k_1 & -k_2-k_3 & k_4 \\ 1 & 1 & 1\\ \end{array}\right] \left[\begin{array}{c} N_0\\ N_1\\ N_2\\ \end{array}\right] = \left[\begin{array}{c} 0\\ 0\\ N_T\\ \end{array}\right] $$

This says that the N vector is any scalar multiple of the cross product of the first two rows. The particular scalar multiple is one that makes the sum to agree with $N_T$

For that matter, there is no particular difficulty finding the inverse of this 3 by 3 coefficient matrix.

Answer 2

@Will Jagy, posting as answer because I want to show the steps, as outlined here. Here goes (apologies for any arithmetical mistakes)

$$ \left[\begin{array}{c, c, c | c} -7 & 4 & 5 & 0\\ 1 & -8 & 9 & 0\\ 6 & 4 & -14 & 0\\ 1 & 1 & 1 & 73\\ \end{array}\right] $$ Scale row 1 by $M_{1,1}$ $$ \left[\begin{array}{c, c, c | c} 1 & \frac{-4}{7} & \frac{-5}{7} & 0\\ 1 & -8 & 9 & 0\\ 6 & 4 & -14 & 0\\ 1 & 1 & 1 & 73\\ \end{array}\right] $$ Subtract row 1 * $M_{n, 1}$ for $n = [2, 3, 4]$ from row $n$ $$ \left[\begin{array}{c, c, c | c} 1 & \frac{-4}{7} & \frac{-5}{7} & 0\\ 0 & \frac{-52}{7} & \frac{68}{7} & 0\\ 0 & \frac{44}{7} & \frac{-93}{7} & 0\\ 0 & \frac{11}{7} & \frac{12}{7} & 73\\ \end{array}\right] $$ Scale row 2 by $M_{2,2}$ $$ \left[\begin{array}{c, c, c | c} 1 & \frac{-4}{7} & \frac{-5}{7} & 0\\ 0 & 1 & \frac{-68}{52} & 0\\ 0 & \frac{44}{7} & \frac{-93}{7} & 0\\ 0 & \frac{11}{7} & \frac{12}{7} & 73\\ \end{array}\right] $$ Subtract row 2 * $M_{n, 2}$ for $n = [1, 3, 4]$ from row $n$ $$ \left[\begin{array}{c, c, c | c} 1 & 0 & \frac{-532}{364} & 0\\ 0 & 1 & \frac{-68}{52} & 0\\ 0 & 0 & \frac{-1844}{364} & 0\\ 0 & 0 & \frac{1372}{364} & 73\\ \end{array}\right] $$ Scale row 3 by $M_{3,3}$ $$ \left[\begin{array}{c, c, c | c} 1 & 0 & \frac{-532}{364} & 0\\ 0 & 1 & \frac{-68}{52} & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & \frac{1372}{364} & 73\\ \end{array}\right] $$ Subtract row 3 * $M_{n, 3}$ for $n = [1, 2, 4]$ from row $n$ $$ \left[\begin{array}{c, c, c | c} 1 & 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 73\\ \end{array}\right] $$ and so the matrix is in reduced row echelon form with an inconsistent final row. Where is the mistake in this "algorithm"?