A vibrating mechanical system has degrees of freedom (DOF) equal to the minimum number of independent coordinates required to describe its configuration. For a lumped translational model, we collect displacements into the coordinate vector \( \mathbf{q}(t) \in \mathbb{R}^n \), where \( n \) is the DOF count.

In this lesson we restrict to the classical linear vibration regime:

• Small motions so geometry does not change appreciably (linear kinematics).
• Linear springs and viscous dampers (force proportional to relative displacement/velocity).
• Lumped parameters so the model is a finite-dimensional ODE (not PDE).

Under these assumptions, the governing equations admit the standard matrix form:

\[ \mathbf{M}\ddot{\mathbf{q} }(t) + \mathbf{C}\dot{\mathbf{q} }(t) + \mathbf{K}\mathbf{q}(t) = \mathbf{f}(t), \quad \mathbf{M},\mathbf{C},\mathbf{K} \in \mathbb{R}^{n\times n},\; \mathbf{f}(t)\in\mathbb{R}^{n}. \]

Here \( \mathbf{M} \) is the mass matrix, \( \mathbf{C} \) the damping matrix, and \( \mathbf{K} \) the stiffness matrix. Later lessons will analyze this structure via eigenvalue problems and modal coordinates; in this lesson our focus is how to derive \( \mathbf{M},\mathbf{C},\mathbf{K} \) correctly from physics.

Consider a linear spring–damper element connecting DOF \( i \) and \( j \) with stiffness \( k_e \) and damping \( c_e \). Let the scalar displacement coordinates be \( q_i(t) \) and \( q_j(t) \). Define the relative displacement and velocity:

\[ \delta_e(t) = q_i(t) - q_j(t), \qquad \dot{\delta}_e(t) = \dot{q}_i(t) - \dot{q}_j(t). \]

For a passive linear spring and viscous damper, the internal element force opposing relative motion is:

\[ s_e(t) = k_e\,\delta_e(t) + c_e\,\dot{\delta}_e(t). \]

The force contributed by this element on DOF \( i \) is \( -s_e(t) \), and on DOF \( j \) is \( +s_e(t) \) (action–reaction). This “equal and opposite” structure is exactly what produces symmetric coupling matrices in assembled MDOF models.

For each mass coordinate \( q_i \), Newton’s second law yields: \( m_i \ddot{q}_i = \sum \text{(external forces)} + \sum \text{(element forces)} \). By writing this balance for all DOFs and collecting terms in \( \mathbf{q},\dot{\mathbf{q} },\ddot{\mathbf{q} } \), we obtain the matrix equation in Section 1.

Consider the standard 2-DOF chain with displacements \( q_1(t) \) and \( q_2(t) \). Mass 1 is connected to ground by \( (k_1,c_1) \), the masses are connected by \( (k_2,c_2) \), and mass 2 is connected to ground by \( (k_3,c_3) \). External forces are \( f_1(t) \) and \( f_2(t) \).

flowchart LR
  G1["Ground"] --- E1["k1,c1"] --- M1["m1 (q1)"] --- E2["k2,c2"] --- M2["m2 (q2)"] --- E3["k3,c3"] --- G2["Ground"]
  F1["f1(t)"] --> M1
  F2["f2(t)"] --> M2
        

Using the element law from Section 2, the internal force in the coupling element is \( s_2(t)=k_2(q_1-q_2)+c_2(\dot{q}_1-\dot{q}_2) \). The ground elements generate \( s_1(t)=k_1 q_1 + c_1 \dot{q}_1 \) and \( s_3(t)=k_3 q_2 + c_3 \dot{q}_2 \).

Equation for mass 1. Forces on mass 1 oppose its displacement/velocity relative to ground and relative to mass 2:

\[ m_1\ddot{q}_1 = f_1(t) - s_1(t) - s_2(t) = f_1(t) - (k_1 q_1 + c_1 \dot{q}_1) - \bigl(k_2(q_1-q_2)+c_2(\dot{q}_1-\dot{q}_2)\bigr). \]

Equation for mass 2. The coupling element acts in the opposite direction on mass 2:

\[ m_2\ddot{q}_2 = f_2(t) - s_3(t) + s_2(t) = f_2(t) - (k_3 q_2 + c_3 \dot{q}_2) + \bigl(k_2(q_1-q_2)+c_2(\dot{q}_1-\dot{q}_2)\bigr). \]

Rearranging each equation into “all system terms on the left” gives:

\[ \begin{aligned} m_1\ddot{q}_1 + (c_1+c_2)\dot{q}_1 - c_2\dot{q}_2 + (k_1+k_2)q_1 - k_2 q_2 &= f_1(t), \\ m_2\ddot{q}_2 - c_2\dot{q}_1 + (c_2+c_3)\dot{q}_2 - k_2 q_1 + (k_2+k_3)q_2 &= f_2(t). \end{aligned} \]

