Matrix Manipulation

In [1]:
import numpy as np
  • We can compute the transpose using $\textbf{np.tranpose}$
  • We can also compute the transpose of a matrix $\textbf{A}$ using $\textbf{A.T}$
  • The function $\textbf{np.allclose}$ can be used to find if the values are equal or not
    • We can specify the tolerance using $\textbf{rtol}$ while comparing floating point values
In [2]:
A = np.array([[3, 2, 1], [5, 7, 4], [9, 6, 8]]); print(A);
B = np.transpose(A); print(B);
C = A.T; print(C);
print(np.allclose(B, C)); # To check whether all the elements of B are equal to C. 

[[3 2 1]
 [5 7 4]
 [9 6 8]]
[[3 5 9]
 [2 7 6]
 [1 4 8]]
[[3 5 9]
 [2 7 6]
 [1 4 8]]
True
  • We can compute the determinant using $\textbf{numpy.linalg.det}$
In [3]:
print(A);
d = np.linalg.det(A);
print(d);

[[3 2 1]
 [5 7 4]
 [9 6 8]]
55.000000000000014
  • Let's see whether $|\mathbf{A}\cdot \mathbf{B}| = |\mathbf{A}|\cdot |\mathbf{B}|$
  • We generate two matrix $\textbf{A}$ and $\textbf{B}$ with random integer values
  • We perform matrix multiplication and store it in $\textbf{C}$
In [4]:
A = np.random.randint(1, 10, size=(3,3));
B = np.random.randint(1, 10, size=(3,3));

C = np.dot(A, B); # Matrix multiplication

d = np.linalg.det(C);
e = np.linalg.det(A)*np.linalg.det(B);

print("LHS:%3.16f\n"%d);
print("RHS:%3.16f\b"%e);

print(np.isclose(d, e, rtol=1e-5));

print("Relerr:%3.16f"%(np.abs(1-d/e)));

LHS:3773.9999999999931788

RHS:3774.000000000001364
True
Relerr:0.0000000000000022

Finding out the inverse

In [5]:
print(A);

[[1 7 5]
 [3 5 9]
 [2 3 8]]
  • We check whether determinant is non-zero
In [6]:
a = np.linalg.det(A); print(a);

-34.00000000000002
  • We can compute the inverse of a matrix using $\textbf{np.linalg.inv}$
In [7]:
B = np.linalg.inv(A); print(B);

[[-0.38235294  1.20588235 -1.11764706]
 [ 0.17647059  0.05882353 -0.17647059]
 [ 0.02941176 -0.32352941  0.47058824]]
  • Let's check whether $A\cdot B = I$
    • We round the values to 1 decimal place
In [8]:
C = np.dot(A, B); print(np.round(C,1));

[[ 1.  0. -0.]
 [-0.  1. -0.]
 [-0.  0.  1.]]
  • Let's check whether $(A^{-1})^T = (A^T)^{-1}$
In [9]:
lhs = np.transpose(np.linalg.inv(A));
rhs = np.linalg.inv(np.transpose(A));
print(lhs); 
print(rhs);
print(np.allclose(lhs, rhs))

[[-0.38235294  0.17647059  0.02941176]
 [ 1.20588235  0.05882353 -0.32352941]
 [-1.11764706 -0.17647059  0.47058824]]
[[-0.38235294  0.17647059  0.02941176]
 [ 1.20588235  0.05882353 -0.32352941]
 [-1.11764706 -0.17647059  0.47058824]]
True
  • Multiplication of matrices is not commutative: $A\cdot B \neq B\cdot A$
In [10]:
A = np.random.randint(1, 10, size=(3,3));
B = np.random.randint(1, 10, size=(3,3));
C = np.dot(A,B); D = np.dot(B, A);
print(C); print(D);
print(np.allclose(C, D));

[[113  66 142]
 [ 79  26  76]
 [ 35  18  40]]
[[ 32  84  33]
 [ 65 111  56]
 [ 38  90  36]]
False

A function to check whether a matrix is symmetric or not

  • A symmteric matrix has the following property: $A_{ij} = A_{ji}$
In [11]:
def check_symm(aditya):
    return(np.allclose(aditya, aditya.T));
In [12]:
A = np.array([[1, 2, 3], [2, 4, 5], [3, 5, 6]]); print(A);

[[1 2 3]
 [2 4 5]
 [3 5 6]]
In [13]:
print(check_symm(A));

True
  • This is useful for debugging purposes. Suppose you are working with stresses then you can check whether the matrices are symmetric or we have made mistakes somewhere

Matrix-Matrix double dot product

  • The double dot product is defined as: $A:B = A_{ij}B_{ij} = \sum_{i, j} A_{ij}B_{ij}$
  • Example viscous dissipation term: $\tau_{ij}S_{ij}$ where $\tau$ is the stress tensor and $S$ is the rate of deformation tensor
In [14]:
print(A); print(B);

[[1 2 3]
 [2 4 5]
 [3 5 6]]
[[3 4 7]
 [8 2 7]
 [4 4 6]]
  • We can compute the double dot product of $\textbf{A}$ and $\textbf{B}$ using the function $\textbf{np.tensordot}$
In [15]:
c = np.tensordot(A, B); print(c);

159
  • We can verify this is correct by using nested for loops to compute the sum again
In [16]:
s = 0;
for i in np.arange(0,3):
    for j in np.arange(0,3):
        s = s + A[i, j]*B[i, j];
