PS4_part1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fast PDE/IE course, Skoltech, Spring 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Problem 1 (70 pts). The $N$-Body Problem\n",
    "\n",
    "The motion of more than two orbiting bodies (the $N$-body problem) may be described by the Newton second law\n",
    "$$\n",
    "\\frac{d^2}{dt^2} {\\bf x}_j = {\\bf g} ({\\bf x}_j),\n",
    "$$\n",
    "where $M_i$ are masses, \n",
    "${\\bf x}_i$ are coordinates of particles and ${\\bf g}$ is the force field,\n",
    "which may be written in terms of the gravitational potential ${\\bf g} \\equiv - \\nabla \\phi$.\n",
    "\n",
    "The $N$-body problem is a challenging problem that had eluded mathematicians\n",
    "since Newton’s time.\n",
    "In honour of his 60th birthday on January 21, 1889, Oscar II, King of Sweden and\n",
    "Norway, established a prize for anyone who could ﬁnd the solution to the problem. The prize was awarded to\n",
    "Poincaré, even though he did not solve the original problem. Only a century later in 1990 was an analytical\n",
    "solution to the $n$-body problem ﬁnally found, in the form of a converging power series.\n",
    "\n",
    "Fortunately, it is signiﬁcantly easier to approximate the solution to this famous problem using\n",
    "numerical techniques.\n",
    "However, from the computational point of view this deceptively simple ODE is greatly complicated by two considerations. Firstly, ${\\bf g}$ is aﬀected by the\n",
    "presence of all particles within the system, and computing the forces upon $N$ particles due their interaction\n",
    "with $N$ particles by direct summation is an $\\mathcal{O}(N^2)$ operation. Secondly, the particles are mobile, and the\n",
    "forces upon the $N$ particles must be repeatedly computed at each new time-step. Considering that galaxies\n",
    "are comprised of billions of stars, the complexity of computations quickly become prohibitive.\n",
    "\n",
    "Fortunately, once accuracy is fixed, fast algorithms can be used. These algorithms are based on the idea of using approximate low-parametric representations of potentials. They yield $\\mathcal{O}(N \\log N)$ (e.g., Barnes–Hut algorithm) or even with $\\mathcal{O}(N)$ complexity (fast multipole method). With the help of these algorithms recent simulations of galactic motion, such as the millenium simulation\n",
    "project by the Max Planck Institute for Astrophysics, have included over $300$ billion particles simulated over\n",
    "the equivalent span of more than $13$ billion years! \n",
    "\n",
    "* (5 pts) Estimate how much time requires to do one matrix-by-vector multiplication for $300$ billion particles with $\\mathcal{O}(N^2)$ and $\\mathcal{O}(N)$ complexity on a supercomputer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The Low-Rank nature of the FMM\n",
    "\n",
    "Recently the low-rank interpretation of the Fast Multipole Method was developed. The $\\mathcal{O}(N \\log N)$ complexity can be obtained with $\\mathcal{H}$-matrices and $\\mathcal{O}(N)$ with $\\mathcal{H^2}$-matrices. In this part of the problem set you will investigate the low-rank nature of the FMM and will implement matrix-by-vector multiplication with the help of $\\mathcal{H}$-matrices.\n",
    "\n",
    "In this example we will investigate low-rank nature of the interaction between well-separated clusters of particles.\n",
    "Let ${\\bf y}_j$, $j=1\\dots N$ be source point charges located randomly within the unit square $[0,1]^2$ and let ${\\bf x}_j$, $j=1\\dots N$ be receivers points randomly located within a speciﬁed target square.\n",
    "\n",
    "#### Newton potential\n",
    "Consider the Newtonial gravitational potential function $\\phi$ as follows:\n",
    "$$\n",
    "\\phi({\\bf x}) = \\frac{1}{|{\\bf x}|}\n",
    "$$\n",
    "The dense coeﬃcient matrix $G$ associated with the potential interaction between ${\\bf x}_m$ and ${\\bf y}_n$ under Newtonian gravity contains an element on the $j$-th row and $k$-th column equal to:\n",
    "$$\n",
    "G_{mn} = \\frac{1}{|{\\bf x}_m - {\\bf y}_n|}\n",
    "$$\n",
    "Compute the relative rank of this matrix, by performing a singular value decomposition \n",
    "$G\\approx U\\Sigma V^T$,\n",
    "so that the relative spectral error of the reconstruction \n",
    "$\\|G- U\\Sigma V^T\\|_2/\\|G\\|_2 \\leqslant \\epsilon$,\n",
    "where $\\epsilon = 10^{-5}$.\n",
    "The relative rank of the matrix is equal to the number of non-zero singular values with spectral norms\n",
    "$\\geqslant \\sigma_1 \\epsilon$.\n",
    "\n",
    "* (5 pts) Use 20 source points in the unit square with lower left corner at the origin, and 20 test points in a target unit square. Produce a plot of the rank of the associated $G$ matrix, as a function of the target square position. In particular, place the target square’s bottom left corner at $[2, p]$ and the top right corner at $[3, p+ 1]$. Vary $p$ from $0$ to $5$ in steps of $0.5$.\n",
    "\n",
    "\n",
    "* (5pts) With the bottom left corner of the target square placed at [2,4] and the top right corner placed at [3,5]. Vary the number of interacting points N from 20 to 200 in steps of 20. Does increasing the number of interactions and the size of the matrix G signiﬁcantly aﬀect its rank?\n",
    "\n",
    "\n",
    "#### High-Frequency issues\n",
    "\n",
    "Suppose that Newtonian gravity were instead governed by the Helmholtz Green’s function as follows:\n",
    "$$\n",
    "\\phi({\\bf x}) = \\frac{e^{ik|{\\bf x}|}}{|{\\bf x}|}.\n",
    "$$\n",
    "The dense coeﬃcient matrix $G$ associated with the interaction between ${\\bf x}_m$ and ${\\bf y}_n$ associated with\n",
    "this new potential function contains an element on the $m$-th row and $n$-th column equal to:\n",
    "$$\n",
    "G_{mn} = \\frac{e^{ik|{\\bf x}_m - {\\bf y}_n|}}{|{\\bf x}_m - {\\bf y}_n|}.\n",
    "$$\n",
    "\n",
    "* (5 pts) Use $20$ source points in the unit square with lower left corner at the origin, and $20$ test points in a target unit square. Let $k = 10$, produce a plot of the rank of the associated $G$ matrix, as a function of the target square position. In particular, place the target square’s bottom left corner at $[2, p]$ and the top right corner at $[3, p+ 1]$. Vary $p$ from $0$ to $5$ in steps of $0.5$.\n",
    "\n",
    "\n",
    "* (5 pts)  Use $20$ source points and $20$ test points, and set the bottom left corner of the target square placed to $[2,4]$ and the top right corner to $[3,5]$. Produce a plot of the rank of the interaction as a function of $k$, varying from $1$ to $100$. What are the diﬃculties presented by having a high value of $k$?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### The $\\mathcal{H}$-matrices\n",
    "Let ${\\bf y_j} $, $j=1\\dots N$ be coordinates of randomly distributed particles in the 2D box $[0,1]^2$. Let reciever points ${\\bf x}_j $ be the same points ${\\bf y_j}$. Consider Newton potential from the previous part and implement matvec operation with the matrix $G$ in the $\\mathcal{H}$-format with $\\mathcal{O}(N \\log N)$.\n",
    "\n",
    "The idea of $\\mathcal{H}$-matrices is quite simple. First, we separate our particles in a quad-tree structure of particle groups. If the group on some tree level is sufficiently far away from another group, we can approximate its gravitational effects by the low-rank approximation with some accuracy $\\epsilon$. If it is not, we compute its influence recursively by using lower tree levels. \n",
    "Note that the group of sources and receivers are separated if the paramenter $\\theta = s / d$ is small enough (here $s$ is the width of the region represented by the internal node, and $d$ is the distance between the body and the node's center-of-mass).\n",
    "\n",
    "\n",
    "* (8 pts) Propose an algorithm for the quadtree construction (here tree for sources = tree for receivers). How much memory the tree storage takes? Can you reach $\\mathcal{O}(N)$ memory for the storage? Propose a way to store the tree in Python and write the program that computes the tree, given the location of the particles. What do you need to store in each node of the tree?\n",
    "\n",
    "\n",
    "* (17 pts) Write $\\mathcal{H}$-matrix by vector multiplication. The program should consist of three parts:\n",
    " 1.  Tree and interaction list construction, given the location of the particles, geometric constant $\\theta$ and $\\epsilon$-truncation level\n",
    " 2.  Filling the information in the tree \n",
    " 3.  Computing the matvec\n",
    "\n",
    "\n",
    "* (10 pts) Compare the results computed by direct evaluation and by $\\mathcal{H}$-matrices. Study the dependence of accuracy and computational cost on the geometric parameter $\\theta$ and $\\epsilon$\n",
    "    \n",
    "    \n",
    "* (5 pts) Explain in a nutshell the $\\mathcal{H}^2$-matrices as compared with the $\\mathcal{H}$-matrices\n",
    "\n",
    "\n",
    "* (5 pts) Explain how does $\\mathcal{H}^2$-matrices relate to the FMM. Give exapmles where low-rank interpretation is preferable"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 160%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}