In matrix form \( \mathbf{M}\ddot{\mathbf{q} }+\mathbf{C}\dot{\mathbf{q} }+\mathbf{K}\mathbf{q}=\mathbf{f}(t) \) with \( \mathbf{q}=[q_1\;q_2]^\top \):

\[ \mathbf{M}=\begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix},\quad \mathbf{C}=\begin{bmatrix} c_1+c_2 & -c_2 \\ -c_2 & c_2+c_3 \end{bmatrix},\quad \mathbf{K}=\begin{bmatrix} k_1+k_2 & -k_2 \\ -k_2 & k_2+k_3 \end{bmatrix},\quad \mathbf{f}(t)=\begin{bmatrix} f_1(t) \\ f_2(t) \end{bmatrix}. \]

Note the physically meaningful pattern: diagonal terms collect all elements attached to a DOF; off-diagonal terms are negative couplings of elements connecting two DOFs.

For an MDOF network of springs and dampers, it is useful to represent each element as connecting two nodes (DOFs). For an element \( e \) connecting DOFs \( i \) and \( j \), define a selector vector \( \mathbf{b}_e \in \mathbb{R}^n \) with \( +1 \) at index \( i \), \( -1 \) at index \( j \), and zeros elsewhere. Then:

\[ \delta_e = q_i - q_j = \mathbf{b}_e^\top \mathbf{q}, \qquad \dot{\delta}_e = \dot{q}_i - \dot{q}_j = \mathbf{b}_e^\top \dot{\mathbf{q} }. \]

The element internal force magnitude is \( s_e = k_e \mathbf{b}_e^\top \mathbf{q} + c_e \mathbf{b}_e^\top \dot{\mathbf{q} } \). The corresponding nodal force contribution is \( \mathbf{f}^{(e)} = -\mathbf{b}_e s_e \). Summing over all elements yields:

\[ \sum_e \mathbf{f}^{(e)} = -\sum_e \mathbf{b}_e\Bigl(k_e \mathbf{b}_e^\top \mathbf{q} + c_e \mathbf{b}_e^\top \dot{\mathbf{q} }\Bigr) = -\Bigl(\sum_e k_e\,\mathbf{b}_e\mathbf{b}_e^\top\Bigr)\mathbf{q} -\Bigl(\sum_e c_e\,\mathbf{b}_e\mathbf{b}_e^\top\Bigr)\dot{\mathbf{q} }. \]

Therefore, for a passive spring–damper network: \( \mathbf{K} = \sum_e k_e\,\mathbf{b}_e\mathbf{b}_e^\top \) and \( \mathbf{C} = \sum_e c_e\,\mathbf{b}_e\mathbf{b}_e^\top \).

Compact matrix factorization. Stack the element vectors into an incidence matrix \( \mathbf{B}\in\mathbb{R}^{m\times n} \) (one row per element, equal to \( \mathbf{b}_e^\top \)). Let \( \mathbf{D}_k=\operatorname{diag}(k_1,\dots,k_m) \) and \( \mathbf{D}_c=\operatorname{diag}(c_1,\dots,c_m) \). Then:

\[ \mathbf{K} = \mathbf{B}^\top \mathbf{D}_k \mathbf{B}, \qquad \mathbf{C} = \mathbf{B}^\top \mathbf{D}_c \mathbf{B}. \]

Proof 1 (Symmetry). Since \( \mathbf{D}_k \) and \( \mathbf{D}_c \) are diagonal and symmetric, we have:

\[ \mathbf{K}^\top = (\mathbf{B}^\top \mathbf{D}_k \mathbf{B})^\top = \mathbf{B}^\top \mathbf{D}_k^\top \mathbf{B} = \mathbf{B}^\top \mathbf{D}_k \mathbf{B} = \mathbf{K}, \]

and similarly \( \mathbf{C}^\top=\mathbf{C} \). Thus passive interconnections produce symmetric stiffness and damping matrices.

Proof 2 (Positive semidefiniteness). Assume all spring constants satisfy \( k_e \ge 0 \). Then for any vector \( \mathbf{v}\in\mathbb{R}^n \):

\[ \mathbf{v}^\top \mathbf{K}\mathbf{v} = \mathbf{v}^\top \mathbf{B}^\top \mathbf{D}_k \mathbf{B}\mathbf{v} = (\mathbf{B}\mathbf{v})^\top \mathbf{D}_k (\mathbf{B}\mathbf{v}) = \sum_{e=1}^{m} k_e \, (\mathbf{b}_e^\top \mathbf{v})^2 \ge 0. \]

Hence \( \mathbf{K}\succeq 0 \). The same argument gives \( \mathbf{C}\succeq 0 \) when \( c_e \ge 0 \). Intuitively, \( \mathbf{v}^\top \mathbf{K}\mathbf{v} \) is twice the strain energy stored by the spring network when the system is displaced by \( \mathbf{v} \).

Mass matrix. For lumped translational DOFs, the mass matrix is typically diagonal: \( \mathbf{M}=\operatorname{diag}(m_1,\dots,m_n) \) with \( m_i > 0 \), implying \( \mathbf{M}\succ 0 \). More general models (e.g., with kinematic couplings or coordinate transformations) can yield symmetric positive definite non-diagonal \( \mathbf{M} \).

