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

Question

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)?

Answer 1

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.

Answer 2

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

Answer 3

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.

Answer 4

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)

Answer 5

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

Answer 6

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.

• Article: An Explicit Sixth-Order Runge-Kutta Formula By H. A. Luther - https://www.ams.org/journals/mcom/1968-22-102/S0025-5718-68-99876-1/S0025-5718-68-99876-1.pdf

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
o1_0  = np.radians(60)
o2_0  = np.radians(60)
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 ((m2l22po12+(m1+m2)*l12po22-2m2l1l2po1po2np.cos(o1-o2))/(2m2l11l22(m1+m2(np.sin(o1-o2))2))) - (m1+m2)gl1np.cos(o1)-m2gl2np.cos(o2)
'''
H = lambda o1,o2,po1,po2:

    ((m2l22po12+(m1+m2)*l12po22-2m2l1l2po1po2np.cos(o1-o2))/(2m2l11l22(m1+m2(np.sin(o1-o2))2))) - (m1+m2)gl1np.cos(o1)-m2gl2np.cos(o2)
def pd(s,t):
    o1,o2,po1,po2 = s
    h1   = ((po1po2np.sin(o1-o2))/(l1l2(m1+m2(np.sin(o1-o2))2)))
    h2   = ((m2l22*po22+(m1+m2)l12po22-2m1l1l2po1po2np.cos(o1-o2))/(2*l12l22(m1+m2(np.sin(o1-o2))2)2))
    edo1 = (l2po1-l1po2np.cos(o1-o2))/(l12l2(m1+m2*(np.sin(o1-o2))2))
    edo2 = ((-m2l2po1np.cos(o1-o2)+(m1+m2)l1po2)/(m2l1l22(m1+m2(np.sin(o1-o2))2)))
    edo3 = -(m1+m2)gl1np.sin(o1)-h1+h2np.sin(2o1-2o2)
    edo4 = -m2gl2np.sin(o2)+h1-h2np.sin(2o1-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 = l1np.sin(o1)
x2 = l1np.sin(o1)+l2np.sin(o2)
y1 = -l1np.cos(o1)
y2 = -l1np.cos(o1)-l2np.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()