Matrices#

Mathematics Methods 1

Defining a matrix#

Slide 3

Remember: when stating the dimension of a matrix (e.g. \(3 \times 3\), \(2 \times 4\), etc) it’s Rows FiRst, Columns SeCond.

Define matrices using numpy.array, and find their dimensions using numpy.shape.

# import necessary modules
import numpy as np

# mA is a SQUARE matrix
mA = np.array([[2, 3, -4],
                [3, -1, 2],
                [4, 2, 2]])

# mB is a NON-SQUARE matrix
mB = np.array([[2, -4, 1, -2],
               [7, 8, 0, 3]])

# So, mA is 3 x 3 matrix, mB is a 2 x 4 matrix.

print("The dimension of A is", np.shape(mA))
print("The dimension of B is", np.shape(mB))

The dimension of A is (3, 3)
The dimension of B is (2, 4)

Matrix-scalar multiplication#

Slide 5

mD = np.matrix([[2, -4, 1, -2],
                [7, 8, 0, 3]])

print("3*D =")
print(3*mD)

3*D =
[[  6 -12   3  -6]
 [ 21  24   0   9]]

Matrix addition and subtraction#

Slide 6, 7

mE = np.array([[2, -4],
               [7, 8]])

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

print("E + F =")
print(mE + mF)

print("E - F =")
print(mE - mF)

E + F =
[[ 6 -4]
 [ 9  7]]
E - F =
[[-2 -4]
 [ 5  9]]

The two matrices being added together must have the same dimension. If they do not, Python will give ValueError.

Adding mE ( \(2 \times 2\) ) and mB ( \(2 \times 4\) ):

print(mE + mB)

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-5f7c0614b2fd> in <module>
----> 1 print(mE + mB)

ValueError: operands could not be broadcast together with shapes (2,2) (2,4)

Matrix-matrix multiplication#

Slide 10

Matrix-matrix multiplication of two numpy arrays is performed with the @ operator.

mE = np.array([[2, -4],
               [7, 8],
               [3, 2]])
mH = np.array([[4, 0, 0, 1],
               [2, -1, 3, -2]])

print("E @ H =")
print(mE @ mH)

Matrix transpose#

Slide 12, 13

Finding the transpose of a matrix using numpy.transpose.

# Non-square matrix
mB = np.array([[2, -4],
               [7, 8],
               [3, 2]])

mB_transpose = np.transpose(mB)
print("The transpose of (non-square) B =")
print(mB_transpose)

# Square matrix
mB = np.array([[2, 3, -4],
               [3, -1, 2],
               [4, 2, 2]])

mB_transpose = np.transpose(mB)
print("The transpose of (square) B =")
print(mB_transpose)

Identity matrix#

Slide 15

Create identity matrix using numpy.identity.

mI = np.identity(3)

print("I=")
print(mI)

mB = np.array([[2, 3, -4],
               [3, -1, 2],
               [4, 2, 2]])

mIB = mI @ mB
print("I @ B =")
print(mIB)

# mIB should be equal to mB
if np.allclose(mB, mIB) == True:
    print("B = I @ B")
else:
    print("Error")

Determinants#

Slide 16

Finding determinants of matrices using numpy.linalg.deg

# 2 x 2 example
mD = np.array([[1, 4],
               [-3, 5]])
sDet = np.linalg.det(mD)
print("det(D) = %.2f" % sDet) 

# 3 x 3 example
mA = np.array([[1, 2, 4],
               [-3, 1, -2],
               [3, 2, 5]])
sDet = np.linalg.det(mA)
print("det(A) = %.2f" % sDet)

Solving systems of linear equations#

Find solutions of a system of linear equations of the form \( A \mathbf{x} = \mathbf{b} \) using numpy.linalg.solve(A, b).

# Set up a system of the form Ax = b
mA = np.array([[2, 3, -4],
               [3, -1, 2],
               [4, 2, 2]])
vb = np.array([10, 3, 8])

# Solve the system of linear equations
vX = np.linalg.solve(mA, vb)

print("The solution to Ax = b is:")
print(vX)

Existence of solutions#

Is the solution above unique?

# Compute the determinant of the square matrix A
sDet = np.linalg.det(mA)
print("Determinant of A = %.3f" % sDet) 
if(sDet != 0):
    print("The system Ax = b has a UNIQUE solution.") 
else:
    print("The system Ax = b has either infinite solutions or no solutions.")

Inverse of a matrix#

Compute the inverse of a matrix using numpy.linalg.inv

mA_inv = np.linalg.inv(mA)
print("The inverse of A is:") 
print(mA_inv)

Note

A is an object with a set of methods (you’ll learn more about these in the Programming for Geophysicists course). Therefore, you could also use the method getI() to compute the inverse by typing

A.getI()

Let us check that A @ A_inv gives the identity matrix I. Note: You may find that some elements of the matrix are around 1e-16 (i.e. \(10^{-16}\) ) instead of exactly zero.

print("A @ A_inv=") 
print(mA @ mA_inv) 

mI = np.identity(3)

if np.allclose(mI, mA @ mA_inv) == True:
    print("A*A_inv gives the identity matrix")
else:
    print("Error")