Solutions_to_Exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [01_Variables_Control_Packages](01_Variables_Control_Packages.ipynb)\n",
    "\n",
    "Compute $2^{100}$:\n",
    "- As fast as possible\n",
    "- As exact as possible"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This will give the wrong result:\n",
    "2^100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# An exact way to compute the answer is using arbitrary-precision integers\n",
    "BigInt(2)^100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Usually faster than arbitrary-precision computations are\n",
    "# floating point operations (might not be true for this simple case, however)\n",
    "2.0^100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "Compute\n",
    "$$ 15! \\qquad 100! \\qquad \\left(\\begin{array}{c} 100 \\\\ 15 \\end{array}\\right) $$\n",
    "with the Julia you know so far."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "res15 = 1\n",
    "for i in 1:15\n",
    "    res15 = res15 * i\n",
    "end\n",
    "\n",
    "res100 = BigInt(1)\n",
    "for i in 1:100\n",
    "    res100 = res100 * i\n",
    "end\n",
    "\n",
    "res85 = BigInt(1)\n",
    "for i in 1:(100 - 15)\n",
    "    res85 = res85 * i\n",
    "end\n",
    "\n",
    "res100 / res85 / res15"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [02_Functions_Types_Dispatch](02_Functions_Types_Dispatch.ipynb)\n",
    "\n",
    "### Exercise\n",
    "- 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$.\n",
    "- Now do the same thing using the confining potential:\n",
    "  $$ V_\\alpha(x) = \\left\\{ \\begin{array}{ll} (|x| - 1)^\\alpha & |x| > 1 \\\\ 0 & \\text{else}\\end{array}\\right.$$\n",
    "  for $\\alpha = 4$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "function acceleration(x)\n",
    "    if x > 1 || x < -1\n",
    "        -4 * sign(x) * (abs(x) - 1)^3\n",
    "    else\n",
    "        0\n",
    "    end\n",
    "end\n",
    "\n",
    "x = 0.\n",
    "v = 1.\n",
    "for i in 1:5\n",
    "    x, v = euler(acceleration, 0.1, x, v)\n",
    "end\n",
    "println(x, \" \", v)\n",
    "\n",
    "#\n",
    "# A shorter, but equivalent way using the \"do\" Syntax\n",
    "#\n",
    "\n",
    "x = 0.\n",
    "v = 1.\n",
    "for i in 1:5\n",
    "    x, v = euler(0.1, x, v) do x\n",
    "        if x > 1 || x < -1\n",
    "            -4 * sign(x) * (abs(x) - 1)^3\n",
    "        else\n",
    "            0\n",
    "        end\n",
    "    end\n",
    "end\n",
    "println(x, \" \", v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "Which of the following type are subtypes of another?\n",
    "Try to guess first and then verify by using the operator `<:`.\n",
    "\n",
    "```julia\n",
    "Float64     AbstractFloat      Integer\n",
    "Number      AbstractArray      Complex\n",
    "Real        Any                Nothing\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following type chains (subtype ``<:`` supertype) are true:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Float64 <: AbstractFloat <: Real <: Number <: Any"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Integer <: Real <: Number <: Any"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Complex <: Number <: Any"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Nothing <: Any"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "AbstractArray <: Any"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [03_Arrays](03_Arrays.ipynb)\n",
    "\n",
    "### Exercise\n",
    "Create the following arrays using Julia code:\n",
    "$$\\left(\\begin{array}{ccccc}\n",
    "   2&2&2&2&2 \\\\\n",
    "   2&2&2&2&2 \\\\\n",
    "   2&2&2&2&2 \\\\\n",
    "   \\end{array}\\right) \\qquad\n",
    "   \\left(\\begin{array}{cccc}\n",
    "   0.1&0.5&0.9&1.3\\\\\n",
    "   0.2&0.6&1.0&1.4\\\\\n",
    "   0.3&0.7&1.1&1.5\\\\\n",
    "   0.4&0.8&1.2&1.6\\\\\n",
    "   \\end{array}\\right)\n",
    "   $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ones(Int, 3, 5) + ones(Int, 3, 5)\n",
    "\n",
    "# or \n",
    "\n",
    "2ones(Int, 3, 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "reshape(collect(1:16), 4, 4) / 10"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Propagating a particle in 2D\n",
    "We want improve our capabilities to perform simulations to two dimensions using the 2D harmonic potential\n",
    "$$ V\\left( \\begin{array}{c} x_1 \\\\ x_2 \\end{array}\\right) = x_1^2 + x_2^2 $$\n",
    "and the acceleration map\n",
    "$$ \\vec{A}_V = - \\nabla V. $$\n",
    "The adapted forward-Euler scheme is\n",
    "$$\\left\\{\\begin{array}{l}\n",
    "\\vec{v}^{(n+1)} = \\vec{v}^{(n)} + \\vec{A}_V(\\vec{x}^{(n)}) \\Delta t\\\\\n",
    "\\vec{x}^{(n+1)} = \\vec{x}^{(n)} + \\vec{v}^{(n)} \\Delta t\\\\\n",
    "\\end{array}\\right. .$$\n",
    "\n",
    "Change the `euler` function we introduced in [02_Functions_Types_Dispatch.ipynb](02_Functions_Types_Dispatch.ipynb)\n",
    "accordingly (Hint: Suprisingly few changes are needed) and run the dynamics for five steps using\n",
    "$$x_n = \\left( \\begin{array}{c} 0 \\\\ 0 \\end{array}\\right) \\qquad v_n = \\left( \\begin{array}{c} 1 \\\\ 0 \\end{array}\\right) .$$\n",
    "Check against your previous implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The previous euler implementation can be used without any additional change:\n",
    "euler(A, Δt, xn, vn) = xn + vn * Δt, vn + A(xn) * Δt\n",
    "\n",
    "# Assuming x to be the vector (x_1 \\\\ x_2), the derivative of V can be written as:\n",
    "A2d(x) = -2x\n",
    "\n",
    "# Therefore propagating 5 steps is equal to:\n",
    "x = [0, 0]\n",
    "v = [1, 0]\n",
    "for i in 1:5\n",
    "    x, v = euler(A2d, 0.1, x, v)\n",
    "end\n",
    "println(x)\n",
    "println(v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# [05_Dancing_Particles](05_Dancing_Particles.ipynb)\n",
    "\n",
    "**Exercise:**\n",
    "\n",
    "A much more stable integrator than the `euler` we used so far is the verlocity Verlet:\n",
    "\n",
    "$$\\left\\{\\begin{array}{l}\n",
    "x^{(n+1)} = x^{(n)} + v^{(n)} \\Delta t + \\frac{A_V(x^{(n)})}{2} \\Delta t^2\\\\\n",
    "v^{(n+1)} = v^{(n)} + \\frac{A_V(x^{(n))} + A_V(x^{(n+1)})}{2} \\Delta t\\\\\n",
    "\\end{array}\\right. $$\n",
    "\n",
    "- 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 \n",
    "  of the dynamics. For our purposes there is no need to keep an eye on that.\n",
    "- How does the previous dynamics look like in this example. Does this algorithm conserve energy (phase-space plot)?\n",
    "- Also look at the Morse potential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# One way to code the verlet function is:\n",
    "function verlet(A, Δt, xn, vn)\n",
    "    An = A(xn)\n",
    "    x_next = xn + vn * Δt + An/2 * Δt^2\n",
    "    v_next = vn + (An + A(x_next)) / 2 * Δt\n",
    "    x_next, v_next\n",
    "end"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "Program the total potential function for a matrix $\\textbf{x}$. A useful function is `norm` from the `LinearAlgebra` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# One solution:\n",
    "function Vtot(Vpair, x)\n",
    "    n_particles = size(x, 2)\n",
    "    accu = zero(eltype(x))   # Get a zero of the appropriate type\n",
    "    for i in 1:n_particles, j in i+1:n_particles\n",
    "        accu += Vpair(norm(x[:, i] .- x[:, j]))\n",
    "    end\n",
    "    accu\n",
    "end"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.3.0",
   "language": "julia",
   "name": "julia-1.3"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}