Analytical solutions are elegant. Real life is messy.

In previous chapters, we enjoyed the luxury of closed-form solutions like $x(t) = X \sin(\omega t)$. But nature rarely gives us clean sine waves. What happens when your bridge is hit by a random earthquake? What if your spring stiffness changes as it compresses (non-linear)? What if the damping isn’t constant?

In these cases, the “perfect” equations of motion become unsolvable by hand.

In this post, based on Chapter 11 of Rao’s Mechanical Vibrations, we trade the calculus textbook for the compiler. We introduce Numerical Integration Methods—the computational algorithms that slice time into tiny steps to solve even the most challenging, chaotic vibration problems.

The Numerical Approach: Discrete Time

Analytical methods try to satisfy the differential equation of motion at every instant of time (like a continuous movie film). Numerical methods, by contrast, satisfy the equation only at discrete time points (like a slide show).

We divide the total duration $T$ into $n$ equal time steps $\Delta t$: $$ t_0=0, \quad t_1=\Delta t, \quad t_2=2\Delta t, \dots $$

Our goal is to jump from one stepping stone ($t_i$) to the next ($t_{i+1}$). Numerical methods generally fall into two categories:

1. Explicit Integration (Looking Back): We calculate the future state $x_{i+1}$ using only the past information ($x_i, x_{i-1}$, etc.). It’s like driving by looking in the rearview mirror.

   • Examples: Finite Difference, Runge-Kutta.
   • Pros: Simple, fast per step.
   • Cons: Can be unstable if $\Delta t$ is too large.

2. Implicit Integration (Looking Forward): The future state $x_{i+1}$ depends on itself (via the equation of motion at $t_{i+1}$). This requires solving a system of simultaneous equations at every step.

   • Examples: Newmark-$\beta$, Wilson-$\theta$.
   • Pros: Unconditionally stable (great for stiff problems).
   • Cons: Computationally expensive per step.

The Finite Difference Method

How do we represent “velocity” and “acceleration” without calculus? We approximate them using simple algebra. This is the core of the Finite Difference Method.

We start with the Taylor Series Expansion, which predicts the future ($x_{i+1}$) and the past ($x_{i-1}$) based on the current state ($x_i$):

$$ x_{i+1} = x_i + \Delta t \dot{x}_i + \frac{(\Delta t)^2}{2}\ddot{x}_i + \dots $$

$$ x_{i-1} = x_i - \Delta t \dot{x}_i + \frac{(\Delta t)^2}{2}\ddot{x}_i - \dots $$

Central Difference Formulas

How do we isolate velocity? If we subtract the second equation from the first, the $x_i$ and $\ddot{x}_i$ terms cancel out, leaving us with a clean expression for velocity.

Similarly, if we add the two equations, the $\dot{x}_i$ terms cancel out, allowing us to solve for acceleration.

These are the Central Difference formulas. They are called “central” because they look at both the left ($i-1$) and right ($i+1$) neighbors to estimate the value at the center ($i$).

The Math (Simplified):

• Velocity ($\dot{x}$): The average slope between the point before and the point after. $$ \dot{x}_i = \frac{x_{i+1} - x_{i-1}}{2\Delta t} $$
• Acceleration ($\ddot{x}$): The “curvature” or rate of change of the slope. $$ \ddot{x}_i = \frac{x_{i+1} - 2x_i + x_{i-1}}{(\Delta t)^2} $$

Solving the Single-Degree-of-Freedom System

To solve the equation of motion:

$$ m\ddot{x}_i + c\dot{x}_i + kx_i = F_i $$

We don’t solve it as a differential equation. Instead, we substitute our algebraic approximations for $\ddot{x}_i$ and $\dot{x}_i$ directly into the equation.

This turns the calculus problem into a pure algebra problem. We then group all the $x_{i+1}$ terms on one side and the known $x_i, x_{i-1}$ terms on the other.

The result is a Recurrence Formula. This is a “recipe” that tells us how to calculate the next position $x_{i+1}$ using only the current ($x_i$) and past ($x_{i-1}$) positions:

$$ x_{i+1} = \frac{F_i - (k - A)x_i - (C - B)x_{i-1}}{A + C} $$

