Linear Algebra 2: Enter The Matrix

Linear Algebra 2: Enter The Matrix#

Functions in QM seen as vectors#

Take any function and evaluate at $N$ points e.g $x=(1,2,3, ...N)$.
We just turned out function into a vector with N components or N-dimensional vector
Example: $f(x)=x^2$ can be represented as $(1, 4, 8)$ where we evaluated function at three points.

Operators in QM are seen as Matrices#

In quantum mechanics, functions that describe quantum states can be treated as vectors in a vector space, and operators that act on these functions are represented as matrices.
For instance, a differentiation operator can be represented by a matrix when acting on discretized versions of the wavefunctions. Let’s illustrate this with the first and second derivative operators.

Example: First Derivative as a Matrix#

Suppose we discretize a function $\psi(x)$ over a set of points $\{x_1, x_2, \dots, x_n\}$ with spacing $h$. We can approximate the first derivative using the finite difference:

\[ \frac{d\psi}{dx} \approx \frac{\psi(x_{i+1}) - \psi(x_i)}{h} \]

This operation can be represented by a matrix. For $n = 5$ grid points, the first derivative matrix $D_1$ in finite difference apporxmation looks like:

\[\begin{split} D_1 = \frac{1}{h} \begin{pmatrix} -1 & 1 & 0 & 0 & 0 \\ 0 & -1 & 1 & 0 & 0 \\ 0 & 0 & -1 & 1 & 0 \\ 0 & 0 & 0 & -1 & 1 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \end{split}\]

For a function represented as a vector $\psi = (\psi(x_1), \psi(x_2), \dots, \psi(x_5))$, applying this matrix gives the derivative at each point (except the boundary).

Numerical example of Applying the Derivative Matrix for $\psi(x) = x^2$#

Let’s consider the function $\psi(x) = x^2$ sampled at 5 points $x = (0, 1, 2, 3, 4)$, giving us:

\[\begin{split} \psi = \begin{pmatrix} 0 \\ 1 \\ 4 \\ 9 \\ 16 \end{pmatrix} \end{split}\]

If the spacing between points is $h = 1$, applying the first derivative matrix $D_1$ gives:

\[\begin{split} D_1 \cdot \psi = \frac{1}{1} \begin{pmatrix} -1 & 1 & 0 & 0 & 0 \\ 0 & -1 & 1 & 0 & 0 \\ 0 & 0 & -1 & 1 & 0 \\ 0 & 0 & 0 & -1 & 1 \\ 0 & 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0 \\ 1 \\ 4 \\ 9 \\ 16 \end{pmatrix} = \begin{pmatrix} 1 \\ 3 \\ 5 \\ 7 \\ 0 \end{pmatrix} \end{split}\]

This approximates the first derivative of $\psi(x) = x^2$ at each of the sampled points:

At $x = 1$, the derivative is approximately 1 (exact value: 2).
At $x = 2$, the derivative is approximately 3 (exact value: 4).
At $x = 3$, the derivative is approximately 5 (exact value: 6).
At $x = 4$, the derivative is approximately 7 (exact value: 8).
The last value is zero due to boundary handling.

This illustrates how the finite difference method approximates the derivative for the function $\psi(x) = x^2$ at discrete points.

import numpy as np

# Define the grid points and function
n = 5
x = np.linspace(0, 4, n)  # Grid points from 0 to 4
h = x[1] - x[0]  # Spacing between points
psi = x**2  # Function psi(x) = x^2

# Define the first derivative matrix D_1
D_1 = (1/h) * np.array([
    [-1, 1, 0, 0, 0],
    [0, -1, 1, 0, 0],
    [0, 0, -1, 1, 0],
    [0, 0, 0, -1, 1],
    [0, 0, 0, 0, 0]
])

# Apply the first derivative matrix to the function vector
derivative_psi = D_1 @ psi

print('derivative of x^2 using finite difference matrix', derivative_psi)  

