| SOCR β« | BPAD Website β« | BPAD GitHub β« |
Biomedical Physics with Applications to Disease (BPAD)
Chapter 4: Complex Particle Systems
BPAD Team + SOCR
January 2024
Some simple physical systems permit reliable deterministic modeling, or static model-based inference. For instance, for a single mass attached to a spring, we can use Newtonβs second law (\(F = m\times a\)) to model the motion of the system. Similarly, for a single atom it is possible to use the quantum-mechanical equivalent of \(F = m\times a\), the SchrΓΆdinger equation whose wavefunction solution captures the dynamics of the periodic atomic motion. Fitting a linear model of an outcome using observed data represents yet another example of a static representation of intrinsically stochastic processes that generate the observations.
The situation of modeling and representing the dynamics of more complex systems involving many particles is much more difficult, and sometimes impossible all together. Consider the example of modeling the number of particles in a cubic millimeter of blood (\(1mm^3\)), which may contain a number of different constituents with different volume, concentration, and
For instance, the hemoglobin concentration in blood depends on many factors (e.g., hormonal levels) and is expected to be around \(15 g/dl\pm 3 g/dl\) representing about \(1.4\times 10^{15}\) hemoglobin molecules in \(1 mm^3\) of blood. Whereas the corresponding concentration and number of water molecules \(H_2O\) in \(1 mm^3\) of blood are \(1 g/cm^3 \pm 0.01g/cm^3\) and \(3.3 \times 10^{19}\), respectively. Of course, the number of oxygen atoms attached to each of the four iron sites in the hemoglobin molecule may vary widely depending on the location of the blood sample, oxygenated and deoxygeneted blood travels outward from and into the lungs, respectively.
Modeling a 3D translational motion of each particle may require three equations representing the the position in each of the basic directions. Suppose at time \(t\) the force acting on a particle is \(F(t)=(F_x(t), F_y(t), F_z(t)):\mathbb{R}^+ \to \mathbb{R}^3\). Then, using Newtonβs second law, \(F=a\times m\), and first-order Taylor series expansion on a time interval \([t, t + \Delta t]\), the velocity of the particle, \(v(t)=(v_x(t), v_y(t), v_z(t))\), changes according to the three equations
\[\left | \begin{aligned} v_x (t + \Delta t) & = v_x (t) + \frac{F_x \Delta t}{m} \\ v_y (t + \Delta t) & = v_y (t) + \frac{F_y \Delta t}{m} \\ v_z (t + \Delta t) & = v_z (t) + \frac{F_z \Delta t}{m} \\ \end{aligned} \right. .\]
The corresponding three equations representing the \(x\), \(y\), and \(z\) changes of the position of the particle, \({\bf{x}}(t)=({\bf{x}}_x(t), {\bf{x}}_y(t), {\bf{x}}_z(t))\), are expressed as a second order Taylor series (in terms of the velocity \(v(t)\) and the acceleration, \(a(t)=\frac{F(t)}{m}\))
\[\left | \begin{aligned} {\bf{x}}_x(t + \Delta t) &= {\bf{x}}_x(t) + v_x (t)\Delta t + \frac{F_x (t)(\Delta t)^2}{2m} \\ {\bf{y}}_x(t + \Delta t) &= {\bf{y}}_x(t) + v_y (t)\Delta t + \frac{F_y (t)(\Delta t)^2}{2m} \\ {\bf{z}}_x(t + \Delta t) &= {\bf{z}}_x(t) + v_z (t)\Delta t + \frac{F_z (t)(\Delta t)^2}{2m} \\ \end{aligned} \right. .\]
For small \(\Delta t\), the second term can be dropped and we can use again a first order approximation for the position \({\bf{x}}(t)\). Solving these equations for each particle requires a number of addition and multiplication operations, e.g., for \(10^{19}\) particles (\(H_2O\) molecules in \(1cm^3\) of blood), the system solution may require over \(10^{20}\) arithmetic operations per time interval. Modern GPU chips can perform less than a petaflop operations, i.e., (\(10^{15}\)) floating point operations per second (FLOPS). This suggests that solving the position calculations solely for the water molecules in a tiny blood sample may require \(10^5\) seconds (over a day).
Computational complexity is just one issue that makes complex system modeling challenging. The physics of the processes represent another limitation that sometimes leads to deterministic chaos, even in relatively simple systems. For instance, identically configured experimental systems that are initialized using slightly different starting conditions may become unstable and generate progressively divergent subsequent behaviors even over a brief period of time. This is also known as the butterfly effect.
The tracing, modeling, forecasting, or prediction of the behavior of such chaotic systems is impossible, and sometimes may not be necessary. In certain situations, the overall properties of a system are more interesting than the specific instantiations or results. For instance, we may be interested in the distribution of the results, e.g., reporting average number of molecules, overall variation of pressure, skewness of the concentration, kurtosis of the speed, and so on.
1 Double Pendulum
A double pendulum is a system of a pair of interconnected pendula attached at their ends of their arms. The dynamics of the double pendulum system chaotic and sensitivity to its initial conditions.
Letβs model this system parameterized by the angles \(\theta_1\) and \(\theta_2\) between each arm and the vertical axis. Both arms are identical and of equal length, \(l\), and mass, \(m\). Assume that the uniform gravity acceleration is \(g\), the arm moment of inertia is \(I = \frac{1}{12}ml^2\) about the armβs central point of mass, and the origin of the Cartesian coordinate system is the base point of suspension of the first pendulum arm. Then, the positions of the centers of mass of each arm, the points \(P_1=(x_1, y_1)\) and \(P_2=(x_2, y_2)\), can be expressed in terms of the parameters \(\theta_i\):
\[{\begin{aligned}x_{1}&={\frac {l}{2}}\sin \theta _{1} \\ y_{1}&=-{\frac {l}{2}}\cos \theta _{1} \\ x_{2}&=l\left(\sin \theta _{1}+{\tfrac {1}{2}}\sin \theta _{2}\right)\\y_{2}&=-l\left(\cos \theta _{1}+{\tfrac {1}{2}}\cos \theta _{2}\right) \end{aligned}}.\]
The system Lagrangian can be expressed as
\[\begin{aligned} L &= \underbrace{ \overbrace{\frac {1}{2}m\left(v_{1}^{2}+v_{2}^{2}\right)}^{linear\ kinetic\ energy\\ of\ the\ center\ of\ mass} +\overbrace{\frac {1}{2}I\left({{\dot {\theta }}_{1}}^{2}+{{\dot {\theta }}_{2}}^{2}\right)}^{rotational\ kinetic\ energy\ around\\ the\ center\ of\ mass\ for\ each\ arm}}_{Total\ kinetic\ energy}- \underbrace{mg\left(y_{1}+y_{2}\right)}_{potential\ energy}\\ &={\tfrac {1}{2}}m\left({{\dot {x}}_{1}}^{2}+{{\dot {y}}_{1}}^{2}+{{\dot {x}}_{2}}^{2}+{{\dot {y}}_{2}}^{2}\right)+{\tfrac {1}{2}}I\left({{\dot {\theta }}_{1}}^{2}+{{\dot {\theta }}_{2}}^{2}\right)-mg\left(y_{1}+y_{2}\right)\end{aligned}.\]
The time derivatives are denoted by \(f'(t)\equiv \dot {f}(t)\) for each time-dependent variable.
In terms of the bivariate parameterization (\(\theta_1,\theta_2\)) of the centers of mass positions of the two arms, the Lagrangian is:
\[L= \frac {1}{6}ml^{2}\left({{\dot {\theta }}_{2}}^{2}+4{{\dot {\theta }}_{1}}^{2}+3{{\dot {\theta }}_{1}}{{\dot {\theta }}_{2}}\cos(\theta _{1}-\theta _{2})\right)+ \frac {1}{2}mgl\left(3\cos \theta _{1}+\cos \theta _{2}\right).\]
To derive the equations of motion for the double pendulum, we will express the two generalized momenta in terms of the angle parameters. Note that the total energy (kinetic plus potential) is conserved, but the momenta are not conserved:
\[\begin{aligned} p_{\theta_{1}} &= {\frac {\partial L}{\partial {{\dot {\theta }}_{1}}}}={\tfrac {1}{6}}ml^{2}\left(8{{\dot {\theta }}_{1}}+3{{\dot {\theta }}_{2}}\cos(\theta _{1}-\theta _{2})\right)\\ p_{\theta _{2}}&={\frac {\partial L}{\partial {{\dot {\theta }}_{2}}}}={\tfrac {1}{6}}ml^{2}\left(2{{\dot {\theta }}_{2}}+3{{\dot {\theta }}_{1}}\cos(\theta _{1}-\theta _{2})\right)\end{aligned}.\]
Solving the above equations for \((\dot{\theta}_{1},\dot{\theta}_{2})\), we get:
\[\begin{aligned} {{\dot {\theta }}_{1}}&={\frac {6}{ml^{2}}}{\frac {2p_{\theta _{1}}-3\cos(\theta _{1}-\theta _{2})p_{\theta _{2}}}{16-9\cos ^{2}(\theta _{1}-\theta _{2})}} \\ {{\dot {\theta }}_{2}}&={\frac {6}{ml^{2}}}{\frac {8p_{\theta _{2}}-3\cos(\theta _{1}-\theta _{2})p_{\theta _{1}}}{16-9\cos ^{2}(\theta _{1}-\theta _{2})}} \end{aligned} .\]
Hence, the four double pendulum equations of motion showing the time evolution of the chaotic system given its current state can be expressed in terms of the time-varying parameters \(\theta_1\) and \(\theta_2\):
\[\begin{aligned} {{\dot {\theta }}_{1}}&={\frac {6}{ml^{2}}}{\frac {2p_{\theta _{1}}-3\cos(\theta _{1}-\theta _{2})p_{\theta _{2}}}{16-9\cos ^{2}(\theta _{1}-\theta _{2})}}\\ {{\dot {\theta }}_{2}}&={\frac {6}{ml^{2}}}{\frac {8p_{\theta _{2}}-3\cos(\theta _{1}-\theta _{2})p_{\theta _{1}}}{16-9\cos ^{2}(\theta _{1}-\theta _{2})}} \\ {{\dot {p}}_{\theta _{1}}}&={\frac {\partial L}{\partial \theta _{1}}}=-{\tfrac {1}{2}}ml^{2}\left({{\dot {\theta }}_{1}}{{\dot {\theta }}_{2}}\sin(\theta _{1}-\theta _{2})+3{\frac {g}{l}}\sin \theta _{1}\right)\\ {{\dot {p}}_{\theta _{2}}}&={\frac {\partial L}{\partial \theta _{2}}}=-{\tfrac {1}{2}}ml^{2}\left(-{{\dot {\theta }}_{1}}{{\dot {\theta }}_{2}}\sin(\theta _{1}-\theta _{2})+{\frac {g}{l}}\sin \theta _{2}\right) \end{aligned}.\]
The simulation below shows the path traversed by the end-point of the
second pendulum arm. The simulation below solves the differential
equations above by numerical integration using the method
odeintr::integrate_sys(). We can explore the effects of
different initial conditions and parameter settings on the system
dynamics, as well as, run the simulation multiple times with the same
initial conditions.
# Setup configuration
library(tidyverse)
library(plotly)
# 1. Constants and initial conditions
G <- 9.807 # acceleration due to gravity, in m/s^2
L1 <- 1.0 # length of pendulum 1 (m)
L2 <- 1.0 # length of pendulum 2 (m)
M1 <- 4.0 # mass of pendulum 1 (kg)
M2 <- 4.0 # mass of pendulum 2 (kg)
parms <- c(L1, L2, M1, M2, G)
# initial conditions
theta1 <- 30 # initial angle theta of pendulum 1 (degree)
w1 <- 5.0 # initial angular velocity of pendulum 1 (degrees per second)
theta2 <- 60 # initial angle theta of pendulum 2 (degree)
w2 <- -5.0 # initial angular velocity of pendulum 2 (degrees per second)
state <- c(theta1, w1, theta2, w2)*pi/180 #convert degree to radians
x1 <- L1*sin(state[1])
y1 <- -L1*cos(state[1])
x2 <- L2*sin(state[3]) + x1
y2 <- -L2*cos(state[3]) + y1
df <- tibble(x1, y1, x2, y2)
plot_ly(df) %>%
# Add arm of first pendulum
add_segments(xend=x1, yend=y1, x=0, y=0, type="scatter", mode="lines",
line = list(width = M1), name="Arm 1") %>%
# Attach the arm of the second pendulum
add_segments(xend=x2, yend=y2, x=x1, y=y1, type="scatter", mode="lines",
line = list(width = M1), name="Arm 2") %>%
# Add theta1
add_text(x = 0.1, y = -0.5, text = "\U1D703<sub>1</sub>", name="\U1D703<sub>1</sub>") %>%
# Add theta2
add_segments(xend=x1, yend=y1-1, x=x1, y=y1, type="scatter", mode="lines",
line = list(width = 1, color="black"), name="Theta2 Angle") %>%
add_text(x = 0.6, y = -1.3, text = "\U1D703<sub>2</sub>", name="\U1D703<sub>2</sub>") %>%
# Add anchor base point
add_trace(x=0, y=0, type="scatter", mode="markers",
marker = list(size = 20), name="Anchor Base") %>%
# add joint point
add_trace(x=x1, y=y1, type="scatter", mode="markers",
marker = list(size = 20), name="Joint Point") %>%
# add last (free) end-point
add_trace(x=x2, y=y2, type="scatter", mode="markers",
marker = list(size = 20), name="Free End Point") %>%
# add first center of mass point (arm 1)
add_trace(x=x1/2, y=y1/2, type="scatter", mode="markers",
marker = list(size = 10), name="Arm 1 Center of Mass") %>%
add_text(x = x1/2+0.4, y =y1/2, text = "(x<sub>1</sub>,y<sub>1</sub>)", name="(x<sub>1</sub>,y<sub>1</sub>)") %>%
# add second center of mass point (arm 2)
add_trace(x=(x1+x2)/2, y=(y1+y2)/2, type="scatter", mode="markers",
marker = list(size = 10), name="Arm 2 Center of Mass") %>%
add_text(x = (x1+x2)/2+0.4, y =(y1+y2)/2, text = "(x<sub>2</sub>,y<sub>2</sub>)", name="(x<sub>2</sub>,y<sub>2</sub>)") %>%
# Add annotation for gravity acceleration
add_annotations(ax = (x1+x2)/2, ay =(y1+y2)/2-0.1, axref = "x", ayref = "y",
arrowhead=2, text= "<i>m g</i>", showarrow=T,
x=(x1+x2)/2, y =(y1+y2)/2-1) %>%
# add rectangular field of view and layout the scene
layout(title = 'Double Pendulum - Initial State',
shapes = list(list(type = "rect",
line = list(color = "blue"),
opacity = 0.1,
x0 = "-2", x1 = "2", xref = "x",
y0 = -3, y1 = 1, yref = "y")),
legend = list(orientation = 'h'),
yaxis = list(scaleanchor = "x"))Numerical integration for solving the differential equations and animating the solutions over a period of 60 seconds (s).
# install.packages("gganimate")
library(tidyverse)
library(gganimate)
# 1. Constants and initial conditions
G <- 9.807 # acceleration due to gravity, in m/s^2
L1 <- 1.0 # length of pendulum 1 (m)
L2 <- 1.0 # length of pendulum 2 (m)
M1 <- 4.0 # mass of pendulum 1 (kg)
M2 <- 4.0 # mass of pendulum 2 (kg)
parms <- c(L1, L2, M1, M2, G)
# initial conditions
theta1 <- 20.0 # initial angle theta of pendulum 1 (degree)
w1 <- 5.0 # initial angular velocity of pendulum 1 (degrees per second)
theta2 <- 180.0 # initial angle theta of pendulum 2 (degree)
w2 <- -5.0 # initial angular velocity of pendulum 2 (degrees per second)
state <- c(theta1, w1, theta2, w2)*pi/180 #convert degree to radians
# Partial derivatives calculator
pdes <- function(state, t){
L1 <- parms[1]
L2 <- parms[2]
M1 <- parms[3]
M2 <- parms[4]
G <- parms[5]
dydx <- rep(0,length(state))
dydx[1] <- state[2]
del_theta <- state[3] - state[1]
den1 <- (M1 + M2)*L1 - M2*L1*cos(del_theta)*cos(del_theta)
dydx[2] <- (M2*L1*state[2]*state[3]*sin(del_theta)*cos(del_theta) +
M2*G*sin(state[3])*cos(del_theta) +
M2*L2*state[4]*state[4]*sin(del_theta) -
(M1 + M2)*G*sin(state[1]))/den1
dydx[3] <- state[4]
den2 <- (L2/L1)*den1
dydx[4] <- (-M2*L2*state[4]*state[4]*sin(del_theta)*cos(del_theta) +
(M1 + M2)*G*sin(state[1])*cos(del_theta) -
(M1 + M2)*L1*state[2]*state[2]*sin(del_theta) -
(M1 + M2)*G*sin(state[3]))/den2
return(dydx)
}
# 2. ODE Integration to find solutions of equations of motion
# start duration and step_size parameters control the time dynamics of the double pendulum
# We can run the pendulum for 60 seconds and record the end-point position at increments of 0.1 seconds.
# install.packages("odeintr")
library(odeintr)
ode_solver <- integrate_sys(pdes, init = state, duration = 60, start = 0, step_size = 0.1)
# 3. Use the theta parameters obtained in step 2 via integration
# compute the spatial coordinates for both pendulums (x_i,y_i)
x1 <- L1*sin(ode_solver[, 2])
y1 <- -L1*cos(ode_solver[, 2])
x2 <- L2*sin(ode_solver[, 4]) + x1
y2 <- -L2*cos(ode_solver[, 4]) + y1
df <- tibble(t = ode_solver[,1], x1, y1, x2, y2, group=1)
# 4. Plot the DF containing the spatial positions of the pendulum arms for time
# and animating the pendulum dynamics
y_offset <- 1.5 # shift upwards the base point (or set to zero)
# pl <- ggplot(df)+
# geom_segment(aes(xend=x1, yend=y1+y_offset), x=0, y=y_offset, size=M1, col="green")+
# geom_segment(aes(xend=x2, yend=y2+y_offset, x=x1, y=y1+y_offset, size=M2, col="red"), show.legend = FALSE)+
# geom_point(size=5, x=0, y=y_offset, col="purple")+
# geom_point(aes(x1, y1+y_offset), col="yellow", size=5)+
# geom_point(aes(x2,y2+y_offset), col="blue", size=5)+
# scale_y_continuous(limits=c(-2,2))+
# scale_x_continuous(limits=c(-2,2))+
# theme(legend.position="none") +
# ggraph::theme_graph()+
# labs(title="{frame_time} s")+
# transition_time(t)
# double_pendulum <- animate(pl, nframes=nrow(df), fps=25)
# double_pendulum
# 5. Trace the path of the end-point of the second pendulum
# install.packages("transformr")
library(transformr)
tmp <- select(df, t, x2, y2)
trail <- tibble(x=c(sapply(1:5,function(x) lag(tmp$x2,x))),
y=c(sapply(1:5,function(x) lag(tmp$y2,x))),
t=rep(tmp$t,5)) %>%
dplyr::filter(!is.na(x))
pl <- ggplot(df)+
geom_path(data=trail, aes(x,y+y_offset, group=t), colour="blue", size=0.5)+
geom_segment(aes(xend=x1, yend=y1+y_offset), x=0, y=y_offset, size=M1, col="green")+
geom_segment(aes(xend=x2, yend=y2+y_offset, x=x1, y=y1+y_offset, size=M2, col="red"), show.legend = FALSE)+
geom_point(size=5, x=0, y=y_offset, col="purple")+
geom_point(aes(x1, y1+y_offset), col="yellow", size=5)+
geom_point(aes(x2,y2+y_offset), col="blue", size=5)+
scale_y_continuous(limits=c(-2,2))+
scale_x_continuous(limits=c(-2,2))+
theme(legend.position="none") +
ggraph::theme_graph()+
labs(title="Time {frame_time} s")+
transition_time(t) +
shadow_mark(colour="grey", size=0.2, exclude_layer = 2:6)
double_pendulum <- animate(pl, nframes=nrow(df), fps=25) #, rewind=TRUE, start_pause=3)
double_pendulum
anim_save("C:/Users/Dinov/Desktop/DoublePendulumAnime.gif", double_pendulum)
# 