** Solving linear equations
:PROPERTIES:
:categories: Linear algebra
:date: 2013/02/27 13:13:06
:updated: 2013/02/27 13:13:06
:END:
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.
#+BEGIN_SRC python
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
#+END_SRC
#+RESULTS:
: [ 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.
#+BEGIN_SRC python
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))
#+END_SRC
#+RESULTS:
: Singular values: [ 27.63016717 21.49453733 1.5996022 ]
: # of independent rows: 3
: Singular values: [ 60.21055203 1.63994657 -0. ]
: # of independent rows: 2
[[http://matlab.cheme.cmu.edu/2011/08/01/solving-linear-equations/][Matlab comparison]]