# Exact derivative of psi = x^2 is 2*x
exact_derivative_psi = 2 * x
print('\nExact derivative of x^2:\n', exact_derivative_psi)

derivative of x^2 using finite difference matrix [1. 3. 5. 7. 0.]

Exact derivative of x^2:
 [0. 2. 4. 6. 8.]

Centered Difference approximation#

A much better approximation is obtained by using centered difference (difference between one point to the left and one point to the right from a give point)
And taking care of boundary points via higher order evaluation of derivatives.

\[ \frac{d\psi}{dx} \approx \frac{\psi(x_{i+1}) - \psi(x_{i-1})}{2h} \]

This operation can be represented by a matrix. For $n = 5$ grid points, with improved boundary handling, the first derivative matrix $D_1$ using the centered difference method looks like:

\[\begin{split} D_1 = \frac{1}{2h} \begin{pmatrix} -3 & 4 & -1 & 0 & 0 \\ -1 & 0 & 1 & 0 & 0 \\ 0 & -1 & 0 & 1 & 0 \\ 0 & 0 & -1 & 0 & 1 \\ 0 & 0 & 1 & -4 & 3 \end{pmatrix} \end{split}\]

The coefficients $-3,4,-1$ come from applying a Taylor series expansion to derive a higher-order accurate one-sided difference approximation for the derivative at the boundary point.

import numpy as np

# Define the grid points and function
n = 5
x = np.linspace(0, 4, n)  # Grid points from 0 to 4
h = x[1] - x[0]  # Spacing between points
psi = x**2  # Function psi(x) = x^2

# Define the centered difference approximation for the first derivative matrix
D_1_centered = (1/(2*h)) * np.array([
    [-3, 4, -1, 0, 0],  # One-sided difference for the first point (improved boundary handling)
    [-1, 0, 1, 0, 0],   # Centered difference for interior points
    [0, -1, 0, 1, 0],
    [0, 0, -1, 0, 1],
    [0, 0, 1, -4, 3]    # One-sided difference for the last point (improved boundary handling)
])

# Apply the centered difference matrix to the function vector
derivative_psi_centered = D_1_centered @ psi

# Print the results
print('First derivative using centered difference approximation:\n', derivative_psi_centered)

# Exact derivative of psi = x^2 is 2*x
exact_derivative_psi = 2 * x
print('\nExact derivative of x^2:\n', exact_derivative_psi)

First derivative using centered difference approximation:
 [0. 2. 4. 6. 8.]

Exact derivative of x^2:
 [0. 2. 4. 6. 8.]

Second Derivative as a Matrix (Centered Difference)#

Similarly, the second derivative can be approximated by a centered difference formula for interior points:

\[ \frac{d^2\psi}{dx^2} \approx \frac{\psi(x_{i+1}) - 2\psi(x_i) + \psi(x_{i-1})}{h^2} \]

At the boundaries, we use one-sided difference formulas to maintain accuracy:

At the first point $x_1$, we approximate using:

\[ \frac{d^2\psi}{dx^2} \bigg|_{x_1} \approx \frac{\psi(x_1) - 2\psi(x_2) + \psi(x_3)}{h^2} \]
At the last point $x_n$, we approximate using:

\[ \frac{d^2\psi}{dx^2} \bigg|_{x_n} \approx \frac{\psi(x_{n-2}) - 2\psi(x_{n-1}) + \psi(x_n)}{h^2} \]

The corresponding second derivative matrix $D_2$ for $n = 5$ grid points, with this improved boundary handling, is:

\[\begin{split} D_2 = \frac{1}{h^2} \begin{pmatrix} 1 & -2 & 1 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 \\ 0 & 1 & -2 & 1 & 0 \\ 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 1 & -2 & 1 \end{pmatrix} \end{split}\]

Applying this matrix to the vector $\psi$ gives the second derivative of the function at each point, including the boundary points with improved accuracy.

second derivative of x^2 by using centered difference 
 [  2.   2.   2.   2. -23.]

Matrix Multiplication#

