The Kitchin Research Group

#Problem set 7.4 #17
import numpy as np

v1 = [0.2, 1.2, 5.3, 2.8, 1.6]
v2 = [4.3, 3.4, 0.9, 2.0, -4.3]

A = np.row_stack([v1, v2])
U, s, V = np.linalg.svd(A)
print s

[ 7.57773162  5.99149259]

You can tell by inspection the rank is 2 because there are no near-zero singular values.

the rank command roughly works in the following way: the matrix is converted to a reduced row echelon form, and then the number of rows that are not all equal to zero are counted. Matlab uses a tolerance to determine what is equal to zero. If there is uncertainty in the numbers, you may have to define what zero is, e.g. if the absolute value of a number is less than 1e-5, you may consider that close enough to be zero. The default tolerance is usually very small, of order 1e-15. If we believe that any number less than 1e-5 is practically equivalent to zero, we can use that information to compute the rank like this.

import numpy as np

A = [[1, 2, 3],
     [0, 2, 3],
     [0, 0, 1e-6]]

U, s, V = np.linalg.svd(A)
print s
print np.sum(np.abs(s) > 1e-15)
print np.sum(np.abs(s) > 1e-5)

[  5.14874857e+00   7.00277208e-01   5.54700196e-07]
3
2

reference: Exercise 2.4 in Chemical Reactor Analysis and Design Fundamentals by Rawlings and Ekerdt.

The following reactions are proposed in the hydrogenation of bromine:

Let this be our species vector: v = [H2 H Br2 Br HBr].T

the reactions are then defined by M*v where M is a stoichometric matrix in which each row represents a reaction with negative stoichiometric coefficients for reactants, and positive stoichiometric coefficients for products. A stoichiometric coefficient of 0 is used for species not participating in the reaction.

import numpy as np

#    [H2  H Br2 Br HBr]
M = [[-1,  0, -1,  0,  2],  # H2 + Br2 == 2HBR
     [ 0,  0, -1,  2,  0],  # Br2 == 2Br
     [-1,  1,  0, -1,  1],  # Br + H2 == HBr + H
     [ 0, -1, -1,  1,  1],  # H + Br2 == HBr + Br
     [ 1, -1,  0,  1,  -1], # H + HBr == H2 + Br
     [ 0,  0,  1, -2,  0]]  # 2Br == Br2

U, s, V = np.linalg.svd(M)
print s
print np.sum(np.abs(s) > 1e-15)

import sympy 
M = sympy.Matrix(M)
reduced_form, inds = M.rref()

print reduced_form

labels = ['H2',  'H', 'Br2', 'Br', 'HBr']
for row in reduced_form.tolist():
    s = '0 = '
    for nu,species in zip(row,labels):
        if nu != 0:
            
            s += ' {0:+d}{1}'.format(int(nu), species)
    if s != '0 = ': print s

[  3.84742803e+00   3.32555975e+00   1.46217301e+00   1.73313660e-16
   8.57422679e-17]
3
[1, 0, 0,  2, -2]
[0, 1, 0,  1, -1]
[0, 0, 1, -2,  0]
[0, 0, 0,  0,  0]
[0, 0, 0,  0,  0]
[0, 0, 0,  0,  0]
0 =  +1H2 +2Br -2HBr
0 =  +1H +1Br -1HBr
0 =  +1Br2 -2Br

6 reactions are given, but the rank of the matrix is only 3. so there are only three independent reactions. You can see that reaction 6 is just the opposite of reaction 2, so it is clearly not independent. Also, reactions 3 and 5 are just the reverse of each other, so one of them can also be eliminated. finally, reaction 4 is equal to reaction 1 minus reaction 3.

There are many possible independent reactions. In the code above, we use sympy to put the matrix into reduced row echelon form, which enables us to identify three independent reactions, and shows that three rows are all zero, i.e. they are not independent of the other three reactions. The choice of independent reactions is not unique.