# What are good particle dynamics ODEs for an introductory scientific computing course?

I'm teaching an introductory course on scientific computing (programming in C/C++) and am looking for application problems which the assignments can be centered around. I'm thinking of ODEs for particle dynamics $$x_i'(t) = u(x_i, y_i, t), \quad y_i'(t) = v(x_i, y_i, t)$$ where particle positions $$x_i(t), y_i(t)$$ are solved using a simple time-stepper given a background velocity field $$(u,v)$$.

I'd like to make the system more coupled/complex to give students a more realistic picture of a scientific computing application, but without introducing more complex mathematical concepts like PDE discretizations. Are there good ways to do so within particle dynamics (for example, making the velocity particle-dependent)?

• While not necessarily "particle" based, chemical reaction kinetics problems provide some interesting ODE dynamics, especially if you try modeling a chemical oscillator. RLC circuits are also pretty interesting, and building up the system matrix introduces some finite element-like concepts without any spatial discretization information. Jul 21, 2022 at 16:23
• Just out of curiosity, if you are interested in the approach others have taken, I really recommend checking out Nick Trefethen's course on scientific computing – podcasts.ox.ac.uk/series/scientific-computing-dphil-students. It is all algorithm focused rather than computing focussed and although he uses matlab as a "toy" language, the dynamics and structures of the algorithms discussed really shine through. He is a SIAM von Neumann Prize winner and I think has a good pedagogical style. (I also really, really liked his course when I took it several decades ago!) Jul 22, 2022 at 13:12
• Thanks! Our Numerical Analysis is very similar to this, and we definitely use several of Trefethen's examples liberally (it was designed by a close collaborator of Trefethen back in the 00s) Jul 22, 2022 at 16:00
• I trust your motives but I was curious why do an into course in C/C++? I guess it would depend on what level this course is at if you can mention that (undergrad or grad and what CS background they have). Jul 24, 2022 at 22:13
• It's an upper level undergrad or early grad, and the C/C++ has been a part of the course for ages (the Julia part is new). I'd like a Julia-focused course, but C/C++ are probably more common tools in most scientific computing research groups (and industry) so we keep them. Jul 25, 2022 at 23:11

For a single particle, interesting dynamics already arise if you are in magnetic and electric fields.

The situation becomes even more interesting if you consider several particles at once and how they interact. An example is the solar system, where particles interact gravitationally. But if you consider charged particles, you can also consider electromagnetic interactions between charged particles through the electric and magnetic fields their respective motion generates.

• Imho. the example with the solar system is very good. You can start by simulating two bodies (sun/earth), and verify our 24h cycle. When that works you can extend to three bodies and blissfully observe the chaos:-) It is also easy to see the importance of numerical precision in a nonlinear problem as the solutions will diverge fast if the initial conditions are altered slightly. One can also track if the numerical method preserves some integral quantities (angular momentum, kinetic+potential energy). Jul 21, 2022 at 6:03
• That's a nice idea! A gravitational system with lots of bodies might be too chaotic, but a small 2 or 3-body example could be used to construct a particle simulation like youtube.com/watch?v=2n7nFjuS98I. My goal was to find an application where increasing the number of particles results in a fluid-like or "emergent" behavior; I think this would work. Thanks! Jul 21, 2022 at 6:39
• When you have the code for the 2/3 body gravitational system, you can easily switch the interacting force to a Lennard-Jones potential which has an equilibrium distance. Then you give the bodies the mass of atoms and build tiny crystals that wiggle and have temperature (oh ... the memories!). As a side effect, you can also teach how to plot the results from .csv via matplotlib/gnuplot etc. Jul 22, 2022 at 8:38
• This is a great idea. A switch to magnetohydrodynamics – moving particles creating a current and thus a B field – is a fairly straightforward extension that produces a horrible set of DEs out the other end that are wonderfully susceptible to rounding errors. Jul 22, 2022 at 13:08
• Dealing with rounding errors in a careful way is a bit outside the scope of this course IMO. The gravitational system is "large enough" computationally and implementationally to give students a sense of how to structure a reasonably efficient computational science project; additional numerical or physical aspects probably belong in a different course. Jul 22, 2022 at 15:54

One very educational example is the Lodka-Volterra system. It can describe the observed effects of predator and prey population levels in many ecological system (foxes & rabbits). High populance of foxes will reduce the number of rabbits and so on :

from Wikipedia:

\begin{align} \frac{dx}{dt} &= \alpha x - \beta x y, \\ \frac{dy}{dt} &= \delta x y - \gamma y, \end{align} where:

• $$x$$ is the number of prey

• $$y$$ is the number of some predator

• $$\frac{dy}{dt}$$ and $$\tfrac{dx}{dt}$$ represent the instantaneous growth rates of the two populations;

• $$t$$ represents time

