In [1]:
from IPython.html.services.config import ConfigManager
from IPython.utils.path import locate_profile
cm = ConfigManager(profile_dir=locate_profile(get_ipython().profile))
cm.update('livereveal', {
'theme': 'sky',
'transition': 'zoom',
'start_slideshow_at': 'selected',
})
Out[1]:
The discretization of the integral equation leads to dense matrices.
The main question is how to compute the matrix-by-vector product,
i.e. the summation of the form:
$$\sum_{j=1}^M A_{ij} q_j = V_i, \quad i = 1, \ldots, N.$$The matrix $A$ is dense, i.e. its element can not be omitted. The complexity is $\mathcal{O}(N^2)$.
Can we make it faster?
The simplest case is the computation of the potentials from the system of charges
$$V_i = \sum_{j} \frac{q_j}{\Vert r_i - r_j \Vert}$$This summation appears in:
It is called the N-body problem .
There is no problem with memory, since you only have two cycles.
In [34]:
import numpy as np
import math
from numba import jit
N = 10000
x = np.random.randn(N, 2);
y = np.random.randn(N, 2);
charges = np.ones(N)
res = np.zeros(N)
@jit
def compute_nbody_direct(N, x, y, charges, res):
for i in xrange(N):
res[i] = 0.0
for j in xrange(N):
dist = (x[i, 0] - y[i, 0]) ** 2 + (x[i, 1] - y[i, 1]) ** 2
dist = math.sqrt(dist)
res[i] += charges[j] / dist
In [35]:
%timeit compute_nbody_direct(N, x, y, charges, res)
One of the most famous N-body computations is the Millenium run
In [23]:
from IPython.display import YouTubeVideo
YouTubeVideo('UC5pDPY5Nz4')
Out[23]:
The particle systems can be used to model a lot of things.
For nice examples, see the website of Ron Fedkiw
In [22]:
from IPython.display import YouTubeVideo
YouTubeVideo('6bdIHFTfTdU')
Out[22]:
The idea of Barnes and Hut was to split the sources into big blocks using the cluster tree
The algorithm is recursive. Let $\mathcal{T}$ be the tree, and $x$ is the point where we need to compute the potential.
The complexity is $\mathcal{O}(\log N)$ for 1 point!
Another popular choice of the tree is the KD-tree
The construction is simple as well:
Split along x-axis, then y-axis in such a way that the tree is balanced (i.e. the number of points in the left child/right child is similar).
The tree is always balanced, but biased towards the coordinate axis.
For each node of the tree, we need to compute its total mass and the center-of-mass. If we do it in a cycle, then the complexity will be $\mathcal{O}(N^2)$ for the tree constuction.
However, it is easy to construct it in a smarter way.
Now you can actually code this (minor things remaining are the bounding box and separation criteria).
Well, several
The original BH method has low accuracy, and is based on the expansion
$$f(x, y) \approx f(x_0, y_0)$$What to do?
Answer: Use higher-order expansions!
$$f(x + \delta x, y + \delta y) \approx f(x, y) + \sum_{k, l=0}^p (D^{k} D^{l} f) \delta ^k \delta y^l \frac{1}{k!} \frac{1}{l!} + \ldots $$For the Coloumb interaction $\frac{1}{r}$ we have the multipole expansion
$$ v(R) = \frac{Q}{R} + \frac{1}{R^3} \sum_{\alpha} P_{\alpha} R_{\alpha} + \frac{1}{6R^5} \sum_{\alpha, \beta} Q_{\alpha \beta} (3R_{\alpha \beta} - \delta_{\alpha \beta}R^2) + \ldots, $$where $P_{\alpha}$ is the dipole moment, $Q_{\alpha \beta}$ is the quadrupole moment (by actually, nothing more than the Taylor series expansion).
For cosmology this problem is so important, so that they have released a special hardware Gravity Pipe for solving the N-body problem
Sidenote, When you Google for "FMM", you will also encounter fast marching method (even in the scikit).
Everyone uses its own in-house software, so a good Python open-source software is yet to be written.
This is also a perfect test for the GPU programming (you can try to take such project in the App Period, by the way).
In [2]:
from IPython.core.display import HTML
def css_styling():
styles = open("./styles/alex.css", "r").read()
return HTML(styles)
css_styling()
Out[2]: