01_Variables_Control_Packages

Compute $2^{100}$:

  • As fast as possible
  • As exact as possible
In [ ]:
# This will give the wrong result:
2^100
In [ ]:
# An exact way to compute the answer is using arbitrary-precision integers
BigInt(2)^100
In [ ]:
# Usually faster than arbitrary-precision computations are
# floating point operations (might not be true for this simple case, however)
2.0^100

Exercise

Compute $$ 15! \qquad 100! \qquad \left(\begin{array}{c} 100 \\ 15 \end{array}\right) $$ with the Julia you know so far.

In [ ]:
res15 = 1
for i in 1:15
    res15 = res15 * i
end

res100 = BigInt(1)
for i in 1:100
    res100 = res100 * i
end

res85 = BigInt(1)
for i in 1:(100 - 15)
    res85 = res85 * i
end

res100 / res85 / res15

02_Functions_Types_Dispatch

Exercise

  • Take timesteps of $\Delta t = 0.1$ and start at $x_n = 0$ and $v_n = 1$. Propagate the dynamics for 5 steps in the harmonic potential $V = x^2$.
  • Now do the same thing using the confining potential: $$ V_\alpha(x) = \left\{ \begin{array}{ll} (|x| - 1)^\alpha & |x| > 1 \\ 0 & \text{else}\end{array}\right.$$ for $\alpha = 4$.
In [ ]:
function acceleration(x)
    if x > 1 || x < -1
        -4 * sign(x) * (abs(x) - 1)^3
    else
        0
    end
end

x = 0.
v = 1.
for i in 1:5
    x, v = euler(acceleration, 0.1, x, v)
end
println(x, " ", v)

#
# A shorter, but equivalent way using the "do" Syntax
#

x = 0.
v = 1.
for i in 1:5
    x, v = euler(0.1, x, v) do x
        if x > 1 || x < -1
            -4 * sign(x) * (abs(x) - 1)^3
        else
            0
        end
    end
end
println(x, " ", v)

Exercise

Which of the following type are subtypes of another? Try to guess first and then verify by using the operator <:.

Float64     AbstractFloat      Integer
Number      AbstractArray      Complex
Real        Any                Nothing

The following type chains (subtype <: supertype) are true:

In [ ]:
Float64 <: AbstractFloat <: Real <: Number <: Any
In [ ]:
Integer <: Real <: Number <: Any
In [ ]:
Complex <: Number <: Any
In [ ]:
Nothing <: Any
In [ ]:
AbstractArray <: Any

03_Arrays

Exercise

Create the following arrays using Julia code: $$\left(\begin{array}{ccccc} 2&2&2&2&2 \\ 2&2&2&2&2 \\ 2&2&2&2&2 \\ \end{array}\right) \qquad \left(\begin{array}{cccc} 0.1&0.5&0.9&1.3\\ 0.2&0.6&1.0&1.4\\ 0.3&0.7&1.1&1.5\\ 0.4&0.8&1.2&1.6\\ \end{array}\right) $$

In [ ]:
ones(Int, 3, 5) + ones(Int, 3, 5)

# or 

2ones(Int, 3, 5)
In [ ]:
reshape(collect(1:16), 4, 4) / 10

Exercise: Propagating a particle in 2D

We want improve our capabilities to perform simulations to two dimensions using the 2D harmonic potential $$ V\left( \begin{array}{c} x_1 \\ x_2 \end{array}\right) = x_1^2 + x_2^2 $$ and the acceleration map $$ \vec{A}_V = - \nabla V. $$ The adapted forward-Euler scheme is $$\left\{\begin{array}{l} \vec{v}^{(n+1)} = \vec{v}^{(n)} + \vec{A}_V(\vec{x}^{(n)}) \Delta t\\ \vec{x}^{(n+1)} = \vec{x}^{(n)} + \vec{v}^{(n)} \Delta t\\ \end{array}\right. .$$

Change the euler function we introduced in 02_Functions_Types_Dispatch.ipynb accordingly (Hint: Suprisingly few changes are needed) and run the dynamics for five steps using $$x_n = \left( \begin{array}{c} 0 \\ 0 \end{array}\right) \qquad v_n = \left( \begin{array}{c} 1 \\ 0 \end{array}\right) .$$ Check against your previous implementation.

In [ ]:
# The previous euler implementation can be used without any additional change:
euler(A, Δt, xn, vn) = xn + vn * Δt, vn + A(xn) * Δt

# Assuming x to be the vector (x_1 \\ x_2), the derivative of V can be written as:
A2d(x) = -2x

# Therefore propagating 5 steps is equal to:
x = [0, 0]
v = [1, 0]
for i in 1:5
    x, v = euler(A2d, 0.1, x, v)
end
println(x)
println(v)

05_Dancing_Particles

Exercise:

A much more stable integrator than the euler we used so far is the verlocity Verlet:

$$\left\{\begin{array}{l} x^{(n+1)} = x^{(n)} + v^{(n)} \Delta t + \frac{A_V(x^{(n)})}{2} \Delta t^2\\ v^{(n+1)} = v^{(n)} + \frac{A_V(x^{(n))} + A_V(x^{(n+1)})}{2} \Delta t\\ \end{array}\right. $$
  • Program this algorithm, taking care that it supports multi-dimensional positions and velocities as well. In practice we would like to avoid recomputing $A_V(x)$ as much as possible, since this is usually the expensive step of the dynamics. For our purposes there is no need to keep an eye on that.
  • How does the previous dynamics look like in this example. Does this algorithm conserve energy (phase-space plot)?
  • Also look at the Morse potential
In [ ]:
# One way to code the verlet function is:
function verlet(A, Δt, xn, vn)
    An = A(xn)
    x_next = xn + vn * Δt + An/2 * Δt^2
    v_next = vn + (An + A(x_next)) / 2 * Δt
    x_next, v_next
end

Exercise

Program the total potential function for a matrix $\textbf{x}$. A useful function is norm from the LinearAlgebra package.

In [ ]:
# One solution:
function Vtot(Vpair, x)
    n_particles = size(x, 2)
    accu = zero(eltype(x))   # Get a zero of the appropriate type
    for i in 1:n_particles, j in i+1:n_particles
        accu += Vpair(norm(x[:, i] .- x[:, j]))
    end
    accu
end