Euler, Runge-Kutta and Dynamic Simulation

Time integration of physical equations of motion using Euler, Verlet, and fourth-order Runge-Kutta schemes.

1. Formulation of Dynamical Systems and Ordinary Differential Equations

In theoretical and computational physics, the temporal evolution of a physical system's state is commonly described by a system of first-order ordinary differential equations (ODEs):

$\frac{d\vec{y}}{dt} = \vec{f}(t, \vec{y}), \quad \vec{y}(t_0) = \vec{y}_0$

where $\vec{y}(t) \in \mathbb{R}^N$ represents the state vector (such as positions $\vec{r}$, velocities $\vec{v}$, currents, or charges) and $\vec{f}: \mathbb{R} \times \mathbb{R}^N \to \mathbb{R}^N$ is the vector field defining the equations of motion. Any higher-order ODE, such as Newton's second law $m \frac{d^2 x}{dt^2} = F(t, x, \dot{x})$, is converted into a first-order system through auxiliary variables:

$\vec{y} = \begin{pmatrix} x \\ v \end{pmatrix}, \quad \frac{d\vec{y}}{dt} = \begin{pmatrix} v \\ \frac{1}{m} F(t, x, v) \end{pmatrix}$

The objective of numerical time integration is to construct a discrete trajectory $\{\vec{y}_n\}$ approximating the true analytical solution $\vec{y}(t_n)$ at discrete times $t_n = t_0 + n h$, where $h = \Delta t$ is the integration step size.

2. Explicit Euler Method and Its Dynamical Limits

The Forward Euler method is derived by truncating the Taylor expansion of $\vec{y}(t+h)$ after the linear term:

$\vec{y}(t+h) = \vec{y}(t) + h \frac{d\vec{y}}{dt} + \frac{h^2}{2} \frac{d^2\vec{y}}{dt^2} + \mathcal{O}(h^3)$

Dropping second-order and higher terms yields the update rule:

$\vec{y}_{n+1} = \vec{y}_n + h \vec{f}(t_n, \vec{y}_n)$

Critical Properties and Weaknesses:

1. Local Truncation Error (LTE): Each single step introduces an error of order $\mathcal{O}(h^2)$.
2. Global Truncation Error (GTE): Over a fixed interval $T = N h$, accumulated errors scale as $\mathcal{O}(h)$, making Euler a first-order method.
3. Phase Space Volume Distortion (Non-Symplectic Nature): For a simple harmonic oscillator $\dot{x} = v$, $\dot{v} = -\omega^2 x$, the step transition matrix is: $\begin{pmatrix} x_{n+1} \\ v_{n+1} \end{pmatrix} = \begin{pmatrix} 1 & h \\ -\omega^2 h & 1 \end{pmatrix} \begin{pmatrix} x_n \\ v_n \end{pmatrix}$ The determinant of this matrix is $\det(\mathbf{M}) = 1 + \omega^2 h^2 > 1$. Because the determinant exceeds unity, phase-space volume continually expands, injecting spurious, unphysical energy into the system at every timestep.

3. Symplectic Integration: Semi-Implicit Euler and Velocity-Verlet

For conservative Hamiltonian systems where $\mathcal{H} = T(\vec{p}) + V(\vec{q})$, symplectic integrators preserve phase-space volume (Liouville's theorem) and conserve an effective shadow Hamiltonian over extremely long integration times.

Symplectic Euler (Euler-Cromer)

Updates velocity first, then immediately uses the new velocity to advance position:

$v_{n+1} = v_n + h \, a(t_n, x_n)$ $x_{n+1} = x_n + h \, v_{n+1}$

The transition determinant satisfies $\det(\mathbf{M}) \equiv 1$, bounding energy fluctuations and preventing secular outward spiral drift.

Velocity-Verlet Algorithm

Widely adopted in molecular dynamics:

$x_{n+1} = x_n + v_n h + \frac{1}{2} a_n h^2$ $v_{n+1} = v_n + \frac{1}{2} (a_n + a_{n+1}) h$

It is strictly time-reversible and exhibits local position accuracy of $\mathcal{O}(h^4)$, achieving $\mathcal{O}(h^2)$ global accuracy for conservative trajectories.

4. Higher-Order Runge-Kutta Methods

Runge-Kutta schemes achieve higher Taylor accuracy by probing the derivative vector field $\vec{f}(t, \vec{y})$ at carefully chosen intermediate points within $[t_n, t_{n+1}]$, avoiding the need for analytical higher derivatives.

Second-Order Runge-Kutta (RK2 / Midpoint Method)

$k_1 = \vec{f}(t_n, \vec{y}_n)$ $k_2 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}_n + \frac{h}{2} k_1\right)$ $\vec{y}_{n+1} = \vec{y}_n + h k_2$ LTE is $\mathcal{O}(h^3)$ and GTE is $\mathcal{O}(h^2)$.

Classical Fourth-Order Runge-Kutta (RK4)

The gold standard for non-stiff ordinary differential equations in physics:

$k_1 = \vec{f}(t_n, \vec{y}_n)$ $k_2 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}_n + \frac{h}{2} k_1\right)$ $k_3 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}_n + \frac{h}{2} k_2\right)$ $k_4 = \vec{f}\left(t_n + h, \vec{y}_n + h k_3\right)$ $\vec{y}_{n+1} = \vec{y}_n + \frac{h}{6}\left(k_1 + 2k_2 + 2k_3 + k_4\right)$

• Local Truncation Error: $\mathcal{O}(h^5)$.
• Global Truncation Error: $\mathcal{O}(h^4)$. Halving the step size ($h \to h/2$) reduces global truncation error by $2^4 = 16$.

5. Adaptive Step Size Control (RK45)

In physical systems exhibiting multiple timescales (such as high-eccentricity planetary orbits with intense accelerations near periapsis), a static step size is either computationally wasteful or numerically unstable.

Embedded Runge-Kutta pairs (e.g., RKF45 or Dormand-Prince DOPRI5) calculate two estimates: an order-4 solution $\vec{y}_{n+1}$ and an order-5 solution $\hat{\vec{y}}_{n+1}$ using shared intermediate stages $k_i$. The local error estimate is:

$\epsilon = \|\hat{\vec{y}}_{n+1} - \vec{y}_{n+1}\|$

If $\epsilon > \text{TOL}$, the step is discarded and repeated with a smaller $h$. If $\epsilon \le \text{TOL}$, the step is committed and the next step is predicted via:

$h_{\text{new}} = S \cdot h \left(\frac{\text{TOL}}{\epsilon}\right)^{1/5}$

where $S \approx 0.84 - 0.9$ is an empirical safety margin.

6. Stiffness and Numerical Stability

A system is considered stiff when internal timescales differ by several orders of magnitude (Jacobian eigenvalues with large negative real components where $\max |\text{Re}(\lambda_i)| / \min |\text{Re}(\lambda_i)| \gg 1$).

For explicit integrators such as RK4, stability dictates step size constraints:

$h < \frac{C}{|\lambda_{\max}|}$

even if the asymptotic macroscopic motion is gentle. Stiff physical problems require implicit methods (such as Radau, BDF, or Gauss-Legendre) to maintain unconditional stability.

Implementation Procedure: From Physical ODEs to RK4 Numerical Integration

Step-by-Step Methodological Protocol

Simulating the dynamic evolution of physical systems governed by ODEs requires a rigorous workflow ensuring numerical stability, bounded error, and conservation of physical laws.

Step 1: State-Space Vectorization and Order Reduction

1. Identify all generalized coordinates and their higher derivatives.
2. Define the consolidated state vector $\vec{y}(t) = [y_1, y_2, \dots, y_N]^T$.
3. Express all first derivatives as closed-form functions of $\vec{y}$ and $t$: $\frac{d\vec{y}}{dt} = \vec{f}(t, \vec{y})$

Step 2: Algorithmic Selection Criteria

• Non-stiff systems with general trajectories: Choose RK4 (fixed or adaptive step).
• Conservative Hamiltonian dynamics over extended epochs: Choose symplectic integrators (Velocity-Verlet or Euler-Cromer) to eliminate secular energy drift.
• Stiff systems with disparate relaxation times: Adopt implicit schemes (e.g., backward differentiation formulas - BDF).

Step 3: RK4 Algorithmic Execution Loop

For each successive interval of step size $h$:

