Mass-spring with RK2

This example is useful because the mass-spring oscillator is simple, but not trivial.
It has a known analytical solution, it is periodic, and it allows us to study both convergence and numerical stability.

1. Physical model

The system consists of a mass m attached to a spring with elastic constant k.
The motion takes place along one spatial direction, which we call the x axis.
The displacement x(t) is measured from the equilibrium position.

The model assumptions are:

• the motion is one-dimensional;
• the mass is treated as a point mass;
• the spring is ideal;
• there is no friction;
• there is no air resistance;
• the spring constant k and the mass m are constant.

The spring force is described by Hooke’s law:

The minus sign means that the spring force is a restoring force: it always points toward the equilibrium position.

Using Newton’s second law:

and substituting the spring force:

Since acceleration is the second derivative of position:

we obtain:

Dividing by m gives:

2. Differential equation of the oscillator

The equation of motion can be written as:

Introducing the natural angular frequency:

the equation becomes:

This is the equation of the harmonic oscillator.

3. First-order system

To solve the problem numerically, the second-order equation is transformed into a first-order system.
I introduce the velocity:

Then:

and:

Therefore the system to integrate numerically is:

The state vector is:

and the right-hand side of the ODE is:

4. Python model function

The following function implements the physical model. The convention used here is f(t,u),
which is the same convention used by many numerical ODE solvers.

import numpy as np

def f(t, u):
    """Physical model"""
    x, v = u

    dxdt = v
    dvdt = -(k / m) * x

    return [dxdt, dvdt]

5. RK2 midpoint method

The midpoint version of the second-order Runge-Kutta method is:

The idea is to first estimate the derivative at the beginning of the interval, then use it to reach a midpoint estimate.
The final update uses the derivative evaluated at that midpoint.

def rk2_midpoint_step(u, model, t, dt):
    """Second-order Runge-Kutta method with midpoint."""
    u = np.asarray(u, dtype=float)

    k1 = np.asarray(model(t, u), dtype=float)

    u_mid = u + (dt / 2) * k1
    t_mid = t + dt / 2

    k2 = np.asarray(model(t_mid, u_mid), dtype=float)

    return u + dt * k2


def integrate_rk2(u0, model, t0, t_end, dt):
    """Time integration loop."""
    n = int(np.ceil((t_end - t0) / dt))

    t = t0 + dt * np.arange(n + 1)
    u = np.zeros((n + 1,) + np.shape(u0), dtype=float)

    u[0] = u0

    for k in range(n):
        u[k + 1] = rk2_midpoint_step(u[k], model, t[k], dt)

    return t, u

6. Initial conditions and parameters

I use the simplest case: the mass starts displaced from equilibrium and is released from rest.

m = 1.0
k = 1.0

x0 = 1.0
v0 = 0.0

u0 = np.array([x0, v0], dtype=float)

t0 = 0.0
t_end = 20.0
dt = 0.05

7. RK2 midpoint trajectory

The first simulation uses only RK2 midpoint. With x_0=1 and v_0=0,
the position starts from its maximum value and then oscillates.

8. Explicit Euler for comparison

To understand the improvement introduced by RK2, I also implement explicit Euler.

Euler uses only the derivative at the current point. This makes it simple, but also less accurate.

def euler_step(u, model, t, dt):
    """Single explicit Euler step."""
    u = np.asarray(u, dtype=float)

    return u + dt * np.asarray(model(t, u), dtype=float)


def integrate_euler(u0, model, t0, t_end, dt):
    """Time integration loop with explicit Euler."""
    n = int(np.ceil((t_end - t0) / dt))

    t = t0 + dt * np.arange(n + 1)
    u = np.zeros((n + 1,) + np.shape(u0), dtype=float)

    u[0] = u0

    for k in range(n):
        u[k + 1] = euler_step(u[k], model, t[k], dt)

    return t, u

9. Comparing Euler and RK2

The discrepancy between the two methods is measured by:

The two position curves can also be plotted directly:

This comparison is useful, but it is not yet a true error analysis, because neither numerical method is the exact solution.
For a rigorous error estimate we compare both methods against the analytical solution.

10. Analytical solution

For the ideal harmonic oscillator, the analytical solution is known:

This allows us to compute the numerical error directly.

def analytical_solution(t, x0, v0, k, m):
    """Analytical solution of the mass-spring harmonic oscillator."""
    omega = np.sqrt(k / m)

    x_exact = x0 * np.cos(omega * t) + (v0 / omega) * np.sin(omega * t)
    v_exact = -x0 * omega * np.sin(omega * t) + v0 * np.cos(omega * t)

    return x_exact, v_exact

11. Convergence

To study convergence, I repeat the simulation for decreasing time steps:

For each time step I compute the maximum error on the position:

The order of convergence is estimated with:

For explicit Euler, we expect first-order convergence:

For RK2 midpoint, we expect second-order convergence:

dt_values = [0.2, 0.1, 0.05, 0.025, 0.0125]

errors_euler = []
errors_rk2 = []

for dt_test in dt_values:
    t_euler, u_euler = integrate_euler(u0, f, t0, t_end, dt_test)
    t_rk2, u_rk2 = integrate_rk2(u0, f, t0, t_end, dt_test)

    x_euler = u_euler[:, 0]
    x_rk2 = u_rk2[:, 0]

    x_exact_euler, _ = analytical_solution(t_euler, x0, v0, k, m)
    x_exact_rk2, _ = analytical_solution(t_rk2, x0, v0, k, m)

    error_euler = np.max(np.abs(x_euler - x_exact_euler))
    error_rk2 = np.max(np.abs(x_rk2 - x_exact_rk2))

    errors_euler.append(error_euler)
    errors_rk2.append(error_rk2)

12. Numerical stability through energy

For an ideal mass-spring oscillator without friction, the total mechanical energy should remain constant:

where:

is the kinetic energy, while:

is the elastic potential energy.

If a numerical method is stable and accurate, the normalized energy E(t)/E(0) should stay close to 1.

def mechanical_energy(u, m, k):
    """Total mechanical energy of the harmonic oscillator."""
    x = u[:, 0]
    v = u[:, 1]

    E = 0.5 * m * v**2 + 0.5 * k * x**2

    return E

The energy plot shows that explicit Euler introduces artificial energy into the system.
RK2 midpoint behaves better, although it is not a perfectly energy-conserving method.

13. Long-time stability

The difference becomes clearer over a longer time interval.

14. Conclusion

This example shows why RK2 midpoint is a clear improvement over explicit Euler.
Both methods solve the same ODE system, but they have different accuracy and stability behaviour.

The convergence test shows that explicit Euler behaves like a first-order method, while RK2 midpoint behaves like a second-order method.
The energy test shows that explicit Euler tends to inject artificial energy into the oscillator, while RK2 midpoint keeps the energy much closer to the correct constant value.

Therefore, the mass-spring oscillator is a good educational test case: it has a simple physical model, a known analytical solution, and it clearly highlights the difference between first-order and second-order numerical integration methods.