print(s)

159
  • Another way of computing the tensor product is by reshaping the matrices to 1D vectors using $\textbf{np.ndarray.flatten}$ and then taking element-wise product and summing them using $\textbf{np.sum}$.
  • We could have also used $\textbf{np.dot}$ after reshaping the matrices
In [17]:
# A vectorized way of computing the tensordot product
C = np.sum(np.ndarray.flatten(A)*np.ndarray.flatten(B))
print(C)

159

Transformation matrix

  • $\mathbf{v}'= \mathbf{Q}\cdot \mathbf{v}$, here $\textbf{Q}$ is the transformation matrix

  • We specify the angle of rotation as $\theta$ and then declare the transformation matrix $\textbf{Q}$ based on that

In [18]:
theta = 30*np.pi/180;
Q = np.array([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]]);
print(Q)

[[ 0.8660254  0.5      ]
 [-0.5        0.8660254]]
  • Some properties of the tranformation matrix:
    • $\det(\mathbf{Q}) = 1$
In [19]:
c = np.linalg.det(Q); print(c);

1.0
  • $\mathbf{Q^T} = \mathbf{Q^{-1}}$
In [20]:
lhs = np.transpose(Q);
rhs = np.linalg.inv(Q);
print(np.allclose(lhs, rhs, rtol=1e-6))

True
  • Identity matrix can be generated using $\textbf{np.eye}$
In [21]:
I = np.eye(2); print(I);

[[1. 0.]
 [0. 1.]]
  • $\mathbf{Q\cdot Q^T} = \mathbf{I}$
In [22]:
print(np.allclose(np.dot(Q, Q.T), I));

True
  • Now we can tranform a vector in original frame to the equivalent vector in rotated frame.
  • If $\textbf{a} = (1,3)$, let's find out $\textbf{a'}$
In [23]:
a = np.array([1, 3])
ap = np.dot(Q, a);
print(ap)

[2.3660254  2.09807621]
  • We now generate a series of transformation matrices by looping over the angle $\theta$ and look at how the elements of $\textbf{a'}$: $\textrm{a'}_1$ and $\textrm{a'}_2$ change
In [24]:
# vary Theta from 0 to 2pi radians
theta_a = np.linspace(0, 2*np.pi, 100); # Theta in radians

# Elements of vector in initial frame: theta = 0
a = np.array([1,3]);

# Generate a matrix to store the transformed elements for correspoding theta
ap1_a = np.zeros(np.shape(theta_a));
ap2_a = np.zeros(np.shape(theta_a));

# Iterate over theta values and assign the values to ap1_a and ap2_a
count = 0;
for theta in theta_a:
    Q = np.array([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]]);
    ap = np.dot(Q, a);
    ap1_a[count] = ap[0]; # 1st element, index = 0
    ap2_a[count] = ap[1]; # 2nd element, index = 1
    count = count + 1;
  • We can now plot the elements as function of theta
In [25]:
import matplotlib.pyplot as plt
plt.rcParams.update({"text.usetex":True})
%config InlineBackend.figure_format='svg';
plt.plot(theta_a, ap1_a, label="$a'_1$");
plt.plot(theta_a, ap2_a, label="$a'_2$");
plt.grid(True);
plt.xlabel("$\\theta$");
plt.ylabel("Components of a'");
plt.legend();
  • What we have here is a parametrized equation for $\textrm{a'}_1$ and $\textrm{a'}_2$ in terms of $\theta$
  • We can plot $\textrm{a'}_1$ vs $\textrm{a'}_2$ and see what the curve looks like
In [26]:
plt.plot(ap1_a, ap2_a);
ax = plt.gca();
ax.set_aspect(1);
plt.plot(a[0], a[1], 'ok'); # initial vector

Transformation of Stress Tensor

  • The stress transformation is give by: $ \sigma' = \textrm{Q} \cdot \sigma \cdot \textrm{Q}^T $
  • Let's take an example of 2D stress
In [27]:
theta = 30*np.pi/180;
Q = np.array([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]]);

# Original state of stress
s = np.array([[50, 30], [30, -20]]);
print(s);

# Transformed matrix
sp = np.dot(np.dot(Q,s), Q.T);
print(sp);

[[ 50  30]
 [ 30 -20]]
[[ 58.48076211 -15.31088913]
 [-15.31088913 -28.48076211]]
  • Let's loop over $ \theta $ and see how state of stress changes
In [28]:
theta_a = np.linspace(0, np.pi, 100); # Theta in radians

sigma_n_a = np.zeros(np.shape(theta_a));
tau_n_a = np.zeros(np.shape(theta_a));

count =0;
for theta in theta_a:
    Q = np.array([[np.cos(theta), np.sin(theta)], [-np.sin(theta), np.cos(theta)]]);
    s = np.array([[50, 30], [30, -20]]);
    sp = np.dot(np.dot(Q,s), Q.T);
    sigma_n_a[count] = sp[0,0];
    tau_n_a[count] = sp[0,1];
    count+=1; # count = count + 1
  • Let's plot the Mohr's circle
In [29]:
plt.plot(sigma_n_a, tau_n_a);
plt.plot(np.trace(s)/2, 0, 'ok'); # center of the circle
plt.xlabel("$\\sigma_n'$");
plt.ylabel("$\\tau_n'$");

ax = plt.gca();
ax.set_aspect(1);
In [ ]: