Solving Linear Systems

Teng-Jui Lin

Content adapted from UW CHEME 375, Chemical Engineering Computer Skills, in Spring 2021.

• Python skills and numerical methods

  • Matrix inversion by numpy.linalg.inv()

  • Solving linear systems by numpy.linalg.solve()

• ChemE applications

  • Balancing chemical equations

Balancing chemical equations

Problem Statement. (FRB 4.62a) Mammalian cells convert glucose to glutamic acid by the reaction

$${\mathrm{C}}_6{\mathrm{H}}_{1}2{\mathrm{O}}_6+a{\mathrm{NH}}_3+b{\mathrm{O}}_2\to p{\mathrm{C}}_5{\mathrm{H}}_9{\mathrm{NO}}_4+q{\mathrm{CO}}_2+r{\mathrm{H}}_2\mathrm{O}$$

Determine the stoichiometric coefficients.

Solution. We establish atomic balance between reactants and products:

$$\begin{array}{r}\begin{array}{rlrl}& \text{C balance:}& 6& =5p+q\\ & \text{H balance:}& 12+3a& =9p+2r\\ & \text{O balance:}& 6+2b& =4p+2q+r\\ & \text{N balance:}& a& =p\end{array}\end{array}$$

Reorganize the equations so that variables are on the left and constant are on the right:

$$\begin{array}{r}\begin{array}{rlrl}& \text{C balance:}& 5p+q& =6\\ & \text{H balance:}& 9p+2r-3a& =12\\ & \text{O balance:}& 4p+2q+r-2b& =6\\ & \text{N balance:}& a-p& =0\end{array}\end{array}$$

Because we have five variables but only four equations, we need to assume an arbitrary basis. Here, we can assume a basis of stoichiometric coefficient of ${\mathrm{C}}_5{\mathrm{H}}_9{\mathrm{NO}}_4$ being one:

$$\begin{array}{rlrl}& \text{Basis:}& p& =1\end{array}$$

Now we have five variables and equations, so we can rewrite them to the form of

$$\mathbf{A}\mathbf{x}=\mathbf{b}$$

where

$$\begin{array}{r}\mathbf{A}=\left[\begin{array}{ccccc}0& 0& 5& 1& 0\\ -3& 0& 9& 0& 2\\ 0& -2& 4& 2& 1\\ 1& 0& -1& 0& 0\\ 0& 0& 1& 0& 0\end{array}\right],\mathbf{b}=\left[\begin{array}{c}6\\ 12\\ 6\\ 0\\ 1\end{array}\right],\mathbf{x}=\left[\begin{array}{c}a\\ b\\ p\\ q\\ r\end{array}\right]\end{array}$$

We can then solve for $\mathbf{x}$ and get the coefficients:

$$\begin{array}{r}\begin{array}{rl}{\mathbf{A}}^{-1}\mathbf{A}\mathbf{x}& ={\mathbf{A}}^{-1}\mathbf{b}\\ \mathbf{x}& ={\mathbf{A}}^{-1}\mathbf{b}\end{array}\end{array}$$

The resulting coefficients are scalable, so we can scale them to whole number for reporting if needed.

Implementation: solving linear systems using numpy.linalg.inv()

In this approach, we use numpy.linalg.inv() to solve the linear system by matrix inversion.

import numpy as np
from numpy.linalg import inv

# define the matrix and vector
A = np.array([[0, 0, 5, 1, 0], 
               [-3, 0, 9, 0, 2], 
               [0, -2, 4, 2, 1], 
               [1, 0, -1, 0, 0],
               [0, 0, 1, 0, 0]])
b = np.array([6, 12, 6, 0, 1])

# calculate the unknown vector using inverse and matrix multiplication
x = inv(A) @ b
x

array([1. , 1.5, 1. , 1. , 3. ])

Therefore, we have the solution

$$\begin{array}{r}\mathbf{x}=\left[\begin{array}{c}a\\ b\\ p\\ q\\ r\end{array}\right]=\left[\begin{array}{c}1\\ 1.5\\ 1\\ 1\\ 3\end{array}\right].\end{array}$$

We have the balanced chemical equation

$${\mathrm{C}}_6{\mathrm{H}}_{12}{\mathrm{O}}_6+{\mathrm{NH}}_3+1.5{\mathrm{O}}_2\to {\mathrm{C}}_5{\mathrm{H}}_9{\mathrm{NO}}_4+{\mathrm{CO}}_2+3{\mathrm{H}}_2\mathrm{O}$$

Implementation: solving linear systems using numpy.linalg.solve()

In this approach, we use numpy.linalg.solve() to directly solve the linear system.

from numpy.linalg import solve

x = solve(A, b)
x

array([1. , 1.5, 1. , 1. , 3. ])

This approach gives the same result and is computationally more efficient.