Anatomy of a Computer Simulation

Acomputer simulation is a way to study dynamical systems and their time evolution using digital calculations. This is of particular use in instances where theoretical calculations or experimental observations appear intractable, complicated, or downright impossible. For instance, the field of computer simulations, and more precisely, scientific simulations, arose out of a need to trace the explosion of nuclear weapons and to study weather evolution. Neither of these are topics one can hope to study via pen-and-paper analysis or through practical experiment (for obvious reasons).

Over the past few years, science has changed dramatically. As researchers, our desire for knowledge used to be sated one of two ways: theoretical modelling, and experimental observation. Both of these approaches left gaps in scientific discourse that simulations have come to fill. Rapid advances in computer technology have also spoiled us; we now seek faster & faster ways to find more & more accurate solutions to increasingly complex problems. Computer simulations are quite common today and serve almost all branches of the scientific community. To the layman, it would appear that simulations are a magical product of much wizardry. But as you will realise through this article, our lives as researchers have become much simpler today.

So What is a Computer Simulation

In the simplest possible terms, a computer simulation utilizes known physical laws and governing equations to trace out how a system would behave given certain initial conditions. In a majority of scientific investigations, this boils down to something as simple as solving Newton’s equation of motion, i.e., F = ma. While not without its limitations, this familiar relation still suffices as the basis of a plethora of studies. And while it may appear rather preposterous that something as simple as a three-variable equation would require giant supercomputers, remember that Newton’s equation is quite dubious. To the untrained eye, it is an expression for force. But consider the example of a pendulum, swinging about its pivot under the effect of gravity, giving us a tiny, harmless-looking equation,

The equation of motion of a pendulum

Do you remember what the solution looks like? Surely, it is just a combination of sine and cosine, right? Not even close. In fact, the solution can’t even be written in the form of elementary functions¹. It involves elliptic integrals and looks something like so:

Funnily enough, despite the scary looking symbols, the solution is actually quite elegant and looks like the image below:

Evolution of angular velocity versus angular displacement for different starting points (represented by the dots).

How to create a simulation

Knowing that computers can make light work of the most complicated of problems, let us now create our own simulation from scratch. Typically, such a job has three parts:

1. The input: this part sets up your problem and tells the computer what it needs to do and, to a certain extent, how it must do it. More often than not, this can be done via python.
2. The brains: the actual moving part which is typically hidden from the common user. It consists of functions that perform the calculation. Most scientific programs use lower-level languages for this part like C++ or Fortran, but python will do just fine for a simple tutorial.
3. The output: involves sifting through, and plotting the data obtained. Python comes to our rescue yet again.

Defining the problem

The nonlinear pendulum problem I showed just now felt like an inelegant example to understand how computer simulations work. However, I will choose to work with something far simpler: the dynamics of a charged particle in uniform electric and magnetic fields. The problem of a nonlinear pendulum, I leave to you as an exercise.

For this article, our problem statement is as follows:

Trace the trajectory of a charged particle in uniform electric and magnetic fields placed at 90⁰ to each other.

The Input Part

The input portion of your simulation has containers to hold information for various physical quantities. It also sets up parameters for the simulation. The following self-explanatory code block serves as a typical input part.

import numpy as np
import matplotlib.pyplot as plt

# We define the controlling parameters here. dt is our time-step, T our final time, Nt, number of time-steps.
# times variable holds the time-steps at which our positions and velocities will be calculated.
dt = 0.01
T  = 20
Nt = int(T/dt) + 1
times = np.linspace(0, T, Nt)

# We define our electric and magnetic fields.
electricField = np.array([0, 5e-2, 0.0])
magneticField = np.array([0, 0.0,  1.0])

# The characteristics of the particle to be tracked go here.
# Note that the velocity we are using is quite unrealisitic but it serves demonstration purpose.
q = 1
m = 1
v0 = np.array([1, 1, 1])
r0 = np.array([0.0,0.0,0.0])

initialState = np.hstack((r0,v0))

The Brains

Charged particles within electromagnetic fields follow the Lorentz force law, i.e., F = q(v × B), where the bold symbols represent vector quantities. The idea behind solving this numerically is to represent the force, which is just mass times the derivative of velocity, in a discrete manner². From thereon, one can use established numerical methods to solve the equation of motion and find velocity and position information at all times. Luckily for us, scipy ships with library functions to perform this step for us.

from scipy.integrate import solve_ivp

# The following function returns the derivative of a matrix that represents the position and velocity of the particle at a given time.
def lorentz(time, current):
  x, y, z, vx, vy, vz = current
  v = np.array([vx, vy, vz])
  a = (q / m) * (electricField + np.cross(v, magneticField))
  return [
        vx,
        vy,
        vz,
        a[0],
        a[1],
        a[2]
  ]

solution = solve_ivp(
    lorentz, (times[0],times[-1]),
    initialState,
    t_eval=times
)

The Postprocessing

Yes, it’s that easy! The solve_ivp method of the scipy library does the heavy lifting while all you have to do is plot the solutions obtained. To that end, here is a snippet to plot just the XY configuration space:

x = solution.y[0]
y = solution.y[1]

plt.plot(x, y)
plt.show()

XY configuration space plotted bare.

Of course, the no frills plot is quite ugly and conveys no real information whatsoever. Luckily, both the matplotlib as well as the seaborn library ship with options to make our plots vastly better. The following snippet adds labels and title to our plot and also colours the trajectory with time.

# Use seaborn to beautify the plot
from seaborn import set_theme, color_palette, cubehelix_palette
set_theme(style='ticks', context='notebook', palette='dark')
cubehelix = cubehelix_palette(as_cmap=True, reverse=True)

# Create figure
fig = plt.figure(figsize=(7,5))
plt.grid(True, which='both', linestyle=':')

plt.scatter(x,y, c=times, cmap=cubehelix)
plt.aspect('equal')

plt.title('XY Configuration Space')
plt.xlabel('x')
plt.ylabel('y')

fig.colorbar(label='Time')

plt.subplots_adjust(right=0.933)
plt.savefig('xy.png')

XY configuration space plotted with seaborn styling and axes labels. So much better!

What Next?

Nothing much, and everything! You can now try to repeat the exercise we just did for the case of a pendulum. You will need to modify the input and post-processing part, while the brains of the operation remain quite similar. Do try it out and comment on what you got.

Comments

In actual research, the “brains” are often highly parallelized programs developed by teams of insanely talented developers. For instance, the group I work with uses a software called SMILEI to probe the world of laser-plasma interactions. The process with using such software can be summarised as:

1. Create input script as per the software’s guidelines.
2. Use high-performance computers (the so-called supercomputers) to run the simulation with the input script provided. This requires submitting what is known as a “job” to the supercomputer.
3. Use various helper functions often provided by the software itself to extract raw data from the simulation output and plot what you need.

The possibilities then are endless and here is just a snippet of what we can create:

Footnotes

¹ The simple pendulum that we study in secondary school does have a sine-cosine solution since we approximate the angle to be quite small, which lets us put sin θ ≅ θ.

² This is called discretization. Basically, since a derivative is simply written as (f(x+h)-f(x)) / h in the limit of h going to 0, one can set h as the time-step and find f(x+h) from f(x) by utilizing the known-form of df/dt. More information here.