To demonstrate the power and flexibility of numerical multi‑degree‑of‑freedom (MDOF) analysis, we use a non‑classical damping matrix and compare the direct time‑history solution with the classical-damping approximation. This could be the case of a base‑isolated building (high damping in the isolator, lower damping in the superstructure), or a structure with story‑level supplemental dampers that do not combine into a Rayleigh form. The same Newmark formulation handles both without modal decomposition.
This page walks through a small 2‑DOF shear-frame example: we report $\mathbf{M}$, $\mathbf{K}$, and $\mathbf{C}$, apply a short sustained load, and show how the coupled equations are integrated step-by-step. The result is a compact reference for the practising engineer who wants to see the algorithm in action.
The model corresponds to a two‑story shear building—for example a short office or laboratory with one lateral degree of freedom per floor. Lumped masses and story stiffnesses are $m_1=m_2=1$, $k_1=20$, $k_2=10$ (units arbitrary but consistent). Mass and stiffness matrices are
Damping is taken as diagonal in physical coordinates, $\mathbf{C}=\operatorname{diag}(c_1,c_2)$, with different effective damping per story—e.g. story 1 at about $2\%$ and story 2 at about $22\%$— as might occur with a soft isolator or a story with supplemental dampers. No single Rayleigh form $\alpha\mathbf{M}+\beta\mathbf{K}$ can match both; the system is strongly non‑classically damped:
The structure is subjected to a sustained unit load at the roof over a short window: $p_2(t)=1$ for $0\le t\le 0.5$ s and $p_2(t)=0$ afterwards, with $p_1(t)\equiv 0$. The plot below shows this loading.
Because $\mathbf{C}$ is diagonal in physical coordinates and not of the form $\alpha \mathbf{M}+\beta \mathbf{K}$, the system is non‑classically damped: the undamped mode shapes do not diagonalize $\mathbf{C}$. This is exactly the situation where direct MDOF time‑history shines—no need to force a Rayleigh fit or to use complex modal analysis. To show the contrast, we still run a modal time‑history using approximate modal damping ratios $\zeta_n$ from the diagonal of $\boldsymbol{\Phi}^T \mathbf{C} \boldsymbol{\Phi}$ (mass‑normalized $\boldsymbol{\Phi}$). Those $\zeta_n$ are only approximate; the reference response is the direct MDOF result.
The undamped modal properties come from $\mathbf{K}\boldsymbol{\phi}_n = \omega_n^2 \mathbf{M}\boldsymbol{\phi}_n$ with mass normalization $\boldsymbol{\phi}_n^T \mathbf{M}\boldsymbol{\phi}_n = 1$. The approximate modal damping ratio is
Because $\mathbf{C}$ is not exactly classical, the true modal equations would contain velocity-coupling terms; the modal time‑history here drops those terms so we can see the impact of the classical-damping approximation.
For the MDOF system $\mathbf{M}\ddot{\mathbf{u}} + \mathbf{C}\dot{\mathbf{u}} + \mathbf{K}\mathbf{u} = \mathbf{p}(t)$, the constant-average-acceleration Newmark scheme with parameters $(\beta,\gamma)=(1/4,1/2)$ defines a set of predictor coefficients $a_0,\dots,a_5$ in terms of the time step $\Delta t$. These are the standard Newmark coefficients used to construct the effective stiffness and effective load at each step.
The effective stiffness matrix and effective load vector at step $i$ are
At each time step, the algorithm solves the coupled 2×2 system $\mathbf{K}_\text{eff}\Delta\mathbf{u}^1 = \mathbf{p}_\text{eff}^1$ to obtain $\mathbf{u}^1$. The velocities and accelerations are then updated from the Newmark formulas.
Below is the core of the algorithm in real Python (constant-average-acceleration: $\beta=1/4$, $\gamma=1/2$). The effective stiffness $\mathbf{K}_\text{eff}$ is formed once; each step solves one linear system and updates displacement, velocity, and acceleration. No modal decomposition is required—so it works for any $\mathbf{C}$.
# Newmark coefficients (β=1/4, γ=1/2)
a0 = 1.0 / (beta * dt**2)
a1 = gamma / (beta * dt)
a2 = 1.0 / (beta * dt)
a3 = 1.0 / (2*beta) - 1.0
a4 = gamma/beta - 1.0
a5 = dt * (gamma/(2*beta) - 1.0)
K_eff = K + a0*M + a1*C # formed once
# Initial acceleration from M ü = p - C u̇ - K u at t=0
udd[0] = np.linalg.solve(M, p_t[0] - C @ ud[0] - K @ u[0])
for i in range(1, n_steps):
p_eff = (p_t[i]
+ M @ (a0*u[i-1] + a2*ud[i-1] + a3*udd[i-1])
+ C @ (a1*u[i-1] + a4*ud[i-1] + a5*udd[i-1]))
u[i] = np.linalg.solve(K_eff, p_eff)
udd[i] = a0*(u[i] - u[i-1]) - a2*ud[i-1] - a3*udd[i-1]
ud[i] = ud[i-1] + dt*((1-gamma)*udd[i-1] + gamma*udd[i])With non‑classical damping, the direct MDOF solution in physical coordinates is the reference; the modal time‑history is an approximation that assumes classical damping (diagonal $\zeta_n$ in modal space). Any difference between the two curves comes from the neglected velocity-coupling terms in the modal equations—a reminder of why numerical MDOF is so useful when damping is not Rayleigh.
For this example the gap is small but visible: the modal approximation captures the main response, while fine details differ. In practice, for base‑isolated or damper‑enhanced structures, running the direct MDOF formulation (as above) gives the engineer a single, consistent reference without modal approximations.
This demo underpins the CE223 isolation dashboard: there we assume classical damping and use modal superposition; here we relax that assumption and show that the same Newmark engine handles both cases—the difference in results is due to the damping model, not the time‑integration method.