• $$α , \beta, \gamma, \delta$$ are positive real parameters describing the interaction of the two species.

Depending on the parameters you may observe different outcomes. It is a nice introduction into the modelling of our ecosystem, and the setting is easy to imagine. As a student, it made me realize just how complex nature is. Every forest or landscape is dominated by the interaction of millions of different species of plants/animals/bacteria etc. (Even if you add a third species to the model, it may turn into a chaotic system).

other applications:

https://en.wikipedia.org/wiki/Pork_cycle

economic applications

The solar system has already been proposed, so I have to fall back onto my next suggestion ! Although it's further away from pure particle dynamics, you may also consider a simple model of a rope, as a series of point masses (your particles) which are linked sequentially by springs. Numerically, this can be challenging to integrate with classical explicit integrators if the spring stiffness is too high though.

You may produce nice animations with that. Here is one I made (in the case of infinitely stiff springs, I will spare the details) :

(see my GitHub)

EDIT: you can even do that in 2D to make a crude model of an elastic material, where a grid of material points are connected with springs (both horizontally and vertically). The problem remains easy to compute (linear and sparse), and allows for arbtrarily high number of points to be used, which may be useful if you want to explore basic code performance tweaks.

• That's a nice example too! As a first step, doing the double pendulum might also be a good start and it, too, has nice dynamics. Jul 21, 2022 at 15:08
• That is a great example! Jul 24, 2022 at 22:09

The Lorenz System, which has its origins in fluid (thermo)dynamics is a good introduction to first-order coupled ODEs with simple specification but non-trivial dynamics (i.e solutions that are chaotic in the face of small adjustments in the choice of initial conditions, or in the parameters $$\sigma$$, $$\rho$$ and $$\beta$$)

The system is specified in terms of arbitrary variables $$x$$ $$y$$ and $$z$$ by the following equations

\begin{align} \frac{\mathrm{d}x}{\mathrm{d}t} &= \sigma (y - x), \\[6pt] \frac{\mathrm{d}y}{\mathrm{d}t} &= x (\rho - z) - y, \\[6pt] \frac{\mathrm{d}z}{\mathrm{d}t} &= x y - \beta z. \end{align}

The trajectory $$(x,y,z)(t)$$ of the solution draws out the "chaos butterfly" image popular in elementary introductions to chaos theory.

(Equations and diagram from https://en.wikipedia.org/wiki/Lorenz_system)

• This is also a very well known system. However it may be harder to envision a more complex coupled system. This seems to exist (I looked up coupled Lorenz systems), but the dynamics and physical meaning of the system might be more complex to grasp for students. Jul 21, 2022 at 18:11

While the pandemic is ravaging, also consider the SIR model which describes an infectious desease spreading through a population. The population is separated into three groups (variables): Susceptible->Infected->(Recovered/dead/immune). At the beginning the number of susceptible persons is near 100%, the more susceptible and not already infected people there are, the faster the virus will spread. Over time, as more people have had the honor, the group of recovered/dead/immune people will grow. It is highly relevant to our time, and there are a lot of extensions that could serve as homework:-)

$$\frac {\mathrm{d}S}{\mathrm{d}t}= \nu \, N - \beta \frac {S \, I} {N} - \mu \, S\;$$

$$\frac{\mathrm{d}I}{\mathrm{d}t}= \beta \frac {S \, I} {N} - \gamma \, I - \mu \, I\;$$

$$\frac {\mathrm{d}R}{\mathrm{d}t}= \gamma \, I - \mu \, R\;$$

https://de.wikipedia.org/wiki/SIR-Modell

• It's an important model, but imo the dynamics are not so interesting on its own. It can be extended with some traveling or community components. Jul 22, 2022 at 9:47
• I do like this model too; I was also considering the zombie version of the SIR model with travel between cities as a possible app. Jul 22, 2022 at 15:55
• I have done a simple SIR type model regarding rabid foxes although it is a simple PDE. You can get the gist of the SIR model but depending on your parameters getting traveling wave solutions. This does sound something similar to your zombie problem (rabid foxes move randomly in space). Jul 24, 2022 at 22:11
• @KyleMandli do you have a link to a derivation? That sounds like a neat system Aug 21, 2022 at 17:46
• I just looked but could not find what I had done I am afraid. I will try to respond with a Jupyter notebook when I create one as I need to teach a course that will require it anyway. Aug 22, 2022 at 23:51

# Suggestion: Double pendulum

All of the above colleagues had great suggestions. In this theme, my proposal would be the double pendulum problem. This system is a classic, simple and complex at the same time, its equations of motion can be obtained by applying the theories of Lagrange, Hamilton or Newton of mechanics, the latter and the least used in the development of the equations of motion of the double pendulum. With a little depth, we can enter chaos theory, and use this system to study this physical phenomenon, because the system has a high sensitive dependence on the initial conditions, and it is very common for programmers to make simulations of two or more pendulums with slight differences in the initial conditions, precisely to observe the different trajectories of the pendulums in the simulations. In the figure below we have the equations of the double pendulum developed with the help of Hamilton's theory of mechanics.