In many physical assemblies, viscous damping is naturally elementwise (dashpots between nodes), which directly produces \( \mathbf{C}=\sum_e c_e \mathbf{b}_e\mathbf{b}_e^\top \). Another common engineering approximation is proportional (Rayleigh) damping:

\[ \mathbf{C} = \alpha\,\mathbf{M} + \beta\,\mathbf{K}, \qquad \alpha \ge 0,\; \beta \ge 0. \]

This form is convenient because it preserves structural symmetry and (as shown in later lessons) often enables modal decoupling. In this lesson, treat it as a modeling assumption that must be justified by system identification or experimental evidence.

Numerical simulation is simplest in first-order form. Define the state vector \( \mathbf{x}(t) = \begin{bmatrix}\mathbf{q}(t) \\ \dot{\mathbf{q} }(t)\end{bmatrix} \in \mathbb{R}^{2n} \). Assuming \( \mathbf{M} \) is invertible, solve for acceleration:

\[ \ddot{\mathbf{q} }(t) = \mathbf{M}^{-1}\Bigl(\mathbf{f}(t) - \mathbf{C}\dot{\mathbf{q} }(t) - \mathbf{K}\mathbf{q}(t)\Bigr). \]

Therefore the first-order system \( \dot{\mathbf{x} } = \mathbf{A}\mathbf{x} + \mathbf{B}\mathbf{f}(t) \) becomes:

\[ \dot{\mathbf{x} }(t) = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix}\mathbf{x}(t) + \begin{bmatrix} \mathbf{0} \\ \mathbf{M}^{-1} \end{bmatrix}\mathbf{f}(t). \]

This conversion uses only linear algebra and the state-space concepts covered previously, and it is the bridge between modeling and simulation.

flowchart TD
  A["Choose DOF coordinates q1..qn and sign conventions"] --> B["Draw free-body diagram for each DOF"]
  B --> C["Write Newton balance: mi*qdd_i = sum forces"]
  C --> D["Collect terms into M qdd + C qd + K q = f(t)"]
  D --> E["Sanity checks: symmetry, units, passivity (K,C PSD)"]
  E --> F["Convert to state space: x = [q; qd], xdot = A x + B f(t)"]
  F --> G["Simulate and validate against expected limiting cases"]
        

Recommended scientific stack for System Dynamics in Python: \( \)

• NumPy (matrix/vector computation)
• SciPy (ODE solvers, linear algebra)
• Matplotlib (visualization)
• python-control (optional; state-space utilities consistent with control coursework)
• SymPy (optional; symbolic checks for small models)

File: Chapter13_Lesson1.py

# Chapter13_Lesson1.py
"""System Dynamics — Chapter 13, Lesson 1
Modeling Two- and Multi-Degree-of-Freedom Mechanical Systems.

This script:
1) Builds M, C, K for a 2-DOF mass–spring–damper system,
2) Converts to first-order state-space form,
3) Simulates free and forced responses using SciPy's solve_ivp.

Dependencies:
  pip install numpy scipy matplotlib
"""

from __future__ import annotations
import numpy as np
from numpy.typing import ArrayLike
from dataclasses import dataclass
from typing import Callable, Tuple
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt


@dataclass(frozen=True)
class TwoDOFParams:
    m1: float = 1.0
    m2: float = 0.8
    k1: float = 200.0  # spring to ground on mass 1
    k2: float = 150.0  # spring between masses
    k3: float = 100.0  # spring to ground on mass 2
    c1: float = 1.5    # damper to ground on mass 1
    c2: float = 1.0    # damper between masses
    c3: float = 0.8    # damper to ground on mass 2


def two_dof_mck(p: TwoDOFParams) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """Return (M, C, K) for the standard 2-DOF chain."""
    M = np.diag([p.m1, p.m2])
    C = np.array([[p.c1 + p.c2, -p.c2],
                  [-p.c2, p.c2 + p.c3]], dtype=float)
    K = np.array([[p.k1 + p.k2, -p.k2],
                  [-p.k2, p.k2 + p.k3]], dtype=float)
    return M, C, K