1. Calculate the initial derivative at the base state: $k_1 = \vec{f}(t_n, \vec{y}_n)$
2. Step half-way into the interval along $k_1$ to sample the midpoint derivative: $\vec{y}^{(1)} = \vec{y}_n + \frac{h}{2} k_1, \quad k_2 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}^{(1)}\right)$
3. Step half-way along $k_2$ to refine the midpoint derivative: $\vec{y}^{(2)} = \vec{y}_n + \frac{h}{2} k_2, \quad k_3 = \vec{f}\left(t_n + \frac{h}{2}, \vec{y}^{(2)}\right)$
4. Advance a full step along $k_3$ to sample the endpoint derivative: $\vec{y}^{(3)} = \vec{y}_n + h k_3, \quad k_4 = \vec{f}(t_n + h, \vec{y}^{(3)})$
5. Form the weighted Simpson-type average: $\vec{y}_{n+1} = \vec{y}_n + \frac{h}{6}(k_1 + 2k_2 + 2k_3 + k_4)$
6. Advance clock $t_{n+1} = t_n + h$.

Step 4: Tracking Invariants and Convergence Verification

At regular diagnostic checkpoints:

• Mechanical Energy Conservation: Track relative drift: $\Delta E(t) = \frac{|E(t) - E(0)|}{|E(0)|}$
• Momenta and Symmetries: Verify that numerical dissipation or truncation does not breach angular or linear momentum conservation.
• Convergence Ratio: Halving step size $h \to h/2$ should suppress global residual errors by a factor close to 16.

Worked Examples: Oscillators and Numerical Integration Methods

Example 1: Harmonic Oscillator - Analytical vs Numerical Benchmarks

Consider an undamped mass-spring system with mass $m = 1\text{ kg}$ and spring constant $k = 1\text{ N/m}$, giving natural angular frequency $\omega = \sqrt{k/m} = 1\text{ rad/s}$. The state-space equations are: $\frac{dx}{dt} = v, \quad \frac{dv}{dt} = -x$ Initial conditions: $x(0) = 1.0\text{ m}$, $v(0) = 0.0\text{ m/s}$. The exact analytical solution is: $x(t) = \cos(t), \quad v(t) = -\sin(t), \quad E(t) = \frac{1}{2}(v^2 + x^2) = 0.5\text{ J}$

We advance the system by a single timestep $h = 0.1\text{ s}$ comparing Forward Euler against RK4:

A. Forward Euler Integration:

$x_1 = x_0 + h v_0 = 1.0 + (0.1)(0.0) = 1.000000\text{ m}$ $v_1 = v_0 + h (-x_0) = 0.0 + (0.1)(-1.0) = -0.100000\text{ m/s}$ Total energy at $t = 0.1\text{ s}$: $E_{\text{Euler}} = \frac{1}{2}[(-0.1)^2 + (1.0)^2] = \frac{1}{2}(0.01 + 1.0) = 0.505000\text{ J}$ Mechanical energy grows by a spurious 1.0% in just one timestep.

B. Fourth-Order Runge-Kutta (RK4):

Intermediate slope evaluations:

1. $k_1^{(x)} = v_0 = 0.0$, $\quad k_1^{(v)} = -x_0 = -1.0$

2. Midpoint trial states with $h/2 = 0.05$: $x^{(1)} = x_0 + 0.05 k_1^{(x)} = 1.0 + 0 = 1.0$ $v^{(1)} = v_0 + 0.05 k_1^{(v)} = 0 - 0.05 = -0.05$ $k_2^{(x)} = v^{(1)} = -0.05, \quad k_2^{(v)} = -x^{(1)} = -1.0$

3. Second midpoint trial states: $x^{(2)} = x_0 + 0.05 k_2^{(x)} = 1.0 + 0.05(-0.05) = 0.9975$ $v^{(2)} = v_0 + 0.05 k_2^{(v)} = 0.0 + 0.05(-1.0) = -0.05$ $k_3^{(x)} = v^{(2)} = -0.05, \quad k_3^{(v)} = -x^{(2)} = -0.9975$

4. Full-step trial states at $h = 0.1$: $x^{(3)} = x_0 + 0.1 k_3^{(x)} = 1.0 + 0.1(-0.05) = 0.9950$ $v^{(3)} = v_0 + 0.1 k_3^{(v)} = 0.0 + 0.1(-0.9975) = -0.09975$ $k_4^{(x)} = v^{(3)} = -0.09975, \quad k_4^{(v)} = -x^{(3)} = -0.9950$

5. Weighted RK4 assembly: $x_1 = 1.0 + \frac{0.1}{6}[0 + 2(-0.05) + 2(-0.05) + (-0.09975)] = 0.99500417\text{ m}$ $v_1 = 0.0 + \frac{0.1}{6}[-1.0 + 2(-1.0) + 2(-0.9975) + (-0.9950)] = -0.09983333\text{ m/s}$