And as one of the results, the following graph of the trajectory curves of the respective masses $$m_{1}$$ and $$m_{2}$$ in the plan x and y can be generated.

# numeric methods

The equations of motion are composed of a set of ordinary differential equations, with initial conditions imposed on the system, basically it is a problem like $$\frac{d\vec y}{dt} = \vec f(t,y)$$ with $$t_{0} = \alpha$$ and $$\vec y(t_{0}) = \vec y_{0}$$.

There are many numerical methods used to solve this type of problem numerically, but the numerical methods of the Runge-Kutta family are widely used in this particular problem, mainly the fourth order Runge-Kutta (RK4). For reasons of curiosity, I have already made implementations in third, fifth, sixth order for comparison between the different methods and plotting the errors committed for comparison between them. In this article, the author demonstrates an explicit sixth order Runge-Kutta formula, which you could be implementing and comparing with the classical fourth order.

If you use the python language, in this case it is better to use the integrators already developed and implemented in scipy.integrate, such as scipy.integrate.odeint or scipy.integrate.solve_ivp, among other options.

# A Python code example of the double pendulum problem

import numpy as np
from scipy.integrate import odeint, RK23, RK45, solve_ivp
import matplotlib.pyplot as plt

m1 = 1.00
m2 = 1.00
l1 = 1.00
l2 = 1.00
g  = 9.81

t0    = 0.0
po1_0 = 0.0
po2_0 = 0.0

t_maximo = 60
dt       = 1.0e-3
N        = 100000

# Hamiltoniana
'''
def H(o1,o2,po1,po2):
return ((m2*l2**2*po1**2+(m1+m2)*l1**2*po2**2-2*m2*l1*l2*po1*po2*np.cos(o1-o2))/(2*m2*l1**1*l2**2*(m1+m2*(np.sin(o1-o2))**2))) - (m1+m2)*g*l1*np.cos(o1)-m2*g*l2*np.cos(o2)
'''
H = lambda o1,o2,po1,po2:\
((m2*l2**2*po1**2+(m1+m2)*l1**2*po2**2-2*m2*l1*l2*po1*po2*np.cos(o1-o2))/(2*m2*l1**1*l2**2*(m1+m2*(np.sin(o1-o2))**2))) - (m1+m2)*g*l1*np.cos(o1)-m2*g*l2*np.cos(o2)

def pd(s,t):
o1,o2,po1,po2 = s
h1   = ((po1*po2*np.sin(o1-o2))/(l1*l2*(m1+m2*(np.sin(o1-o2))**2)))
h2   = ((m2*l2**2*po2**2+(m1+m2)*l1**2*po2**2-2*m1*l1*l2*po1*po2*np.cos(o1-o2))/(2*l1**2*l2**2*(m1+m2*(np.sin(o1-o2))**2)**2))
edo1 = (l2*po1-l1*po2*np.cos(o1-o2))/(l1**2*l2*(m1+m2*(np.sin(o1-o2))**2))
edo2 = ((-m2*l2*po1*np.cos(o1-o2)+(m1+m2)*l1*po2)/(m2*l1*l2**2*(m1+m2*(np.sin(o1-o2))**2)))
edo3 = -(m1+m2)*g*l1*np.sin(o1)-h1+h2*np.sin(2*o1-2*o2)
edo4 = -m2*g*l2*np.sin(o2)+h1-h2*np.sin(2*o1-2*o2)
return np.array([edo1,edo2,edo3,edo4])

t  = np.arange(t0, t_maximo, dt) # Range para o tempo, de t0 á t_maximo com o passo dt
s0 = np.array([o1_0, o2_0, po1_0, po2_0], dtype=np.float64) # Vetor de condições iniciais para os ângulos e momentos canônicos conjugados aos ângulos

# Método scipy.integrate.odeint(f, s0, t0)
s = odeint(pd,s0,t)

o1, o2   = s[:,0], s[:,1] # Ângulos theta1 e theta2 (radianos)
po1, po2 = s[:,2], s[:,3] # Momentos canônicos conujulgados aos ângulos (radianos/s)

x1 = l1*np.sin(o1)
x2 = l1*np.sin(o1)+l2*np.sin(o2)
y1 = -l1*np.cos(o1)
y2 = -l1*np.cos(o1)-l2*np.cos(o2)

plt.figure()
plt.plot(x1,y1,'g.',x2,y2,'b.',0,0,'yx',linewidth = 0.1)
plt.xlabel("x(m)")
plt.ylabel("y(m)")
plt.grid()
plt.title(r'Posições das massas $$m_{1}(kg)$$ e $$m_{2}(kg)$$')
plt.show()