def state_space_from_mck(M: np.ndarray, C: np.ndarray, K: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
    """Build continuous-time state-space matrices for
        M qdd + C qd + K q = f(t).

    State x = [q; qd] ∈ R^{2n}.
    Then xdot = A x + B f, with f ∈ R^{n}.
    """
    n = M.shape[0]
    Minv = np.linalg.solve(M, np.eye(n))
    Z = np.zeros((n, n))
    I = np.eye(n)

    A = np.block([
        [Z, I],
        [-Minv @ K, -Minv @ C]
    ])
    B = np.vstack([Z, Minv])
    return A, B


def simulate(
    A: np.ndarray,
    B: np.ndarray,
    force: Callable[[float], np.ndarray],
    t_span: Tuple[float, float],
    x0: ArrayLike,
    t_eval: np.ndarray | None = None,
) -> solve_ivp:
    """Simulate xdot = A x + B f(t)."""
    x0 = np.asarray(x0, dtype=float).reshape(-1)

    def rhs(t: float, x: np.ndarray) -> np.ndarray:
        return A @ x + B @ force(t)

    sol = solve_ivp(rhs, t_span, x0, t_eval=t_eval, rtol=1e-8, atol=1e-10)
    if not sol.success:
        raise RuntimeError(sol.message)
    return sol


def main() -> None:
    p = TwoDOFParams()
    M, C, K = two_dof_mck(p)
    A, B = state_space_from_mck(M, C, K)

    # ---- Free response (f(t)=0) ----
    def f_free(t: float) -> np.ndarray:
        return np.zeros(2)

    t = np.linspace(0.0, 8.0, 2001)
    x0 = np.array([0.02, -0.01, 0.0, 0.0])  # [q1,q2,q1dot,q2dot]
    sol_free = simulate(A, B, f_free, (t[0], t[-1]), x0, t_eval=t)

    q1, q2 = sol_free.y[0], sol_free.y[1]

    # ---- Forced response (harmonic force on mass 1) ----
    F0 = 5.0
    w = 12.0  # rad/s
    def f_forced(t: float) -> np.ndarray:
        return np.array([F0 * np.cos(w * t), 0.0])

    sol_forced = simulate(A, B, f_forced, (t[0], t[-1]), x0, t_eval=t)
    q1f, q2f = sol_forced.y[0], sol_forced.y[1]

    # ---- Plot ----
    plt.figure()
    plt.plot(t, q1, label="q1 free")
    plt.plot(t, q2, label="q2 free")
    plt.xlabel("t [s]")
    plt.ylabel("displacement [m]")
    plt.grid(True)
    plt.legend()

    plt.figure()
    plt.plot(t, q1f, label="q1 forced")
    plt.plot(t, q2f, label="q2 forced")
    plt.xlabel("t [s]")
    plt.ylabel("displacement [m]")
    plt.grid(True)
    plt.legend()

    plt.show()


if __name__ == "__main__":
    main()
      

Common C++ libraries used for MDOF dynamics in control/mechanical contexts:

• Eigen (linear algebra; header-only)
• Boost::odeint (optional; ODE integration)

File: Chapter13_Lesson1.cpp

// Chapter13_Lesson1.cpp
// System Dynamics — Chapter 13, Lesson 1
// Modeling Two- and Multi-Degree-of-Freedom Mechanical Systems.
//
// Build M, C, K for a 2-DOF mass–spring–damper system, convert to state space,
// and integrate with a simple fixed-step RK4.
//
// Dependencies:
//   - Eigen (header-only): https://eigen.tuxfamily.org
//
// Compile (example):
//   g++ -O2 -std=c++17 Chapter13_Lesson1.cpp -I /path/to/eigen -o ch13_l1
//
#include <iostream>
#include <iomanip>
#include <Eigen/Dense>

using Eigen::MatrixXd;
using Eigen::VectorXd;

struct TwoDOFParams {
    double m1 = 1.0, m2 = 0.8;
    double k1 = 200.0, k2 = 150.0, k3 = 100.0;
    double c1 = 1.5,   c2 = 1.0,   c3 = 0.8;
};

void two_dof_mck(const TwoDOFParams& p, MatrixXd& M, MatrixXd& C, MatrixXd& K) {
    M = MatrixXd::Zero(2,2);
    M(0,0) = p.m1; M(1,1) = p.m2;

    C = MatrixXd(2,2);
    C << (p.c1 + p.c2), -p.c2,
         -p.c2, (p.c2 + p.c3);

    K = MatrixXd(2,2);
    K << (p.k1 + p.k2), -p.k2,
         -p.k2, (p.k2 + p.k3);
}

void state_space_from_mck(const MatrixXd& M, const MatrixXd& C, const MatrixXd& K,
                          MatrixXd& A, MatrixXd& B) {
    const int n = static_cast<int>(M.rows());
    MatrixXd Minv = M.inverse();

    MatrixXd Z = MatrixXd::Zero(n,n);
    MatrixXd I = MatrixXd::Identity(n,n);

    A = MatrixXd::Zero(2*n, 2*n);
    A.block(0,0,n,n) = Z;
    A.block(0,n,n,n) = I;
    A.block(n,0,n,n) = -Minv * K;
    A.block(n,n,n,n) = -Minv * C;

    B = MatrixXd::Zero(2*n, n);
    B.block(n,0,n,n) = Minv;
}

// Fixed-step RK4 for xdot = A x + B f(t)
VectorXd rk4_step(const MatrixXd& A, const MatrixXd& B,
                  const VectorXd& x, double t, double h,
                  const std::function<VectorXd(double)>& force) {

    auto rhs = [&](double tau, const VectorXd& z) {
        return A * z + B * force(tau);
    };

    VectorXd k1 = rhs(t, x);
    VectorXd k2 = rhs(t + 0.5*h, x + 0.5*h*k1);
    VectorXd k3 = rhs(t + 0.5*h, x + 0.5*h*k2);
    VectorXd k4 = rhs(t + h, x + h*k3);
    return x + (h/6.0) * (k1 + 2.0*k2 + 2.0*k3 + k4);
}

int main() {
    TwoDOFParams p;
    MatrixXd M, C, K;
    two_dof_mck(p, M, C, K);

    MatrixXd A, B;
    state_space_from_mck(M, C, K, A, B);

    // Free response: f(t)=0, initial x=[q1,q2,q1dot,q2dot]
    auto f_free = [](double /*t*/) {
        VectorXd f(2); f << 0.0, 0.0; return f;
    };

    double t0 = 0.0, tf = 8.0, h = 0.001;
    int N = static_cast<int>((tf - t0)/h);

    VectorXd x(4);
    x << 0.02, -0.01, 0.0, 0.0;

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "t,q1,q2\n";
    double t = t0;
    for (int i = 0; i <= N; ++i) {
        if (i % 200 == 0) { // print every 0.2 s
            std::cout << t << "," << x(0) << "," << x(1) << "\n";
        }
        x = rk4_step(A, B, x, t, h, f_free);
        t += h;
    }
    return 0;
}
      

Java ecosystems commonly used for MDOF simulation:

• Apache Commons Math (linear algebra + ODE integrators)
• EJML (efficient matrix computations; optional alternative)

File: Chapter13_Lesson1.java

// Chapter13_Lesson1.java
// System Dynamics — Chapter 13, Lesson 1
// Modeling Two- and Multi-Degree-of-Freedom Mechanical Systems.
//
// This example uses Apache Commons Math for:
//   - Linear algebra (RealMatrix)
//   - ODE integration (DormandPrince853Integrator)
//
// Maven dependency (pom.xml):
//   <dependency>
//     <groupId>org.apache.commons</groupId>
//     <artifactId>commons-math3</artifactId>
//     <version>3.6.1</version>
//   </dependency>
//
// Run (example):
//   javac -cp commons-math3-3.6.1.jar Chapter13_Lesson1.java
//   java  -cp .:commons-math3-3.6.1.jar Chapter13_Lesson1
//
import org.apache.commons.math3.linear.*;
import org.apache.commons.math3.ode.FirstOrderDifferentialEquations;
import org.apache.commons.math3.ode.nonstiff.DormandPrince853Integrator;
import org.apache.commons.math3.ode.sampling.FixedStepHandler;
import org.apache.commons.math3.ode.sampling.StepNormalizer;

public class Chapter13_Lesson1 {

    static class TwoDOFParams {
        double m1 = 1.0, m2 = 0.8;
        double k1 = 200.0, k2 = 150.0, k3 = 100.0;
        double c1 = 1.5,   c2 = 1.0,   c3 = 0.8;
    }

    static RealMatrix diag(double a, double b) {
        double[][] d = { {a, 0.0},{0.0, b} };
        return MatrixUtils.createRealMatrix(d);
    }

    static RealMatrix twoDOF_C(TwoDOFParams p) {
        double[][] c = {
            {p.c1 + p.c2, -p.c2},
            {-p.c2, p.c2 + p.c3}
        };
        return MatrixUtils.createRealMatrix(c);
    }

    static RealMatrix twoDOF_K(TwoDOFParams p) {
        double[][] k = {
            {p.k1 + p.k2, -p.k2},
            {-p.k2, p.k2 + p.k3}
        };
        return MatrixUtils.createRealMatrix(k);
    }

    static class StateSpaceODE implements FirstOrderDifferentialEquations {
        private final RealMatrix A;
        private final RealMatrix B;

        StateSpaceODE(RealMatrix A, RealMatrix B) {
            this.A = A;
            this.B = B;
        }

        // Force input f(t) in R^2
        private RealVector force(double t) {
            // Free response: set to zero. Change here for forcing.
            return new ArrayRealVector(new double[]{0.0, 0.0});
        }

        @Override
        public int getDimension() {
            return A.getRowDimension();
        }

        @Override
        public void computeDerivatives(double t, double[] x, double[] xDot) {
            RealVector xv = new ArrayRealVector(x, false);
            RealVector fv = force(t);
            RealVector rhs = A.operate(xv).add(B.operate(fv));
            for (int i = 0; i < xDot.length; i++) {
                xDot[i] = rhs.getEntry(i);
            }
        }
    }

    public static void main(String[] args) {
        TwoDOFParams p = new TwoDOFParams();
        RealMatrix M = diag(p.m1, p.m2);
        RealMatrix C = twoDOF_C(p);
        RealMatrix K = twoDOF_K(p);

        // Build A, B for x = [q; qdot], xdot = A x + B f
        int n = 2;
        RealMatrix Minv = new LUDecomposition(M).getSolver().getInverse();

        RealMatrix Z = MatrixUtils.createRealMatrix(n, n);
        RealMatrix I = MatrixUtils.createRealIdentityMatrix(n);

        RealMatrix A = MatrixUtils.createRealMatrix(2*n, 2*n);
        A.setSubMatrix(Z.getData(), 0, 0);
        A.setSubMatrix(I.getData(), 0, n);
        A.setSubMatrix(Minv.multiply(K).scalarMultiply(-1.0).getData(), n, 0);
        A.setSubMatrix(Minv.multiply(C).scalarMultiply(-1.0).getData(), n, n);

        RealMatrix B = MatrixUtils.createRealMatrix(2*n, n);
        B.setSubMatrix(Minv.getData(), n, 0);

        FirstOrderDifferentialEquations ode = new StateSpaceODE(A, B);

        double t0 = 0.0, tf = 8.0;
        double[] x0 = new double[]{0.02, -0.01, 0.0, 0.0};
        double[] xEnd = new double[x0.length];

        DormandPrince853Integrator integrator =
            new DormandPrince853Integrator(1e-4, 1e-2, 1e-10, 1e-10);

        // Print q1, q2 every 0.2 s
        integrator.addStepHandler(new StepNormalizer(0.2, new FixedStepHandler() {
            @Override
            public void handleStep(double t, double[] x, double[] xDot, boolean isLast) {
                System.out.printf(java.util.Locale.US, "t=%.2f  q1=% .6f  q2=% .6f%n", t, x[0], x[1]);
            }
        }));

        integrator.integrate(ode, t0, x0, tf, xEnd);
        System.out.printf(java.util.Locale.US, "Done. Final state: q1=% .6f  q2=% .6f%n", xEnd[0], xEnd[1]);
    }
}
      

MATLAB toolchain alignment with this lesson:

• ODE suite (ode45/ode15s) for time simulation
• Simulink (State-Space block for continuous-time linear models)
• Control System Toolbox (optional; state-space utilities)

File: Chapter13_Lesson1.m

% Chapter13_Lesson1.m
% System Dynamics — Chapter 13, Lesson 1
% Modeling Two- and Multi-Degree-of-Freedom Mechanical Systems.
%
% This script:
%   1) Defines M, C, K for a 2-DOF system,
%   2) Converts to first-order state-space,
%   3) Simulates free and forced response using ODE45,
%   4) Shows how to use the State-Space block in Simulink (programmatically).

