Building a leapfrog integrator

Building a leapfrog integrator#

Malena Rice (Yale University)

Description: Write a simple N-body integrator in Python using the leapfrog formulation.

Intended Audience: Intermediate undergraduate

tags: N-body integration,functions, modeling

Requirements: numpy, matplotlib

Last Updated: August 22, 2025

Learning Objectives

  1. Implement a simple physical model in Python from scratch.

  2. Demonstrate the importance of resolution settings in numerical simulations.

  3. Visualize the outputs of a model to provide intuition for a system’s behavior.

Introduction#

In this exercise, you will build a Python function that can solve the N-body problem in 2D as a function of time. The acceleration between two bodies in the N-body problem, decomposed into x- and y-components, is given as

\[a_x = \frac{Gm(dx)}{r^3}\]

$\(a_y = \frac{Gm(dy)}{r^3}\)$.

Taking the first step of your N-body integrator#

We will initialize our system with the Sun at coordinate (0,0) and with the Earth at coordinate (0, 1), with distance coordinates given in AU and with the Earth moving in the positive x-direction.

Exercise 1

  1. Initialize a system with the Sun at coordinate (0,0) and with the Earth at coordinate (0, 1), where distance coordinates are given in AU.

  2. Calculate the acceleration between the Sun and the Earth at this starting point. Demonstrate that the two masses are accelerated toward each other by gravity, and that there is no x-component to this initial acceleration.

To run our code, we will need to set some initial conditions. To keep the center of mass static in the 2-body problem, we need to conserve momentum between our two bodies:

\[m_1 v_1 + m_2 v_2 = 0\]

Exercise 2

  1. Based on momentum conservation, calculate the velocity of the Sun around the Earth, in km/s, for a typical Earth orbital speed of 29.8 km/s.

The general “kick-drift-kick” form of a leapfrog integrator is formulated as follows:

\[v_{i+1/2} = v_i + a_i \frac{\Delta t}{2}\]
\[x_{i+1} = x_i + v_{i + 1/2} \Delta t \]
\[v_{i+1} = v_{i+1/2} + a_{i+1}\frac{\Delta t}{2} \]

This set of three functions corresponds to one “step” forward in the leapfrog integration routine.

Exercise 3

  1. Implement these lines of code within a for loop in Python as follows:

    for i in range(N):

    first line of code to advance to \(v_{i+1/2}\)
    second line of code to advance to \(x_{i+1}\)
    third line of code to calculate the new acceleration at \(x_{i+1}\)
    fourth line of code to advance to \(v_{i+1}\)

  2. To check that your code works, start at the same initial point described in Question 1. Advance your code a single step forward, using a time step of 1 day. Print the position and velocity that you obtain after 1 time step. Is the particle moving in the direction that you expect, given the initial conditions and the acceleration from Question 1?

Hint

Hint: You will need to make sure that your inputs have consistent units.

Integrating an orbit#

At this point, you have most of the setup for your numerical integration. The last point that you need to implement is a tracker: a way to save your values when you take more than one time step.

Exercise 4

  1. Create arrays that will track your positions and velocities at each time step.

Everything should now be in place to run your leapfrog integrator. In the next exercise, we will run the integrator and visualize its outputs to demonstrate that the integrator is functioning properly.

Exercise 5

  1. Integrate your program forward in time 5 years, using 100,000 evenly spaced time steps.

  2. Plot how the positions of the Sun and the Earth change over this time.

Sufficiently high resolution in numerical simulations is critical to ensure that the relevant physical processes are captured. In this final exercise, we will show how our integrator’s outputs depend upon the timestep selected.

Exercise 6

  1. Implement the same algorithm, but adjust your time steps so that you only use 50 evenly spaced steps over the course of 5 years. How does your result differ from your previous integration using 100,000 time steps? Why does this result look so different?

Contents