Exact analytical solution at $t = 0.1\text{ s}$: $x(0.1) = \cos(0.1) = 0.995004165\text{ m}, \quad v(0.1) = -\sin(0.1) = -0.099833417\text{ m/s}$ Absolute position error in RK4 is $|0.99500417 - 0.995004165| \approx 5 \times 10^{-9}\text{ m}$, maintaining 8 significant figures of precision.

Example 2: Large-Angle Nonlinear Pendulum

Consider a physical pendulum with length $L = 1\text{ m}$ and gravitational acceleration $g = 9.8\text{ m/s}^2$: $\frac{d\theta}{dt} = \omega, \quad \frac{d\omega}{dt} = -\frac{g}{L} \sin\theta = -9.8\sin\theta$ Initial angle at horizontal release: $\theta(0) = \frac{\pi}{2}\text{ rad}$ ($90^\circ$), $\omega(0) = 0\text{ rad/s}$. Timestep: $h = 0.05\text{ s}$.

Evaluating RK4 stages:

• $k_1^{(\theta)} = \omega_0 = 0$

• $k_1^{(\omega)} = -9.8 \sin(\pi/2) = -9.8000$

• $\theta^{(1)} = \pi/2 + 0.025(0) = 1.570796$

• $\omega^{(1)} = 0 + 0.025(-9.8) = -0.245000$

• $k_2^{(\theta)} = -0.245000$

• $k_2^{(\omega)} = -9.8 \sin(1.570796) = -9.8000$

• $\theta^{(2)} = \pi/2 + 0.025(-0.245000) = 1.564671$

• $\omega^{(2)} = 0 + 0.025(-9.8000) = -0.245000$

• $k_3^{(\theta)} = -0.245000$

• $k_3^{(\omega)} = -9.8 \sin(1.564671) = -9.799015$

• $\theta^{(3)} = \pi/2 + 0.05(-0.245000) = 1.558546$

• $\omega^{(3)} = 0 + 0.05(-9.799015) = -0.489951$

• $k_4^{(\theta)} = -0.489951$

• $k_4^{(\omega)} = -9.8 \sin(1.558546) = -9.796062$

State update: $\theta_1 = \frac{\pi}{2} + \frac{0.05}{6}[0 + 2(-0.245) + 2(-0.245) + (-0.489951)] = 1.558547\text{ rad}$ $\omega_1 = 0 + \frac{0.05}{6}[-9.8 + 2(-9.8) + 2(-9.799015) + (-9.796062)] = -0.489951\text{ rad/s}$ The full non-perturbative trigonometric restoring torque is preserved with fourth-order temporal fidelity.

Common Pitfalls in Time Integration and Dynamical Simulation

1. Using Forward Euler for Long-Term Conservative Hamiltonian Dynamics

• Pitfall: Applying the explicit forward Euler scheme to integrate celestial orbits or molecular systems over thousands of oscillations.
• Consequence: Because the discrete state update matrix determinant exceeds unity ($\det(\mathbf{M}) = 1 + \omega^2 h^2 > 1$), total mechanical energy drifts upward exponentially, causing artificial outward orbital ejection.
• Remedy: Utilize symplectic integrators (Euler-Cromer or Velocity-Verlet) for Hamiltonian dynamics, or high-order RK4 for short trajectories.

2. Conflating Local Truncation Error with Global Truncation Error

• Pitfall: Assuming an algorithm with local truncation error $\mathcal{O}(h^p)$ will exhibit global trajectory deviation scaling as $\mathcal{O}(h^p)$.
• Consequence: Over total simulation duration $T$, integrating requires $N = T/h$ successive steps. Accumulating $N$ local error contributions of order $\mathcal{O}(h^{p+1})$ yields a global truncation error of $\mathcal{O}(h^p)$. In RK4, local error is $\mathcal{O}(h^5)$, but global trajectory error is $\mathcal{O}(h^4)$.
• Remedy: Perform rigorous convergence tests tracking global end-state error as step size $h$ is halved.

3. Overlooking System Stiffness