clear; clc;

%% Parameters
m1 = 1.0; m2 = 0.8;
k1 = 200.0; k2 = 150.0; k3 = 100.0;
c1 = 1.5;   c2 = 1.0;   c3 = 0.8;

%% Matrices (M qdd + C qd + K q = f)
M = diag([m1, m2]);
C = [c1+c2, -c2;
     -c2,   c2+c3];
K = [k1+k2, -k2;
     -k2,   k2+k3];

%% State-space: x = [q; qd], xdot = A x + B f
n = 2;
A = [zeros(n), eye(n);
     -M\K,    -M\C];
B = [zeros(n); M\eye(n)];

%% Free response: f(t)=0
force_free = @(t) [0; 0];

rhs = @(t, x) A*x + B*force_free(t);

tspan = [0 8];
x0 = [0.02; -0.01; 0; 0];

opts = odeset('RelTol',1e-8,'AbsTol',1e-10);
[t, x] = ode45(rhs, tspan, x0, opts);

q1 = x(:,1); q2 = x(:,2);

figure; plot(t, q1, t, q2); grid on;
xlabel('t [s]'); ylabel('displacement [m]');
legend('q1 free', 'q2 free');

%% Forced response: harmonic force on mass 1
F0 = 5.0; w = 12.0;
force_forced = @(t) [F0*cos(w*t); 0];

