Solving linear equations
Posted February 27, 2013 at 01:13 PM | categories: linear algebra | tags:
Given these equations, find [x1, x2, x3]
\begin{eqnarray} x_1 - x_2 + x_3 &=& 0 \\ 10 x_2 + 25 x_3 &=& 90 \\ 20 x_1 + 10 x_2 &=& 80 \end{eqnarray}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
Copyright (C) 2013 by John Kitchin. See the License for information about copying.