Numerical Linear Algebra Lab#
Iterative Solvers: Update Jacobi to Gauss-Seidel#
We learned in the lecture about the Jacobi interation method for solving \(\mathbf{x}\) in the linear equation \(A \mathbf{x} = \mathbf{b}\). The algorithm states that by applying the update rule:
\(x_i^{(k+1)}\) will converge to the solution \(\mathbf{x}\).
We implemeneted and test it in python:
def jacobi_iteration(A, b, max_iter=1000, tol=1e-8):
"""
Solve A x = b using Jacobi iteration.
A is assumed to be square with non-zero diagonal.
"""
n = A.shape[0]
x = np.zeros(n)
for k in range(max_iter):
x_old = np.copy(x)
for i in range(n):
# Sum over off-diagonal terms
s = 0.0
for j in range(n):
if j != i:
s += A[i,j] * x_old[j]
x[i] = (b[i] - s) / A[i, i]
# Check for convergence
if np.linalg.norm(x - x_old, ord=np.inf) < tol:
print(f"Jacobi converged in {k+1} iterations.")
return x
print("Jacobi did not fully converge within max_iter.")
return x
import numpy as np
A = np.array([[4.0, -1.0, 0.0],
[-1.0, 4.0, -1.0],
[0.0, -1.0, 4.0]], dtype=float)
b = np.array([6.0, 6.0, 6.0], dtype=float)
x_jacobi = jacobi_iteration(A, b)
print("Jacobi solution:", x_jacobi)
print("Direct solve comparison:", np.linalg.solve(A, b))
Jacobi converged in 20 iterations.
Jacobi solution: [2.14285714 2.57142857 2.14285714]
Direct solve comparison: [2.14285714 2.57142857 2.14285714]
For the first hands-on, please update the above Jacobi algorithm and turn it into a working Gauss Seidel algorithm.
Recalling from the lecture that the update step in the Gauss Seidel algorithm is:
Note that \(\mathbf{x}^{(k+1)}\) is already partially updated (for \(j < i\)).
# HANDSON: update the above Jacobi algorithm to a working Gauss Seidel algorithm
def gauss_seidel_iteration(A, b, max_iter=1000, tol=1e-8):
"""
Solve A x = b using the Gauss-Seidel iteration.
A is assumed to be square with non-zero diagonal.
"""
pass
x_gs = gauss_seidel_iteration(A, b)
print("Jacobi solution:", x_jacobi)
print("Gauss-Seidel solution:", x_gs)
print("Direct solve comparison:", np.linalg.solve(A, b))
Jacobi solution: [2.14285714 2.57142857 2.14285714]
Gauss-Seidel solution: None
Direct solve comparison: [2.14285714 2.57142857 2.14285714]
Coupled Harmonic Oscillators#
Harmonic oscillators are a classic problem in physics, often forming the basis of more complex vibrational analyses. In this part, you will:
Derive or reference the analytical solution for two coupled oscillators.
Numerically solve the same system (using an eigenvalue approach).
Generalize to \(n\) (and even \(n \times n\)) coupled oscillators, visualizing the mode shapes.
Two Coupled Oscillators–Analytical Solution#
Consider two masses \(m\) connected by three springs (constant \(k\)), arranged in a line and connected to two walls:
|--k--[m]--k--[m]--k--|
If each mass can move only horizontally, the equations of motion form a \(2 \times 2\) eigenvalue problem.
Let:
\(x_1(t)\) be the horizontal displacement of Mass 1 from its equilibrium position.
\(x_2(t)\) be the horizontal displacement of Mass 2.
We assume small oscillations, so Hooke’s law applies linearly.
Mass 1 experiences:
A restoring force \(-k \,x_1\) from the left wall spring.
A coupling force from the middle spring: if \(x_2 > x_1\), that spring pulls Mass 1 to the right; if \(x_1 > x_2\), it pulls Mass 1 to the left. The net contribution is \(-k (x_1 - x_2)\).
Summing forces (Newton’s second law) gives:
(38)#\[\begin{align} m \, \ddot{x}_1 = -k\,x_1 - k\,(x_1 - x_2). \end{align}\]Mass 2 experiences:
A restoring force \(-k\,x_2\) from the right wall spring.
The coupling force from the middle spring: \(-k(x_2 - x_1)\).
Hence,
(39)#\[\begin{align} m \, \ddot{x}_2 = -k\,x_2 - k\,(x_2 - x_1). \end{align}\]
Rewrite each equation:
We can write \(\mathbf{x} = \begin{pmatrix}x_1 \\ x_2\end{pmatrix}\) and express the system as:
where
Equivalently,
This is a second-order linear system describing small oscillations.
We look for solutions of the form
where \(\mathbf{x}(0)\) is the initial condition and \(\omega\) is the (angular) oscillation frequency.
Plugging into \(m\,\ddot{\mathbf{x}} + K\,\mathbf{x} = 0\) gives:
Nontrivial solutions exist only if
which is the eigenvalue problem for \(\omega^2\).
Explicitly, \(K - m\,\omega^2 I\) is:
The determinant is:
Setting this to zero results
Hence, we get two solutions for \(\omega^2\):
\(\omega_+^2\): taking the \(+\) sign:
(50)#\[\begin{align} 2k - m\,\omega_+^2 = k \quad \Longrightarrow \quad \omega_+^2 = \frac{k}{m} \quad \Longrightarrow \quad \omega_1 = \sqrt{\frac{k}{m}}. \end{align}\]\(\omega_-^2\): taking the \(-\) sign:
(51)#\[\begin{align} 2k - m\,\omega_-^2 = -\,k \quad \Longrightarrow \quad \omega_-^2 = \frac{3k}{m} \quad \Longrightarrow \quad \omega_2 = \sqrt{\frac{3k}{m}}. \end{align}\]
For each of the normal modes:
Lower Frequency \(\omega_+ = \sqrt{k/m}\): Plug \(\omega_+^2 = k/m\) back into \((K - m\,\omega_+^2 I)\mathbf{X} = 0\). For instance,
(52)#\[\begin{align} \begin{pmatrix} 2k - k & -k \\ -k & 2k - k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} k & -k \\ -k & k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \mathbf{0}. \end{align}\]This implies \(x_1 = x_2\). Physically, the in-phase mode has both masses moving together.
Higher Frequency \(\omega_- = \sqrt{3k/m}\):
(53)#\[\begin{align} \begin{pmatrix} 2k - 3k & -k \\ -k & 2k - 3k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \begin{pmatrix} -k & -k \\ -k & -k \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \end{pmatrix} = \mathbf{0}. \end{align}\]This yields \(x_1 = -x_2\). Physically, the out-of-phase mode has the two masses moving in opposite directions.
We can compute the position of these coupled oscillators according to the analytical solutions.
# Physical constants
m = 1.0 # mass
k = 1.0 # spring constant
# Frequencies for two normal modes
omegap = np.sqrt(k/m) # in-phase
omegam = np.sqrt(3*k/m) # out-of-phase
# Initial conditions
x1_0 = 0
x2_0 = 0.5
# The analytical solution:
def X_analytic(t):
xp_0 = (x1_0 + x2_0) / 2
xm_0 = (x1_0 - x2_0) / 2
xp = xp_0 * np.cos(omegap * t)
xm = xm_0 * np.cos(omegam * t)
return xp + xm, xp - xm
Plot multiple frames:
from pathlib import Path
import matplotlib.pyplot as plt
import matplotlib.animation as animation
def mkplots(X, t_max=10, n_frames=201):
T = np.linspace(0, t_max, n_frames)
Path("plots").mkdir(parents=True, exist_ok=True)
for i, t in enumerate(T):
x1, x2 = X(t)
plt.plot([-2,-1+x1,1+x2,2], [0,0,0,0], 'o-')
plt.xlim(-2,2)
plt.savefig(f"plots/{i:04}.png")
plt.close()
mkplots(X_analytic)
Can combine them into a movie using ffmpeg:
! rm -rf movie.mpg && ffmpeg -i plots/%04d.png -qmax 2 movie.mpg
/usr/bin/sh: 1: ffmpeg: not found
Two Coupled Oscillators–Semi-Analytical/Numerical Solution#
Instead of solving the coupled oscillator problem analytically, we can at least solve the eigenvalue part of the problem numerically.
# HANDSON: Step 1. rewrite the analytical solution in matrix form
# Physical constants
m = 1.0 # mass
k = 1.0 # spring constant
# Frequencies for two normal modes
Omega = np.array([...]) # this should become a numpy array
# Initial conditions
X0 = np.array([...]) # this should become a numpy array
# The analytical solution in matrix notation:
def X_matrix(t):
M0 = ... # apply an transformation to rewrite the transformation in terms of eigenvectors
M = M0 * np.cos(Omega * t)
return ... # apply an inverse transformation to rewrite the modes in terms of x1 and x2
# HANDSON: Step 2. Replace the manual solutions of eigenvalues Omega and the transform by calling `np.linalg.eig()`
# Coupling matrix
K = np.array([
[...],
[...],
])
# Initial conditions
X0 = np.array([...])
# The analytical solution in matrix notation:
def X_matrix(t):
...
Omega = ...
M0 = ...
M = M0 * np.cos(Omega * t)
return ...
# HANDSON: Step 3. Generalize the solution to work for arbitrary number of coupled oscillators
# Coupling matrix
K = np.array([
[...],
[...],
])
# Initial conditions
X0 = np.array([...])
# The analytical solution in matrix notation:
def X_matrix(t):
...
Omega = ...
M0 = ...
M = M0 * np.cos(Omega * t)
return ...
# Turn the outputs into a movie
# mkplots(X)