(Where A, B, C are constants based on mass $m$, damping $c$, and time step $\Delta t$).

The “Starting” Problem

Look at the diagram above. To calculate the very first step $x_1$ (at $t_1$), the formula asks for $x_0$ (which we know) and $x_{-1}$ (which doesn’t exist!).

Since the method isn’t self-starting, we must create a fictitious “Ghost Point” $x_{-1}$ by working backwards from the initial velocity $v_0$ and acceleration $a_0$: $$ x_{-1} = x_0 - \Delta t \dot{x}_0 + \frac{(\Delta t)^2}{2} \ddot{x}_0 $$

Stability

Stability: Why does it blow up?

The Central Difference method is conditionally stable. This isn’t just a mathematical quirk; it has a physical meaning.

If your time step $\Delta t$ is too large, the numerical approximation adds “artificial energy” to the system faster than the physical damping can dissipate it. The simulation acts like a system with negative damping, causing the amplitude to grow uncontrollably to infinity.

To stay safe, the time step must be smaller than a specific fraction of the system’s natural period $\tau_n$:

$$ \Delta t < \frac{\tau_n}{\pi} \approx 0.318 \tau_n $$

Other Numerical Methods

While the Central Difference method is simple and intuitive, it’s not the only game in town. Two other heavyweights dominate the field:

The Scout: Runge-Kutta (RK4)

• Concept: Instead of blindly jumping forward based on one slope, RK4 sends out 4 “scouts” to test the slope at different points in the interval. It then takes a weighted average to make a highly accurate move.
• Best For: General-purpose problems where high accuracy is needed (e.g., orbital mechanics, long-duration simulations).
• Type: Explicit (Conditionally Stable).

The Anchor: Implicit Methods (Newmark-$\beta$, Wilson-$\theta$)

• Concept: Explicit methods calculate the future based on the equilibrium at the current time $t_i$ (which is known). Implicit methods are more demanding: they insist that the equation of motion must be satisfied at the future time $t_{i+1}$.
• The Trade-off: Since we don’t know $x_{i+1}$ yet, this creates a complex system of equations that must be solved simultaneously at every step. But the reward is unconditional stability.
• Best For: Large-scale structural dynamics (e.g., earthquakes on skyscrapers). You can take massive time steps without crashing the simulation.
• Type: Implicit (Unconditionally Stable).

Julia Example: Central Difference Method

Modern dynamics often relies on Julia for its speed and clean syntax. Here is the complete implementation of the Central Difference Method for an SDOF system.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33

# Central Difference Method for SDOF System in Julia
using Plots

# 1. System Parameters
m, c, k = 1.0, 0.2, 10.0  
F0, ω_f = 5.0, 2.0

# 2. Time Settings
dt = 0.01     
t = 0:dt:20
n = length(t)
x = zeros(n)
F = F0 .* cos.(ω_f .* t)

# 3. Initialization
val_0, vel_0 = 0.0, 0.0
acc_0 = (F[1] - c*vel_0 - k*val_0) / m

# 4. Starting Step (Ghost Point / Taylor Series)
x[1] = val_0
x[2] = val_0 + vel_0*dt + 0.5*acc_0*dt^2

# 5. Loop (Constants pre-calculated for speed)
A = m + c*dt/2
B = k*dt^2 - 2*m
C = m - c*dt/2

for i in 2:n-1
    x[i+1] = (F[i]*dt^2 - B*x[i] - C*x[i-1]) / A
end

plot(t, x, title="SDOF Response", label="x(t)", lw=2)
savefig("central_difference_result.svg")

Summary

We have crossed the bridge from the elegance of exact calculus to the brute force of numerical arithmetic. While we lose perfect continuity, we gain the power to solve any system—no matter how complex the force or how non-linear the physics.

Here is your toolkit for stepping through time:

What’s Next?

Our single-degree-of-freedom journey is complete. But real machines aren’t just one mass and one spring; they are continuums.

In MV-12, we will scale up. We will apply these numerical tools to Finite Element Analysis (FEM), discretizing not just time, but space itself, to simulate the vibration of complex, real-world structures.

References

Rao, S. S. (2018). Mechanical Vibrations (6th ed.). Pearson.