In this chapter we focus on developing the particle dynamics code. First we install required packages:
import Pkg; Pkg.add(["Plots", "Printf", "Zygote"])
using LinearAlgebra
In many day-to-day tasks there is not yet the one canonical way to do it in Julia. This also concerns plotting and visualisation of data. Right now quite a few different packages with their own APIs and feature sets exist. To give you some examples:
What we will use in this course is Plots. I like it because it is a metapackage unifying some of the above under a common interface. Flipping the backend is just a single function call, so you can change between them as needed. Also it will probably not disappear any time soon ;). We load Plots
first and explicitly switch to the default GR
backend.
using Plots
Plots.gr()
Plotting is then done with the plot
function. Additional data series can be added using plot!
and attributes can be changed via keywords arguments or standalone functions like xaxis!
, title!
and so on.
Vho(r, a=0.5) = (r - a)^2 # Shifted harmonic oscillator
dVho(r, a=0.5) = 2r -2a # Derivative
r = collect(-5:0.1:5) # Create series of x-values
p = plot(r, Vho.(r), label="Vho a=0.5")
plot!(p, r, dVho.(r), label="∂Vho")
xaxis!(p, "relative radial distance")
Of course scatter plots can be made as well ...
scatter(r, Vho.(r), label="Vho", marker=:cross)
scatter!(r, dVho.(r), label="∂Vho", marker=:circle)
xaxis!("relative radial distance")
... and in fact there are plenty of things you can do to style your plot further, see attributes and layouts for some ideas if you are curious.
Recall that for our dynamics with a given potential, say $$ V\left( \begin{array}{c} x_1 \\ x_2 \end{array}\right) = x_1^2 + x_2^2 $$ we needed the acceleration $$ \vec{A}_V = - \nabla V. $$
For the harmonic oscillator this is kind of ok, but for more complicated potentials $V$, taking the derivative is error prone and most importantly quite boring. The solution is the great package Zygote
, which is capable of adjoint-mode automatic differentiation. Without going into details, this means that the Julia compiler takes the derivatives for us and the result is (usually almost) as fast as if we did it ourselves. Let's try this in 1D on Vho
as defined above:
using Zygote
r = collect(-5:0.1:5)
p = plot(r, Vho.(r), label="Vho a=0.5")
xaxis!(p, "relative radial distance")
plot!(p, r, Vho'.(r), label="∂Vho") # Notice the dash
With ease we add the second derivative to the plot:
plot!(p, r, Vho''.(r), label="∂∂Vho")
This is extremely helpful for more complicated potentials, for example the Morse potential, which is a common model for a chemical bond: $$ V_\text{Morse} = D (1 - exp(-\alpha (x-x_0)))^2 - D$$
Vmorse(r; re=0.7, α=1.3, D=10) = D*(1 - exp(-α * (r - re)))^2 - D
r = collect(0:0.1:7)
p = plot(r, Vmorse.(r), label="Vmorse", ylim=(-25, 25))
xaxis!(p, "Radial distance")
plot!(p, r, Vmorse'.(r), label="Vmorse'")
plot!(p, r, Vmorse''.(r), label="Vmorse''")
By the way: This works for higher dimensions and more complicated expressions, too, we will use this in a sec.
Now we return to our forward-euler implementation, which works for both 1D and higher dimensions:
euler(A, Δt, xn, vn) = xn + vn * Δt, vn + A(xn) * Δt
where $A_V$ (or in the code $A$) was the negative gradient of a potential $V$. With plotting at hand, let us actually see the dynamics. We define a plot function to plot the potential and animate the particle over time:
using Printf
# The arguments after the ; are so-called keyword arguments, they can be omitted
# and then the default value after the = is used
function plot_dynamics_line(V, Δt, n_steps; x0=0, v0=randn(), integrator=euler,
xlim=(-5, 5), ylim=(0, 10))
t, x, v = 0, x0, v0 # Initialise state
# Compute potential values (for plotting)
xrange = xlim[1]:0.1:xlim[2]
Vrange = V.(xrange)
# @gif is a macro to "record" an animation of the dynamics,
# each loop iteration gives a frame
@gif for i in 1:n_steps
# Propagate dynamics (notice the derivative)
x, v = integrator(x -> -V'(x), Δt, x, v)
t += Δt
# Plot potential
p = plot(xrange, Vrange, label="Potential", xlim=xlim, ylim=ylim)
# Plot the particle and its velocity (as an arrow)
timestr = @sprintf "%.2f" t # Format time as a nice string
scatter!(p, [x], [V(x)], label="Particle at t=$timestr")
plot!(p, [(x, V(x)), (x + 0.5v, V(x))], arrow=1.0,
label="particle velocity / 2")
end
end
# Now let's actually see it ....
plot_dynamics_line(Vho, 0.1, 200)
Wow ... that's strange. Our plot points at the well-known problem that energy is not conserved for a simple forward Euler scheme. We can also diagnose this using a phase-space plot, where the time evolution of particle position $x$ and particle velocity $v$ is shown. This should be a closed curve if energy is conserved, so let's see ...
function plot_phase_space(A, Δt, n_steps; x0=0, v0=randn(), integrator=euler,
xlim=(-7, 7), ylim=(-7, 7))
x, v = x0, v0
p = plot([x], [v], xlim=xlim, ylim=ylim, label="") # Plot first position
xaxis!(p, "position")
yaxis!(p, "velocity")
@gif for N in 1:n_steps
x, v = integrator(A, Δt, x, v)
push!(p, x, v) # Add new positions to the plot ...
end
end
plot_phase_space(x -> -Vho'(x), 0.1, 200)
Exercise:
A much more stable integrator than the euler
we used so far is the verlocity Verlet:
# solution
# Some example plots and parameters:
# plot_dynamics_line(Vho, 0.1, 200, integrator=verlet, ylim=(0, 2.5), xlim=(-1, 3))
# plot_phase_space(x -> -Vho'(x), 0.1, 200, integrator=verlet, xlim=(-2, 2), ylim=(-2, 2))
# plot_dynamics_line(Vmorse, 0.03, 200, integrator=verlet, xlim=(-1, 4), ylim=(-10, 5), x0=2.0)
Now we want to simulate multiple identical partices in 2D. For a system of $N$ particles in 2D, we define the particle positions as the matrix $$ \textbf{x} = \left( \begin{array}{cccc} x_{1x} & x_{2x} & \cdots & x_{Nx} \\ x_{1y} & x_{2y} & \cdots & x_{Ny} \end{array} \right) = \left( \vec{x}_1 \vec{x}_2 \cdots \vec{x}_N \right). $$ with the individual particle vectors as columns. We assume our particles interact via the idential pair potential $V_\text{pair}(r)$ depending only on particle distance $r$. The total potential is therefore $$ V_\text{tot}(\textbf{x}) = \sum_{i=1}^N \sum_{j=i+1}^N V_\text{pair}(\| \vec{x}_i - \vec{x}_j \|). $$
Program the total potential function for a matrix $\textbf{x}$. A useful function is norm
from the LinearAlgebra
package.
using LinearAlgebra
# You're code here
Now we want to plot the dynamics in a plane. In the following function the acceleration is computed using automatically generated derivatives.
function plot_dynamics_plane(Vpair, Δt, n_steps; x0=randn(2, 2), v0=zero(x0), integrator=verlet,
xlim=(-3, 3), ylim=(-3, 3))
# Total acceleration via automatic differentiation
V(x) = Vtot(Vpair, x)
Atot(x) = -V'(x)
t, x, v = 0, x0, v0 # Initialise state
@gif for i in 1:n_steps
# Propagate dynamics
x, v = integrator(Atot, Δt, x, v)
t += Δt
timestr = @sprintf "%.2f" t # Format time
# Plot the particles and their velocities
p = scatter(x[1, :], x[2, :], xlim=xlim, ylim=ylim, label="Particles at t=$timestr")
label = "Velocity / 4"
for (xi, vi) in zip(eachcol(x), eachcol(v))
plot!(p, [Tuple(xi), Tuple(xi + 0.25vi)], arrow=1.0, colour=:red, label=label)
label = ""
end
end
end
Δt = 0.07
n_steps = 300
x0 = [[0.; 0.] [1.; 0.] [-1.; 0.0]]
plot_dynamics_plane(Vmorse, Δt, n_steps; x0=x0)
Δt = 0.07
n_steps = 300
x0 = [[0.; 0.] [1.; 0.] [-1.; 0.15]]
plot_dynamics_plane(Vmorse, Δt, n_steps; x0=x0)
Δt = 0.05
n_steps = 300
x0 = [[0.; 1.] [1.; 0.] [-1.; 0] [0; -1.2]]
plot_dynamics_plane(Vmorse, Δt, n_steps; x0=x0, ylim=(-10, 10))
Δt = 0.05
n_steps = 300
x0 = 4randn(2, 10)
plot_dynamics_plane(Vmorse, Δt, n_steps; x0=x0, xlim=(-10, 10), ylim=(-10, 10))