For today’s recreational coding exercise, we will simulate a spring network: an array of nodes connected by springs that obey Hooke’s Law. The spring system will fall down due to gravity and bounce off the floor. We use a particle algorithm similar to our N-body simulation.

You may find the accompanying Python code on github.

Before we begin, below is a gif of what running our simulation looks like:

Animation of our spring network simulation
Animation of our spring network simulation

Force Calculation

We will assume a system of N particles (nodes), indexed by i=1,…,N. Each node has a:

  • velocity v = [ vx, vyᵢ, vzᵢ


For today’s recreational coding exercise, we learn how to simulate dilute gas with the Direct Simulation Monte Carlo (DSMC) method. This approach is useful when gas is weakly collisional. We will look at the Rayleigh Problem: the response of rarified gas between two plates after one suddenly starts moving.

You may find the accompanying Python code on github.

But first, below is a gif of what running our simulation looks like:

Rarefied Gas

Rarefied flows cannot be accurately simulated using fluid (Navier-Stokes) equations. By rarified or ‘dilute’, we mean that the mean free path λ of a molecule is of the same…


For today’s recreational coding exercise, we simulate active matter, i.e., swarming. Such a system may describe a flock of birds or a school of fish. We will look at how very simple rules may lead to the emergence of self-ordered motions.

You may find the accompanying Python code on github.

But first, below is a gif of what running our simulation looks like:

Viscek model for flocking behavior

We will describe a famous minimal model for active matter called the Viscek model (1995). Despite its simplicity, the model is able to describe universal properties of swarming behavior.

The model consists of N moving particles, indexed…


For today’s recreational coding exercise, we simulate fluid flow past a cylinder using the Lattice Boltzmann method. This is a really cool and simple technique to simulate fluid flow: instead of evolving the fluid (Navier-Stokes) equations directly, microscopic particles on a lattice are simulated with streaming and collision processes. The power of the method comes from reducing the high-dimensionality of the microscopic physics onto the lattice, which has limited degrees of freedom.

You may find the accompanying Python code on github.

But first, below is a gif of what running our simulation looks like:

Fluid dynamics on a Lattice

We will begin with a microscopic…


For today’s recreational coding exercise, we look at a simple way to create volume renderings to visualize 3D simulation datacubes. This technique is incredibly useful when you have space-filling data you would like to visualize. Such data shows up often in astrophysical datasets but also in other areas of computer graphics and medical data (CT scans and MRIs).

You may find the accompanying Python code on github.

Before diving in, below is a gif of what running the rendering algorithm on some data may look like:

Datacube

We assume that the data is an Nx × Ny × Nz datacube of…


For today’s recreational coding exercise, we will look at quantum mechanical systems, in particular, the Schrodinger-Poisson equations. We will create a simulation for the evolution of a wavefunction under its self-potential. Such a system may describe certain superfluids/Bose-Einstein condensate or exotic dark matter.

You may find the accompanying Python code on github.

Before diving in, below is a gif of what running our simulation looks like:

Schrodinger-Poisson Equations

We will consider the evolution of the wavefunction ψ(x) under the Schrodinger equation:


For today’s recreational coding exercise, we will simulate the Kelvin-Helmholtz Instability with the Finite Volume method. We will consider a compressible fluid with a high density stream moving in opposite direction of the background. The velocity shear induces a famous instability that is seen sometimes in clouds as well as Jupiter’s Great Red Spot.

You may find the accompanying Python code on github.

Before we begin, below is a gif of what running our simulation looks like:

Kelvin-Helmholtz Instability
Kelvin-Helmholtz Instability

Finite Volume Method

We will describe the finite volume method to simulate an ideal compressible fluid. Extensions of the method exist for the simulation of other…


For today’s recreational coding exercise, we will simulate a star with smoothed-particle hydrodynamics (SPH). We will start with some initial condition and relax it into a stable stellar structure and measure the star’s density as a function of radius. The SPH method has many applications in astrophysics and elsewhere, and is in fact a general one to simulate all types of fluids.

You may find the accompanying Python code on github.

Before we begin, below is a gif of what running our simulation looks like:

SPH simulation
SPH simulation

SPH Method

We will represent the fluid as a collection of N point particles, indexed by i=1,…,N


For today’s recreational coding exercise, we will investigate plasma physics with particle-in-cell (PIC) simulations. We will create a simulation of two electron beams passing through each other in opposite directions. An interesting phenomenon occurs, called the Two-Stream Instability.

You may find the accompanying Python code on github.

Before we begin, below is a gif of what running our simulation looks like:

Plasma PIC simulation of the Two-Stream Instability
Plasma PIC simulation of the Two-Stream Instability

Particle-In-Cell (PIC) Method

We will consider a one-dimensional system of electrons in an unmagnetized uniform medium of ions. The electrons will be described by N particles, indexed by i, each having a

  • velocity vᵢ


For today’s recreational coding exercise, we will look at the gravitational N-body problem. We will create a simulation of a dynamical system of particles interacting with each other gravitationally. Such a system may describe the orbits of planets in the Solar System or stars in our Galaxy.

You may find the accompanying Python code on github. (And if you prefer to use Matlab, please see my Matlab version of the article)

But first, below is a gif of what running our simulation looks like:

N-body simulation
N-body simulation

Force Calculation

We will assume a system of N point particles, indexed by i=1,…,N. Each particle has a:

Philip Mocz

Computational Astrophysicist @Princeton, sharing intro tutorials on creating your own computer simulations! Harvard ’12 (A.B), ’17 (PhD). Connect with me @PMocz

Get the Medium app

A button that says 'Download on the App Store', and if clicked it will lead you to the iOS App store
A button that says 'Get it on, Google Play', and if clicked it will lead you to the Google Play store