# 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