4.2 Linear Algebra

NumPy provides robust support for linear algebra operations through its numpy.linalg module. This is essential for mechanics, where vectors and matrices are the primary language.

Marimo notebook

Check out the Marimo notebook for interactive examples of linear algebra operations in mechanics: Marimo Linear Algebra Notebook.

Vector Operations

Dot Product

The dot product is a fundamental operation for projections and calculating work.

import numpy as np
from mechanicskit import la, ltx, labeled

a = np.array([1, 2, 3])
b = np.array([4, 5, 6])

# Dot product
dot_prod = np.dot(a, b)
# OR using the @ operator
dot_prod_op = a @ b

print(f"Dot product: {dot_prod}")

Dot product: 32

Cross Product

The cross product is used for calculating moments and normal vectors.

# Cross product
cross_prod = np.cross(a, b)
print(f"Cross product: {cross_prod}")

Cross product: [-3  6 -3]

Outer Product

The outer product is useful for constructing matrices from vectors.

# Outer product
outer_prod = np.outer(a, b)
ltx(r"\mathbf{a} \otimes \mathbf{b}  = ", a, " \\otimes ", b, " = ", outer_prod)

\[ \mathbf{a} \otimes \mathbf{b} = \begin{bmatrix}1 \\ 2 \\ 3\end{bmatrix} \otimes \begin{bmatrix}4 \\ 5 \\ 6\end{bmatrix} = \begin{bmatrix}4 & 5 & 6 \\ 8 & 10 & 12 \\ 12 & 15 & 18\end{bmatrix} \]

Norms

Calculating the magnitude (length) of a vector.

# L2 Norm (Euclidean length)
magnitude = np.linalg.norm(a)
print(f"Magnitude of a: {magnitude}")

Magnitude of a: 3.7416573867739413

Matrix Operations

Transpose

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

ltx(r"\mathbf{A} = ", A)

\[ \mathbf{A} = \begin{bmatrix}1 & 2 \\ 3 & 4\end{bmatrix} \]

AT = A.T

ltx(r"\mathbf{A}^\mathsf T = ", AT)

\[ \mathbf{A}^\mathsf T = \begin{bmatrix}1 & 3 \\ 2 & 4\end{bmatrix} \]

Matrix Multiplication

It is crucial to distinguish between element-wise multiplication (*) and matrix multiplication (@ or np.matmul).

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

# Element-wise multiplication
print(f"Element-wise (A * B):\n{A * B}")

# Matrix multiplication
print(f"Matrix multiplication (A @ B):\n{A @ B}")

Element-wise (A * B):
[[ 5 12]
 [21 32]]
Matrix multiplication (A @ B):
[[19 22]
 [43 50]]

Determinants

The determinant is a scalar value that can be computed from the elements of a square matrix. It provides important information about the matrix, such as whether it is invertible.

# Calculate determinant
det_A = np.linalg.det(A)
print(f"Determinant of A: {det_A}")

Determinant of A: -2.0000000000000004

Solving Linear Systems

Solving systems of equations \(Ax = b\).

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

# Solve for x
x = np.linalg.solve(A, b)
print(f"Solution x: {x}")

# Verify solution
print(f"Check A@x: {A @ x}")

Solution x: [2. 3.]
Check A@x: [9. 8.]

Eigenvalues and Eigenvectors

Crucial for principal stress analysis, stability analysis, and vibrations.

# Define a symmetric matrix (e.g., stress tensor)
S = np.array([[10.0, 2.0], [2.0, 5.0]])

# Calculate eigenvalues and eigenvectors
eigenvalues, eigenvectors = np.linalg.eig(S)

print(f"Eigenvalues: {eigenvalues}")
print(f"Eigenvectors:\n{eigenvectors}")

Eigenvalues: [10.70156212  4.29843788]
Eigenvectors:
[[ 0.94362832 -0.33100694]
 [ 0.33100694  0.94362832]]

Rotations

Rotations in 3D space can be represented using rotation matrices.

theta = np.radians(45) # Convert degrees to radians
c, s = np.cos(theta), np.sin(theta)

# Rotation matrix around Z-axis
R_z = np.array([
    [c, -s, 0],
    [s,  c, 0],
    [0,  0, 1]
])

v = np.array([1, 0, 0])
v_rotated = R_z @ v
print(f"Rotated vector: {v_rotated}")

Rotated vector: [0.70710678 0.70710678 0.        ]

Rotations in 2D - a special case

Rotating in the 2D plane is just a special case of the 3D rotation around the z-axis. We can see from the matrix \(\mathbf R_z\) that we have (in 2D)

\[ \mathbf{R}=\begin{bmatrix}c & -s\\ s & c \end{bmatrix} \]

a counter-clockwise rotation. We verify this by rotating the \([1,0]\) vector 90 degrees…

import numpy as np
theta = np.radians(90) # Convert degrees to radians
c, s = np.cos(theta), np.sin(theta)
R = np.array([
    [c, -s],
    [s,  c],
])

v = np.array([1, 0,])
v_rotated = R @ v
print(f"Rotated vector: {v_rotated.round(2)}")

Rotated vector: [0. 1.]