The Kitchin Research Group

Given these equations, find [x1, x2, x3]

reference: Kreysig, Advanced Engineering Mathematics, 9th ed. Sec. 7.3

When solving linear equations, we can represent them in matrix form. The we simply use numpy.linalg.solve to get the solution.

import numpy as np
A = np.array([[1, -1, 1],
              [0, 10, 25],
              [20, 10, 0]])

b = np.array([0, 90, 80])

x = np.linalg.solve(A, b)
print x
print np.dot(A,x)

# Let us confirm the solution.
# this shows one element is not equal because of float tolerance
print np.dot(A,x) == b

# here we use a tolerance comparison to show the differences is less
# than a defined tolerance.
TOLERANCE = 1e-12
print np.abs((np.dot(A, x) - b)) <= TOLERANCE

[ 2.  4.  2.]
[  2.66453526e-15   9.00000000e+01   8.00000000e+01]
[False  True  True]
[ True  True  True]

It can be useful to confirm there should be a solution, e.g. that the equations are all independent. The matrix rank will tell us that. Note that numpy:rank does not give you the matrix rank, but rather the number of dimensions of the array. We compute the rank by computing the number of singular values of the matrix that are greater than zero, within a prescribed tolerance. We use the numpy.linalg.svd function for that. In Matlab you would use the rref command to see if there are any rows that are all zero, but this command does not exist in numpy. That command does not have practical use in numerical linear algebra and has not been implemented.

import numpy as np
A = np.array([[1, -1, 1],
              [0, 10, 25],
              [20, 10, 0]])

b = np.array([0, 90, 80])

# determine number of independent rows in A we get the singular values
# and count the number greater than 0.
TOLERANCE = 1e-12
u, s, v = np.linalg.svd(A)
print 'Singular values: {0}'.format(s)
print '# of independent rows: {0}'.format(np.sum(np.abs(s) > TOLERANCE))

# to illustrate a case where there are only 2 independent rows
# consider this case where row3 = 2*row2.
A = np.array([[1, -1, 1],
              [0, 10, 25],
              [0, 20, 50]])

u, s, v = np.linalg.svd(A)

print 'Singular values: {0}'.format(s)
print '# of independent rows: {0}'.format(np.sum(np.abs(s) > TOLERANCE))

Singular values: [ 27.63016717  21.49453733   1.5996022 ]
# of independent rows: 3
Singular values: [ 60.21055203   1.63994657  -0.        ]
# of independent rows: 2

Matlab comparison

Copyright (C) 2013 by John Kitchin. See the License for information about copying.

org-mode source