Matrix multiplicaiton as row by column product

\[\bf C = AB\]

\[C_{ij} = \sum_k A_{ik} B_{kj}\]

Powers of matrix Multiplying same matrix we end up with different powers. E.g second derivative matrix can be obtained by multiplying two first order derivative matrices

\[A \cdot A = A^2\]

Identity matrix I is a special matrix with ones on diagonal and zeroes everywhere else. Mutliplying by I leaves matrix intact

\[I\cdot A = A\cdot I = A\]

Inverse of a Matrix A is a matrix $A^{-1} multiplication with which returns identity matrix:

\[A^{-1}\cdot A = I\]

Why matrix mulitplication order matters#

In linear algebra, matrix multiplication is generally not commutative, meaning that for two matrices $A$ and $B$, it’s usually the case that:

\[ A B \neq B A \]

To illustrate why this happens, let’s consider an example with 2D matrices.

Example: Non-commuting Matrices#

Let’s define two 2D matrices $A$ and $B$:

\[\begin{split} A = \begin{pmatrix} 1 & 2 \\ 0 & 1 \end{pmatrix}, \quad B = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \end{split}\]

We will compute $AB$ and $BA$ to show that they are not equal.

Compute $AB$#

\[\begin{split} AB = \begin{pmatrix} 1 & 2 \\ 0 & 1 \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} = \begin{pmatrix} (1 \cdot 0 + 2 \cdot 1) & (1 \cdot 1 + 2 \cdot 0) \\ (0 \cdot 0 + 1 \cdot 1) & (0 \cdot 1 + 1 \cdot 0) \end{pmatrix} = \begin{pmatrix} 2 & 1 \\ 1 & 0 \end{pmatrix} \end{split}\]

Compute $BA$#

\[\begin{split} BA = \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix} \begin{pmatrix} 1 & 2 \\ 0 & 1 \end{pmatrix} = \begin{pmatrix} (0 \cdot 1 + 1 \cdot 0) & (0 \cdot 2 + 1 \cdot 1) \\ (1 \cdot 1 + 0 \cdot 0) & (1 \cdot 2 + 0 \cdot 1) \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & 2 \end{pmatrix} \end{split}\]

Compare $AB$ and $BA$#

We now have:

\[\begin{split} AB = \begin{pmatrix} 2 & 1 \\ 1 & 0 \end{pmatrix}, \quad BA = \begin{pmatrix} 0 & 1 \\ 1 & 2 \end{pmatrix} \end{split}\]

Clearly, $AB \neq BA$.

Why Don’t Matrices Commute?#

The reason matrix multiplication is not commutative generally lies in how matrix multiplication is performed. Matrix multiplication involves combining rows of the first matrix with columns of the second matrix.
These operations usually yield different results because the effect of applying transformations in different orders changes the outcome.
For example: One matrix might represent a shear transformation and another a rotation. Applying the shear first and then the rotation will generally give a different result than applying the rotation first and then the shear.

Special Cases#

There are some special cases where matrices do commute:

Diagonal matrices: If $A$ and $B$ are both diagonal matrices, then $AB = BA$ because each entry in the diagonal simply scales the corresponding entry in the other matrix.
Identity matrix: The identity matrix $I$ commutes with any matrix $A$, i.e., $IA = AI = A$, because multiplying by the identity matrix leaves the original matrix unchanged.

Show code cell source Hide code cell source

import numpy as np
import matplotlib.pyplot as plt

vector = np.array([1,1])

# Define two commuting matrices (scaling matrices)
A = np.array([[2, 0], [0, 0.5]])
B = np.array([[3, 0], [0, 0.75]])

# Compute the products A*B and B*A
AB = A @ B
BA = B @ A

# Apply the transformations
AB_vector = AB @ vector
BA_vector = BA @ vector

# Redefine non-commuting matrices for more distinct transformations
C = np.array([[0, -1], [1, 0]])  # 90-degree rotation matrix
D = np.array([[1, 2], [0, 1]])  # Horizontal shear