• Pitfall: Integrating coupled dynamics with rapid relaxation modes using explicit RK4 with a step size chosen purely based on slow macroscopic motion.
• Consequence: High-frequency Jacobian eigenvalues require $h < 2.8/|\lambda_{\max}|$ solely for numerical stability. Exceeding this boundary triggers catastrophic exponential blowup to infinity or NaN values.
• Remedy: Evaluate the Jacobian spectrum and employ implicit integrators (such as BDF or Radau) for stiff regimes.

4. Failing to Advance Time $t$ in Non-Autonomous Intermediate Stages

• Pitfall: For explicitly time-dependent vector fields $\vec{f}(t, \vec{y})$, evaluating intermediate stages $k_2, k_3, k_4$ using the initial base time $t_n$.
• Consequence: The convergence order of the RK4 integrator degrades from fourth-order down to second or first order.
• Remedy: Ensure the intermediate evaluations use the correct times: $k_2, k_3$ at $t_n + h/2$ and $k_4$ at $t_n + h$.

5. Neglecting Physical Invariants and Conservation Checks

• Pitfall: Declaring a simulation physically valid merely because the executable completes without runtime crashes or NaN outputs.
• Consequence: Unphysical secular drift in energy, linear momentum, or angular momentum corrupts the trajectory unnoticed.
• Remedy: Implement runtime tracking and telemetry plotting relative drift $\Delta E(t)/E(0)$ and verifying Noether conservation laws.

Physical Applications of Numerical Integration Methods

1. Celestial Mechanics and the Gravitational Three-Body Problem

In astrodynamics and planetary astrophysics, the trajectories of $N$ bodies under mutual Newtonian gravity have no general closed-form analytical solution: $m_i \frac{d^2 \vec{r}_i}{dt^2} = \sum_{j \ne i} \frac{G m_i m_j}{|\vec{r}_j - \vec{r}_i|^3} (\vec{r}_j - \vec{r}_i)$ Integration using adaptive RK4 or high-order symplectic algorithms enables calculating deep-space trajectories (such as gravity-assist flybys for Voyager and Cassini missions) and predicting chaotic resonances in the Kuiper and asteroid belts.

2. Molecular Dynamics and Statistical Ensembles

In computational materials science and biophysics, millions of interacting atoms evolve under classical force fields (Lennard-Jones, Coulombic): $m_i \ddot{\vec{r}}_i = -\vec{\nabla}_i \sum_{j < k} V(|\vec{r}_j - \vec{r}_k|)$ The Velocity-Verlet algorithm serves as the bedrock of molecular dynamics: its symplectic geometry and exact time-reversibility ensure strict energy conservation across microcanonical (NVE) ensembles over millions of femtosecond timesteps without unphysical numerical heating.

3. Chaotic Dynamics and Atmospheric Convection (The Lorenz Attractor)

Edward Lorenz's idealized system of thermal fluid convection couples three nonlinear state variables: $\frac{dx}{dt} = \sigma (y - x), \quad \frac{dy}{dt} = x(\rho - z) - y, \quad \frac{dz}{dt} = x y - \beta z$ High-precision RK4 integration uncovers the fractal geometry of the strange attractor and enables numerical estimation of maximal Lyapunov exponents, characterizing deterministic chaos and fundamental limits of meteorological predictability.

4. Nonlinear Circuits and the Van der Pol Oscillator

In electronic engineering and plasma physics, self-sustaining relaxation oscillators follow the Van der Pol formulation: $\ddot{x} - \mu (1 - x^2) \dot{x} + x = 0$ When the damping parameter $\mu \gg 1$, dynamics become severely stiff, necessitating adaptive-step RK or implicit numerical solvers to capture abrupt limit-cycle voltage collapses.

5. Relativistic Charged Particle Beam Tracking in Accelerators