rhs2 = @(t, x) A*x + B*force_forced(t);
[t2, x2] = ode45(rhs2, tspan, x0, opts);

figure; plot(t2, x2(:,1), t2, x2(:,2)); grid on;
xlabel('t [s]'); ylabel('displacement [m]');
legend('q1 forced', 'q2 forced');

%% Simulink note (programmatic model creation)
% You can simulate the same state-space system with Simulink's State-Space block.
% The block expects (A,B,C,D). For example, to output q1 and q2:
% Csim = [eye(2), zeros(2)];
% Dsim = zeros(2,2);
%
% To create a minimal Simulink model automatically:
try
    mdl = 'Chapter13_Lesson1_Simulink';
    if bdIsLoaded(mdl); close_system(mdl, 0); end
    new_system(mdl); open_system(mdl);

    add_block('simulink/Continuous/State-Space', [mdl '/Plant'], 'Position', [120 90 260 160]);
    set_param([mdl '/Plant'], 'A', 'A', 'B', 'B', 'C', 'Csim', 'D', 'Dsim');

    add_block('simulink/Sources/From Workspace', [mdl '/Force'], 'Position', [20 100 90 130]);
    set_param([mdl '/Force'], 'VariableName', 'force_ts');

    add_block('simulink/Sinks/To Workspace', [mdl '/y'], 'Position', [300 105 360 135]);
    set_param([mdl '/y'], 'VariableName', 'yout', 'SaveFormat', 'Array');

    add_line(mdl, 'Force/1', 'Plant/1');
    add_line(mdl, 'Plant/1', 'y/1');

    % Build a time series force input: [f1(t) f2(t)]
    tgrid = linspace(0, 8, 2001)';
    f1 = F0*cos(w*tgrid);
    f2 = zeros(size(tgrid));
    force_ts = timeseries([f1 f2], tgrid);

    Csim = [eye(2), zeros(2)];
    Dsim = zeros(2,2);

    set_param(mdl, 'StopTime', '8');
    sim(mdl);

    fprintf('Simulink ran. yout contains [q1 q2].\n');