# Compute the products C*D and D*C (non-commuting case)
CD = C @ D
DC = D @ C

# Apply the transformations
CD_vector = CD @ vector
DC_vector = DC @ vector

# Create subplots for commuting and non-commuting cases
fig, axs = plt.subplots(2, 2, figsize=(12, 12))

# Plot 1: A * B (commuting matrices)
axs[0, 0].axhline(0, color='black',linewidth=0.5)
axs[0, 0].axvline(0, color='black',linewidth=0.5)
axs[0, 0].quiver(0, 0, vector[0], vector[1], angles='xy', scale_units='xy', scale=1, color='blue', label='Original Vector')
axs[0, 0].quiver(0, 0, AB_vector[0], AB_vector[1], angles='xy', scale_units='xy', scale=1, color='red', label='A * B Vector')
axs[0, 0].set_xlim(-10, 10)
axs[0, 0].set_ylim(-10, 10)
axs[0, 0].set_title('Commuting: A * B')
axs[0, 0].legend()
axs[0, 0].grid(True)
axs[0, 0].set_aspect('equal', adjustable='box')

# Plot 2: B * A (commuting matrices)
axs[0, 1].axhline(0, color='black',linewidth=0.5)
axs[0, 1].axvline(0, color='black',linewidth=0.5)
axs[0, 1].quiver(0, 0, vector[0], vector[1], angles='xy', scale_units='xy', scale=1, color='blue', label='Original Vector')
axs[0, 1].quiver(0, 0, BA_vector[0], BA_vector[1], angles='xy', scale_units='xy', scale=1, color='green', label='B * A Vector')
axs[0, 1].set_xlim(-10, 10)
axs[0, 1].set_ylim(-10, 10)
axs[0, 1].set_title('Commuting: B * A')
axs[0, 1].legend()
axs[0, 1].grid(True)
axs[0, 1].set_aspect('equal', adjustable='box')

# Plot 3: C * D (non-commuting matrices)
axs[1, 0].axhline(0, color='black',linewidth=0.5)
axs[1, 0].axvline(0, color='black',linewidth=0.5)
axs[1, 0].quiver(0, 0, vector[0], vector[1], angles='xy', scale_units='xy', scale=1, color='blue', label='Original Vector')
axs[1, 0].quiver(0, 0, CD_vector[0], CD_vector[1], angles='xy', scale_units='xy', scale=1, color='purple', label='C * D Vector')
axs[1, 0].set_xlim(-10, 10)
axs[1, 0].set_ylim(-10, 10)
axs[1, 0].set_title('Non-Commuting: C * D')
axs[1, 0].legend()
axs[1, 0].grid(True)
axs[1, 0].set_aspect('equal', adjustable='box')

# Plot 4: D * C (non-commuting matrices)
axs[1, 1].axhline(0, color='black',linewidth=0.5)
axs[1, 1].axvline(0, color='black',linewidth=0.5)
axs[1, 1].quiver(0, 0, vector[0], vector[1], angles='xy', scale_units='xy', scale=1, color='blue', label='Original Vector')
axs[1, 1].quiver(0, 0, DC_vector[0], DC_vector[1], angles='xy', scale_units='xy', scale=1, color='orange', label='D * C Vector')
axs[1, 1].set_xlim(-10, 10)
axs[1, 1].set_ylim(-10, 10)
axs[1, 1].set_title('Non-Commuting: D * C')
axs[1, 1].legend()
axs[1, 1].grid(True)
axs[1, 1].set_aspect('equal', adjustable='box')

# Display the plots
plt.show()

../_images/3c7fa93574c11262b4cfab341c9511b3c00b9cf98e5a8478da9ad4ad39b41c2a.png

Conclusion#

By treating functions as vectors and operators as matrices, quantum mechanics becomes a powerful framework for solving complex problems. Operations like differentiation, integration, and multiplication by physical observables can be represented as matrix operations, which leads to efficient numerical methods in computational quantum mechanics.