In high-energy synchrotron accelerators (such as CERN's LHC), bunches of relativistic charged particles orbit magnetic rings for hours: $\frac{d\vec{p}}{dt} = q(\vec{E} + \vec{v} \times \vec{B})$ Symplectic tracking integrators preserve phase-space emittance across billions of turns, preventing false predictions of beam degradation induced by non-physical numerical dissipation.

Practice Problems: Dynamical Simulation and Numerical Time Integration

Problem 1: Radioactive Decay and Taylor Polynomial Scaling in RK4

A radioactive isotope decays according to: $\frac{dN}{dt} = -\lambda N, \quad \lambda = 0.5\text{ s}^{-1}$ with initial quantity $N(0) = 1000.0\text{ nuclei}$.

1. Analytically derive the single-step amplification factor $N_{n+1} = G(z) N_n$ with $z = -\lambda h$ under RK4.
2. For $h = 0.2\text{ s}$, compute $N_1$ after the first step and evaluate the relative error against the exact analytical benchmark $N(0.2) = 1000 e^{-0.1}$.

Step-by-Step Solution:

1. The RK4 slopes evaluate to:
   • $k_1 = -\lambda N_n$
   • $k_2 = -\lambda (N_n + \frac{h}{2} k_1) = -\lambda N_n (1 - \frac{\lambda h}{2})$
   • $k_3 = -\lambda (N_n + \frac{h}{2} k_2) = -\lambda N_n [1 - \frac{\lambda h}{2} + \frac{(\lambda h)^2}{4}]$
   • $k_4 = -\lambda (N_n + h k_3) = -\lambda N_n [1 - \lambda h + \frac{(\lambda h)^2}{2} - \frac{(\lambda h)^3}{4}]$ Combining into the RK4 update with $z = -\lambda h$: $G(z) = 1 + z + \frac{z^2}{2!} + \frac{z^3}{3!} + \frac{z^4}{4!}$ This matches the degree-4 Taylor expansion of $e^z$ term by term.
2. Given $\lambda = 0.5$ and $h = 0.2$, parameter $z = -0.1$: $G(-0.1) = 1 - 0.1 + \frac{0.01}{2} - \frac{0.001}{6} + \frac{0.0001}{24} = 0.90483750$ $N_1 = 1000 \times 0.90483750 = 904.83750$ Exact benchmark: $N(0.2) = 1000 e^{-0.1} = 904.83742$ Absolute error: $\Delta N = 8 \times 10^{-5}$ The relative error remains below $0.00001\%$.

Problem 2: Duffing Oscillator State Evolution

The undamped bistable Duffing oscillator obeys: $\ddot{x} - x + x^3 = 0$ Reduce this ODE to a two-dimensional first-order state system and calculate the RK4 slope vectors $k_1$ and $k_2$ for initial state $x_0 = 0.5\text{ m}$, $v_0 = 1.0\text{ m/s}$ with step size $h = 0.05\text{ s}$.

Step-by-Step Solution:

1. State system: $\vec{y} = \begin{pmatrix} x \\ v \end{pmatrix}$, with vector field: $\frac{dx}{dt} = v, \quad \frac{dv}{dt} = x - x^3$
2. Stage 1: $k_1^{(x)} = v_0 = 1.0$ $k_1^{(v)} = x_0 - x_0^3 = 0.5 - (0.5)^3 = 0.375$
3. Midpoint trial coordinates ($h/2 = 0.025$): $x^{(1)} = x_0 + 0.025 k_1^{(x)} = 0.5 + 0.025(1.0) = 0.525$ $v^{(1)} = v_0 + 0.025 k_1^{(v)} = 1.0 + 0.025(0.375) = 1.009375$
4. Stage 2: $k_2^{(x)} = v^{(1)} = 1.009375$ $k_2^{(v)} = x^{(1)} - (x^{(1)})^3 = 0.525 - (0.525)^3 = 0.380297$

Problem 3: Kapitza Inverted Pendulum and Timescale Constraints

A pendulum driven by rapid vertical base motion $y_b(t) = a \cos(\Omega t)$ satisfies: $\ddot{\theta} = -\left(\frac{g}{L} - \frac{a \Omega^2}{L}\cos(\Omega t)\right)\sin\theta$ For physical parameters $L = 0.2\text{ m}$, $g = 9.8\text{ m/s}^2$, $a = 0.01\text{ m}$, and driving frequency $\Omega = 100\pi\text{ rad/s}$ ($50\text{ Hz}$):

1. Identify the fastest dynamic timescale in the system.
2. Determine the maximum allowable integration step $h$ to ensure at least 20 steps per driving cycle.

Step-by-Step Solution:

1. Fast external excitation period: $T_{\text{fast}} = \frac{2\pi}{\Omega} = \frac{1}{50} = 0.02\text{ s}$
2. Preserving accurate resolution of the high-frequency inertial force requires: $h \le \frac{T_{\text{fast}}}{20} = \frac{0.02\text{ s}}{20} = 0.001\text{ s} = 1\text{ ms}$ Selecting $h = 0.02\text{ s}$ would alias the cosine modulation to a constant, completely eliminating the time-averaged effective potential well that stabilizes the inverted upright equilibrium.