catch ME
    fprintf('Simulink section skipped (reason: %s).\n', ME.message);
end
      

Mathematica is particularly effective for small-to-medium DOF symbolic verification and high-accuracy numerical integration via \( \) built-in linear algebra and NDSolve.

File: Chapter13_Lesson1.nb

(* Chapter13_Lesson1.nb
   System Dynamics — Chapter 13, Lesson 1
   Modeling Two- and Multi-Degree-of-Freedom Mechanical Systems.

   This notebook:
     1) Defines M, C, K for a 2-DOF mass–spring–damper system,
     2) Builds the first-order state-space,
     3) Uses NDSolve to simulate free response.
*)

Notebook[{
  Cell["Chapter 13 — Lesson 1: MDOF Modeling (2-DOF Example)", "Section"],

  Cell["1. Define parameters and build M, C, K", "Subsection"],
  Cell[BoxData @ ToBoxes @
    HoldForm[
      m1 = 1.0; m2 = 0.8;
      k1 = 200.0; k2 = 150.0; k3 = 100.0;
      c1 = 1.5; c2 = 1.0; c3 = 0.8;

      M = DiagonalMatrix[{m1, m2}];
      C = { {c1 + c2, -c2}, {-c2, c2 + c3} };
      K = { {k1 + k2, -k2}, {-k2, k2 + k3} };
    ],
  "Input"],

  Cell["2. Convert to state-space: x = [q; qdot]", "Subsection"],
  Cell[BoxData @ ToBoxes @
    HoldForm[
      n = 2;
      Z = ConstantArray[0, {n, n}];
      I = IdentityMatrix[n];
      Minv = Inverse[M];

      A = ArrayFlatten[{ {Z, I}, {-Minv.K, -Minv.C} }];
      B = ArrayFlatten[{ {Z}, {Minv} }];
      (* Force input f(t) in R^2 *)
    ],
  "Input"],

  Cell["3. Free response simulation with NDSolve", "Subsection"],
  Cell[BoxData @ ToBoxes @
    HoldForm[
      f[t_] := {0, 0};

      (* State x(t) = {q1, q2, q1dot, q2dot} *)
      x0 = {0.02, -0.01, 0.0, 0.0};

      eq = x'[t] == A.x[t] + B.f[t];
      ic = x[0] == x0;

      sol = NDSolve[{eq, ic}, x, {t, 0, 8}][[1]];
      q1[t_] := (x[t] /. sol)[[1]];
      q2[t_] := (x[t] /. sol)[[2]];

      Plot[{q1[t], q2[t]}, {t, 0, 8},
           PlotLegends -> {"q1(t)", "q2(t)"},
           AxesLabel -> {"t [s]", "displacement [m]"}]
    ],
  "Input"]
}]
      

Problem 1 (2-DOF Matrix Assembly): For the 2-DOF system in Section 3, derive \( \mathbf{M},\mathbf{C},\mathbf{K} \) explicitly in terms of \( m_1,m_2,k_1,k_2,k_3,c_1,c_2,c_3 \) and write the coupled ODEs.

Solution: From Section 3, the ODEs are:

\[ \begin{aligned} m_1\ddot{q}_1 + (c_1+c_2)\dot{q}_1 - c_2\dot{q}_2 + (k_1+k_2)q_1 - k_2 q_2 &= f_1(t), \\ m_2\ddot{q}_2 - c_2\dot{q}_1 + (c_2+c_3)\dot{q}_2 - k_2 q_1 + (k_2+k_3)q_2 &= f_2(t). \end{aligned} \]

Thus:

\[ \mathbf{M}=\begin{bmatrix} m_1 & 0 \\ 0 & m_2 \end{bmatrix},\quad \mathbf{C}=\begin{bmatrix} c_1+c_2 & -c_2 \\ -c_2 & c_2+c_3 \end{bmatrix},\quad \mathbf{K}=\begin{bmatrix} k_1+k_2 & -k_2 \\ -k_2 & k_2+k_3 \end{bmatrix}. \]

Problem 2 (Passivity Implies \( \mathbf{K}\succeq 0 \)): Suppose a network contains only linear springs with \( k_e \ge 0 \). Prove that the assembled stiffness matrix \( \mathbf{K} \) is symmetric positive semidefinite.

