ESE - Numerical Methods I: Systems of Linear Equations



In [1]:

    
%pylab inline
import numpy as np

# scipy linear algebra module
import scipy.linalg as scal









    



Populating the interactive namespace from numpy and matplotlib

Problems

1.)

By evaluating the determinant, classify the following matrices as singular, ill conditioned, or well conditioned (using `scal.det`)

\begin{equation} A_{1} = \begin{pmatrix} 1 & 2 & 3 \\ 2 & 3 & 4 \\ 3 & 4 & 5 \end{pmatrix} \end{equation}\begin{equation} A_{2} = \begin{pmatrix} 2.11 & -0.80 & 1.77 \\ -1.84 & 3.03 & 1.21 \\ -1.57 & 5.25 & 4.30 \end{pmatrix} \end{equation}\begin{equation} A_{3} = \begin{pmatrix} 2 & -1 & 0 \\ -1 & 2 & -1 \\ 0 & -1 & 2 \end{pmatrix} \end{equation}\begin{equation} A_{4} = \begin{pmatrix} 4 & 3 & -1 \\ 7 & -2 & 3 \\ 5 & -18 & 13 \end{pmatrix} \end{equation}

NB: You can use `scal.norm` and the ratio of determinant to norm for this.



In [2]:

    
A1 =  np.array([[1,2,3],[2,3,4],[3,4,5]])
A2 =  np.array([[2.11,-0.80,1.77],[-1.84,3.03,1.21],[-1.57,5.25,4.3]])
A3 =  np.array([[2,-1,0],[-1,2,-1,],[0,-1,2]])
A4 =  np.array([[4,3,-1],[7,-2,3],[5,-18,13]])

def conditioning(A):
    d = scal.det(A)  # determinant
    n = scal.norm(A) # frobenius matrix norm
    cn = d/n
    if cn == 0.:
        cs = 'singular'
    elif cn < 0.1:
        cs = 'ill conditioned'
    else:
        cs = 'well conditioned'
    print 'determinant is {0:.2e}, matrix norm is {1:.2e} and matrix is'.format(d,n), cs

conditioning(A1)
conditioning(A2)
conditioning(A3)
conditioning(A4)









    



determinant is 2.22e-16, matrix norm is 9.64e+00 and matrix is ill conditioned
determinant is 5.99e-01, matrix norm is 8.41e+00 and matrix is ill conditioned
determinant is 4.00e+00, matrix norm is 4.00e+00 and matrix is well conditioned
determinant is 0.00e+00, matrix norm is 2.46e+01 and matrix is singular

2.)

Solve the linear system $A x = b$ for x (using `scal.solve`).

\begin{equation} A = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix} \end{equation}\begin{equation} b = \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix} \end{equation}

check the solution by assuring that np.dot(A, x)-b is numerically equal to 0., (using `np.allclose`)



In [3]:

    
A = np.array([[1,2,3], [4,5,6], [7,8,9]])
b = np.array([1,2,3])

x = scal.solve(A, b)
print x
np.allclose(np.dot(A, x)-b, 0.)









    



[-0.33333333  0.66666667  0.        ]






    Out[3]:





True

3.)

Find $L$ and $U$ for

\begin{equation} A = \begin{pmatrix} 4 & -1 & 0 \\ -1 & 4 & -1 \\ 0 & -1 & 4 \end{pmatrix} \end{equation}

(using `lu`)



In [4]:

    
A = np.array([[1,-1,0], [-1,4,-1], [0,-1,4]])

P,L,U = scal.lu(A)
np.dot(L,U) # matrix multiplication









    Out[4]:





array([[ 1., -1.,  0.],
       [-1.,  4., -1.],
       [ 0., -1.,  4.]])

4.)

Solve the linear system $A x = b$ for x (using the LU factorization `scal.lu_factor` and LU solver `scal.lu_solve`).

\begin{equation} A = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix} \end{equation}\begin{equation} b = \begin{pmatrix} 1 \\ 2 \\ 3 \end{pmatrix} \end{equation}

check the solution by assuring that np.dot(A, x)-b is numerically equal to 0., (using `np.allclose`)



In [5]:

    
A = np.array([[1,2,3], [4,5,6], [7,8,9]])
b = np.array([1,2,3])

LU, P = scal.lu_factor(A)
print LU
x = scal.lu_solve((LU,P), b)
print x
np.allclose(np.dot(A, x)-b, 0.)









    



[[  7.00000000e+00   8.00000000e+00   9.00000000e+00]
 [  1.42857143e-01   8.57142857e-01   1.71428571e+00]
 [  5.71428571e-01   5.00000000e-01   1.11022302e-16]]
[-0.33333333  0.66666667  0.        ]






    Out[5]:





True

5.) We now know the solution to

\begin{equation} A = \begin{pmatrix} 2 & 3 & -1 \\ 3 & 2 & -5 \\ 1 & 4 & -1 \end{pmatrix} \end{equation}\begin{equation} b = \begin{pmatrix} 1 \\ -1 \\ 2 \end{pmatrix} \end{equation}

\begin{equation} x = \begin{pmatrix} 1.84375 \\ 0.40652 \\ 1.46875 \end{pmatrix} \end{equation}

Without using functions from the linalg module (but e.g. np.dot), complete the code below so that the result is True.



In [6]:

    
A = np.array([[2,-3,-1], [3,2,-5], [1,4,-1]],dtype=np.float) 
b = np.array([1,-1,2],dtype=np.float)

def lu_dec(A):
    """
    Returns the matrix A where [A] = [L\U], i.e., [U] in the upper triangle and 
    the nondiagonal terms of [L] in the lower triangle.
    """
    n = len(A)
    for k in range(0,n-1):
        for i in range(k+1,n):
            if A[i,k] != 0.0:
                lam = A[i,k]/A[k,k]
                A[i,k+1:n] = A[i,k+1:n] - lam*A[k,k+1:n]
                A[i,k] = lam
    return A

def lu_sol(LU,b):
    """
    Returns the solution vector x as from [L][U]{x} = b, where [a] = [L\U] is the matrix
    returned from lu_dec.
    """
    n = len(LU)
    for k in range(1,n):
        b[k] = b[k] - dot(LU[k,0:k],b[0:k])
    for k in range(n-1,-1,-1):
        b[k] = (b[k] - dot(LU[k,k+1:n],b[k+1:n]))/LU[k,k]
    return b

LU = lu_dec(A)
x = lu_sol(LU, b)
np.allclose(x, [1.84375,0.40625,1.46875])









    Out[6]:





True



In [6]: