Coupled Spring-Mass

Consider \(N\) identical masses attached to identical springs on both sides. The springs at the extreme ends are attached to rigid supports. Suppose one of the mass \(m_i\) is given a small oscillation. This will set all the masses into oscillation. The equation of motion of the \(i\)th mass can be written as \begin{equation} m_i\frac{d^2x_i}{dt^2} = F_{left} + F_{right} + F_i \end{equation} For small oscillation, the force acting on a mass may be considered due to restoring force due to spring. Here, \begin{align} &F_{left} = -k\left(x_i - x_{i-1}\right)\\ &F_{right} = -k\left(x_i - x_{i+1}\right) \end{align} where \(x_i\), \(x_{i-1}\), and \(x_{i+1}\) are the displaments of \(i\)th, \((i-1)\)th and \((i+1)\)th masses from the equilibrium postion. So, the equation of motion can be written as \begin{equation} m_i\frac{d^2x_i}{dt^2} = -k\left(x_i - x_{i-1}\right) + -k\left(x_i - x_{i+1}\right) + F_i \end{equation} There are exceptions for the first and last masses, since, on the left of the first mass and to the right of the last mass, there are rigid supports. So, the equation of motion for the first and last masses are \begin{align} &m_1\frac{d^2x_1}{dt^2} = -2kx_1 - kx_2 + F_1\\ &m_N\frac{d^2x_N}{dt^2} = -kx_{n-1} - 2kx_N + F_N \end{align}

Coupled Spring-Mass System (3 masses)

Consider three identical masses are attached to identical springs on either side and the extreme ends of the springs are attached to rigid walls. Suppose the first mass is given a small displacement by an external force which is harmonic in nature. Let this force be \(F=Ce^{i\omega t}\). This will set all the masses into oscillation. Let \(x_1\), \(x_2\), and \(x_3\) are respectively displacements of the springs from equlibrium position. The equation of motion of the masses can be written as follows -

\[m\frac{d^2x}{dt^2} = F_{left} + F_{rest} + F_{ext}\]

where \(F_{left}\) and \(F_{rest}\) stands for restoring forces acting on a mass due to spring on its left and right. Since the oscillation is small, the restoring force can be given by Hook’s law. \(F_{ext}\) any external driving force which sets the masses into oscillation.

1st mass

\[m\frac{d^2x_1}{dt^2} = -kx_1 - k(x_1 - x_2) + F_1\]

after rearrangement

\[m\frac{d^2x_1}{dt^2} = -2kx_1 + kx_2\]

---

2nd mass

\[m\frac{d^2x_2}{dt^2} = -k(x_2 - x_1) - k(x_2 - x_3)\]

after rearrangement

\[m\frac{d^2x_1}{dt^2} = -kx_1 - 2kx_2 + kx_3\]

---

3rd mass

\[m\frac{d^2x}{dt^2} = -k(x_3 - x_2) - kx_3\]

after rearrangement

\[m\frac{d^2x}{dt^2} = kx_2 - 2kx_3\]

---

Since the masses are set into oscillation by an external force which is harmonic in nature, the displacement of each mass can be written as \(x(t)=xe^{i\omega t}\). Substituting this in the above equations, we get :nbsphinx-math:`begin{align}

&- momega^2 x_1 = -2k x_1 + kx_2 + C \ &- m omega^2 x_2 = k x_1 -2 k x_2 + k x_3\ &- m omega^2 x_3 = k x_2 -2k x_3

end{align}` This after simplification becomes \begin{align} &(2k - m \omega^2) x_1 - k x_2 = C \\ &-k x_1 + (2k - m \omega^2)x_2 - k x_3 = 0 \\ &- k x_2 + (2k - m \omega^2) x_3 = 0 \end{align} writing \(2k - m\omega^2 = \alpha\), we can write the system of equations in the matrix form as \begin{equation} \begin{pmatrix} \alpha & -k & 0 \\ -k & \alpha & -k \\ 0 & -k & \alpha \end{pmatrix} \begin{pmatrix} x_1\\ x_2\\ x_3 \end{pmatrix} = \begin{pmatrix} C \\ 0 \\ 0 \end{pmatrix} \end{equation} In the matrix form, this can be written as \begin{equation} AX=b \end{equation}

Gauss Elimination Method

def GaussElimination(A, b):
    x = np.copy(b)
    n = len(b)
    #Elimination phase
    for k in range(n-1):             # Fixed row/ Column no.
        for i in range(k+1, n):      # Row no. to be transformed
            if A[i, k] != 0:
                factor = A[i, k]/A[k, k]
                A[i, : ] = A[i, : ] - factor*A[k, : ]
                b[i] = b[i] - factor*b[k]

    # Back substitution
    for k in range(n-1, -1, -1):
        x[k] = (b[k] - np.dot(A[k, k+1:n], x[k+1:n]))/A[k,k]

    return x

import numpy as np

# Constants
N = 3
k = 6
C = 1
omega = 2
m = 1
alpha = 2*k - m*omega**2

A = np.eye(N, k=-1)*(-k) + np.eye(N)*alpha + np.eye(N, k=1)*(-k)
b = np.array([C, 0, 0], dtype=float)

x = GaussElimination(A, b)
print(x)

[-0.4375 -0.75   -0.5625]

def GaussSeidel(A, b, maxiter=10, tol=1e-6):
    row, col = A.shape
    x = np.zeros_like(b)
    for iter in range(maxiter):
        for i in range(row):
            x_old = x.copy()
            j = i
            x[i] = (b[i] - np.sum(A[i, :j]*x[:j]) - np.sum(A[i, j+1:]*x[j+1:]))/A[i,i]
            err = np.abs(x_old -x)
            if all(errval < tol for errval in err) == True:
                break
    return x

x = GaussSeidel(A, b, maxiter=10, tol=1e-6)
print(x)

[-0.4375 -0.75   -0.5625]