Solution: Using Section 4, \( \mathbf{K}=\mathbf{B}^\top\mathbf{D}_k\mathbf{B} \) with \( \mathbf{D}_k=\operatorname{diag}(k_e) \). Symmetry follows from transposition: \( \mathbf{K}^\top=\mathbf{B}^\top\mathbf{D}_k^\top\mathbf{B}=\mathbf{K} \). For any \( \mathbf{v} \):

\[ \mathbf{v}^\top\mathbf{K}\mathbf{v} =(\mathbf{B}\mathbf{v})^\top \mathbf{D}_k (\mathbf{B}\mathbf{v}) =\sum_e k_e(\mathbf{b}_e^\top\mathbf{v})^2 \ge 0, \]

hence \( \mathbf{K}\succeq 0 \).

Problem 3 (Second-Order to First-Order): Given \( \mathbf{M}\ddot{\mathbf{q} }+\mathbf{C}\dot{\mathbf{q} }+\mathbf{K}\mathbf{q}=\mathbf{f}(t) \) with invertible \( \mathbf{M} \), derive \( \dot{\mathbf{x} }=\mathbf{A}\mathbf{x}+\mathbf{B}\mathbf{f}(t) \) for \( \mathbf{x}=[\mathbf{q};\dot{\mathbf{q} }] \).

Solution: Let \( \mathbf{x}=\begin{bmatrix}\mathbf{q}\\\dot{\mathbf{q} }\end{bmatrix} \). Then \( \dot{\mathbf{x} }=\begin{bmatrix}\dot{\mathbf{q} }\\\ddot{\mathbf{q} }\end{bmatrix} \) and \( \ddot{\mathbf{q} }=\mathbf{M}^{-1}(\mathbf{f}-\mathbf{C}\dot{\mathbf{q} }-\mathbf{K}\mathbf{q}) \). Therefore:

\[ \dot{\mathbf{x} } = \begin{bmatrix} \mathbf{0} & \mathbf{I} \\ -\mathbf{M}^{-1}\mathbf{K} & -\mathbf{M}^{-1}\mathbf{C} \end{bmatrix}\mathbf{x} + \begin{bmatrix} \mathbf{0}\\ \mathbf{M}^{-1} \end{bmatrix}\mathbf{f}(t). \]

Problem 4 (3-DOF Chain Stiffness Matrix): Three masses have coordinates \( q_1,q_2,q_3 \). Springs connect (i) ground–1 with \( k_1 \), (ii) 1–2 with \( k_2 \), (iii) 2–3 with \( k_3 \), (iv) 3–ground with \( k_4 \). Write the stiffness matrix \( \mathbf{K}\in\mathbb{R}^{3\times 3} \).

Solution: Each spring adds \( +k \) to each attached diagonal and \( -k \) to the coupling off-diagonal. Therefore:

\[ \mathbf{K}= \begin{bmatrix} k_1+k_2 & -k_2 & 0 \\ -k_2 & k_2+k_3 & -k_3 \\ 0 & -k_3 & k_3+k_4 \end{bmatrix}. \]

Problem 5 (Coordinate Constraint Elimination): Suppose an MDOF model has coordinate partition \( \mathbf{q}=\begin{bmatrix}\mathbf{q}_f\\\mathbf{q}_c\end{bmatrix} \), where \( \mathbf{q}_c \) are constrained to zero (fixed). Show how to obtain reduced matrices for the free coordinates \( \mathbf{q}_f \).

Solution: Partition matrices conformably: \( \mathbf{M}=\begin{bmatrix}\mathbf{M}_{ff}&\mathbf{M}_{fc}\\\mathbf{M}_{cf}&\mathbf{M}_{cc}\end{bmatrix} \) and similarly for \( \mathbf{C},\mathbf{K} \). With \( \mathbf{q}_c=\dot{\mathbf{q} }_c=\ddot{\mathbf{q} }_c=\mathbf{0} \), the second-order equation reduces to:

\[ \mathbf{M}_{ff}\ddot{\mathbf{q} }_f + \mathbf{C}_{ff}\dot{\mathbf{q} }_f + \mathbf{K}_{ff}\mathbf{q}_f = \mathbf{f}_f(t), \]

i.e., delete constrained rows/columns and keep the free–free blocks (under direct coordinate fixing).

We built MDOF vibration models directly from Newton’s laws and linear element constitutive relations, producing the standard matrix dynamics \( \mathbf{M}\ddot{\mathbf{q} }+\mathbf{C}\dot{\mathbf{q} }+\mathbf{K}\mathbf{q}=\mathbf{f}(t) \). Using an incidence representation, we proved that passive spring/damper networks yield symmetric positive semidefinite \( \mathbf{K} \) and \( \mathbf{C} \). Finally, we converted the second-order model to first-order state space for simulation and implemented the pipeline in five computational environments. The next lesson uses these same matrices to define natural frequencies and normal modes via eigenvalue problems.