Matrix Manipulation

  • We can compute the transpose of a matrix $\textbf{A}$ using $\textbf{A'}$
In [1]:
A= [3, 2, 1; 5, 7, 4; 9, 6, 8];
disp(A);

   3   2   1
   5   7   4
   9   6   8
In [2]:
B=A';
disp(B)

   3   5   9
   2   7   6
   1   4   8
In [3]:
disp(A);
d = det(A);
disp(d)

   3   2   1
   5   7   4
   9   6   8
 55.000
  • 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 = randi([0, 10], [3,3]);
disp('The matrix A is:')
disp(A)
B = randi([0, 10], [3,3]);
disp('The matrix B is:')
disp(B)
C=A*B;
d=det(C);
e=det(A)*det(B);
fprintf("LHS: %f\n", d);
fprintf("RHS: %f\n", e);
fprintf("Relerr: %f", abs(1-d/e));

The matrix A is:
   7   1   1
   3   7   3
   6   5   5
The matrix B is:
   9   1   5
   5   8   6
   0   8   1
LHS: -19140.000000
RHS: -19140.000000
Relerr: 0.000000

Finding out the Inverse of a Matrix

In [5]:
disp(A)

   7   1   1
   3   7   3
   6   5   5
  • We check whether determinant is non-zero
In [6]:
a = det(A); disp(a);

 116.00
  • We can compute the inverse of a matrix using $\textbf{inv(A)}$
In [7]:
B = inv(A); disp(B)

   1.7241e-01   1.0408e-17  -3.4483e-02
   2.5862e-02   2.5000e-01  -1.5517e-01
  -2.3276e-01  -2.5000e-01   3.9655e-01
  • Let's check whether $A\cdot B = I$
In [8]:
C=A*B; disp(round(C))

   1   0  -0
   0   1  -0
  -0  -0   1
  • Let's check whether $(A^{-1})^T = (A^T)^{-1}$
In [9]:
format short g
LHS= (inv(A))';
LHS = round(10000 * LHS) / 10000;
RHS= inv(A');
RHS = round(10000 * RHS) / 10000;
disp(LHS)
disp(RHS)
isequal(LHS,RHS)
compare= (LHS==RHS)

  0.1724  0.0259  -0.2328
  0  0.25  -0.25
  -0.0345  -0.1552  0.3966
  0.1724  0.0259  -0.2328
  0  0.25  -0.25
  -0.0345  -0.1552  0.3966
ans = 1
compare =

  1  1  1
  1  1  1
  1  1  1
  • Multiplication of matrices is not commutative: $A\cdot B \neq B\cdot A$
In [10]:
A = randi([0, 10], [3,3]);
B = randi([0, 10], [3,3]);
C=A*B; 
disp('The matrix C is:')
disp(C)
D=B*A;
disp('The matrix D is:')
disp(D)
compare= (C==D)

The matrix C is:
  57  31  55
  45  69  87
  111  73  101
The matrix D is:
  69  63  57
  107  123  51
  59  45  35
compare =

  0  0  0
  0  0  0
  0  0  0

A function to check whether a matrix is symmetric or not

  • A symmteric matrix has the following property: $A_{ij} = A_{ji}$
In [11]:
function check_symm(Matrix)
    if (isequal(Matrix, Matrix')==1);
    fprintf("This is a symmetric matrix");
    else
    fprintf("This is not a symmetric matrix");
    end
    end
In [12]:
A= [1, 2, 3; 2, 4, 5; 3, 5, 6];
B= [3, 3, 9; 9, 4, 9; 8, 8, 3];
disp("The matrix A is")
disp(A);
disp("The matrix B is")
disp(B);

The matrix A is
  1  2  3
  2  4  5
  3  5  6
The matrix B is
  3  3  9
  9  4  9
  8  8  3
In [13]:
check_symm(A)

This is a symmetric matrix
  • 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
In [14]:
check_symm(B)

This is not a symmetric matrix

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 [15]:
disp("The matrix A is")
disp(A);
disp("The matrix B is")
disp(B);

The matrix A is
  1  2  3
  2  4  5
  3  5  6
The matrix B is
  3  3  9
  9  4  9
  8  8  3
  • We can compute the double dot product of $\textbf{A}$ and $\textbf{B}$ by taking the element-wise dot product and summing the terms
In [16]:
C = sum(dot(A,B))

C = 197
  • We can verify this is correct by using nested for loops to compute the sum again
In [17]:
function double_dot(A,B)
   S=0;
    for i=1:1:3
        for j=1:1:3
            S = S + A(i,j)*B(i,j);
        end
    end
    fprintf("The double-dot product is: %d\t", S);
end
In [18]:
double_dot(A,B)

The double-dot product is: 197

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 [19]:
theta = 30*pi/180;
Q = [cos(theta), sin(theta); -sin(theta), cos(theta)];
disp(Q)

  0.86603  0.5
  -0.5  0.86603
  • Some properties of the tranformation matrix:
    • $\det(\mathbf{Q}) = 1$
In [20]:
c = det(Q); disp(c);

1
  • $\mathbf{Q^T} = \mathbf{Q^{-1}}$
In [21]:
lhs = Q';
rhs = inv(Q);
if (isequal(LHS,RHS)==1)
fprintf("TRUE \n");
end
compare= (LHS==RHS)

TRUE 
compare =

  1  1  1
  1  1  1
  1  1  1
  • Identity matrix can be generated using $\textbf{eye}$
In [22]:
I = eye(2); disp(I);

Diagonal Matrix

  1  0
  0  1
  • $\mathbf{Q\cdot Q^T} = \mathbf{I}$
In [23]:
compare (Q*Q'==I)

ans =

  1
  1
  • 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 [24]:
a = [1; 3];
ap = Q*a;
disp(ap)

  2.366
  2.0981
  • 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 [25]:
# vary Theta from 0 to 2pi radians
theta_a = linspace(0, 2*pi, 100); # Theta in radians

# Elements of vector in initial frame: theta = 0
a = [1, 3]';

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

# Iterate over theta values and assign the values to ap1_a and ap2_a
for i=1:length(theta_a)
    theta = theta_a(i);
    Q = [cos(theta), sin(theta); -sin(theta), cos(theta)];
    ap = Q*a;
    ap1_a(i) = ap(1); # 1st element, index = 1
    ap2_a(i) = ap(2); # 2nd element, index = 2
end
  • We can now plot the elements as function of theta
In [26]:
set(0, "defaultlinelinewidth", 5);
set (0, "defaulttextfontname", "TimesNewRoman")
set (0, "defaulttextfontsize", 20)
set (0, "DefaultAxesFontName", "TimesNewRoman")
set(0, 'DefaultAxesFontSize', 20)

plot(theta_a, ap1_a, "-;a'_1;");
hold on;
plot(theta_a, ap2_a, "-;a'_2;");
grid on;
xlabel("\\theta");
ylabel("Components of a'");
h = legend();
legend (h, "location", "southwest");
legend boxoff;

set (h, "fontsize", 16,"fontname", "TimesNewRoman");

hold off;
  • 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 [27]:
plot(ap1_a, ap2_a);
hold on;
plot(a(1), a(2), 'ok'); # initial vector
daspect([1 1]);
hold off;

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 [28]:
theta = 30*pi/180;
Q =[cos(theta), sin(theta);-sin(theta), cos(theta)];

# Original state of stress
s = [50, 30; 30, -20];
disp(s)

# Transformed matrix
sp =(Q*s)*Q';
disp(sp)

  50  30
  30  -20
  58.481  -15.311
  -15.311  -28.481
  • Let's loop over $ \theta $ and see how state of stress changes
In [29]:
theta_a = linspace(0, pi, 100); # Theta in radians

sigma_n_a = zeros(size(theta_a));
tau_n_a = zeros(size(theta_a));

for i=1:length(theta_a)
    theta = theta_a(i);
    Q =[cos(theta), sin(theta);-sin(theta), cos(theta)];
    s = [50, 30; 30, -20];
    sp =(Q*s)*Q';
    sigma_n_a(i) = sp(1,1);
    tau_n_a(i) = sp(1,2);
end
  • Let's plot the Mohr's circle
In [30]:
plot(sigma_n_a, tau_n_a);
hold on;
plot(trace(s)/2, 0, 'ok'); # center of the circle
xlabel("\\sigma_n'");
ylabel("\\tau_n'");

daspect([1 1]);
hold off;
In [ ]: