N-Body Simulations


An N-body simulation approximates the motion of particles that interact with one another through some type of physical force, such as gravitational or electric forces.

The most obvious way to calculate these forces is to directly integrate each force acting on each body. This is known as the "brute force" method and has the major drawback having N² time complexity where N is the number of particles.

For the most accurate results, the brute force method is the preferred method. However, as the number of particles doubles, the computational time of the algorithm quadruples. Unfortunately, this is very slow when working with a system with many particles.

The most common way to improve this computational time is to use an algorithm known as the Barnes-Hut algorithm. It does so by sorting each particle into leaves of a quadtree data structure, and then decreasing the number of force calculations per particle by "neglecting" other particles that are sufficiently far away from it.

The groups of "neglected" particles that are far enough away from the individual particle are treated as one mass and the force calculation is done based on that group's center of mass. This turns out to be a reasonable approximation to make, as gravitational and electric forces drop off in intensity with 1/r².

For the Barnes-Hut algorithm, each body calculates a force with on average log(N) other bodies in the tree. In effect, the overall time for one time-step in the Barnes-Hut algorithm is proportional to Nlog(N), a substantial improvement over the brute force method.

Further Reading:


        <!DOCTYPE html>
        <meta charset="utf-8">
        <body>
            <div style="text-align: center" class="svgContainer"></div>
        <script src="d3.v3.min.js"></script>
        <script>
        
            var width = 960,
                height = 500,
                svg,
                nodes,
                node,
                force,
                forceStrength = 0.4,
                np = 100;
        
            function render(n, str){
        
                d3.select("svg").remove();
        
                var chargeStrength = parseFloat(str);
                var numOfParticles = parseInt(n);
        
                nodes = d3.range(numOfParticles).map(function(i){
                    return {
                        index: i,
                        radius: Math.random() * 4 + 3,
                        color: "steelblue"
                    }
                });
        
                svg = d3.select(".svgContainer").append("svg")
                    .attr("width", width)
                    .attr("height", height)
                    .style("border", "1px solid black");
        
                node = svg.selectAll(".node")
                    .data(nodes)
                        .enter().append("circle")
                        .attr("class", "node")
                        .attr("cx", function(d) { return d.x })
                        .attr("cy", function(d) { return d.y })
                        .attr("r", function(d) { return d.radius })
                        .style("fill", function(d) { return d.color; })
                        .style("stroke", "black");
            
                force = d3.layout.force()
                    .nodes(nodes)
                    .size([width, height])
                    .gravity(0)
                    .charge(chargeStrength)
                    .friction(1)
                    .on("tick", tick)
                    .start();
            }
            
        
            function tick(e){
        
                var q = d3.geom.quadtree(nodes);
        
                var k = 0.2;
                nodes.forEach(function(d,i){
                    q.visit(collide(nodes[i]));
                });
        
                node.attr("cx", function(d) { return d.x; })
                    .attr("cy", function(d) { return d.y; })
                    .style("fill", function(d) { return d.color; });
                
                force.resume();
            }
        
            function collide(node){
                var r = node.radius,
                    nx1 = node.x - r,
                    nx2 = node.x + r,
                    ny1 = node.y - r,
                    ny2 = node.y + r;
                return function(quad, x1, y1, x2, y2) {
                    if (quad.point && (quad.point !== node)) {
                        var x = node.x - quad.point.x,
                            y = node.y - quad.point.y,
                            l = Math.sqrt(x * x + y * y),
                            r = node.radius + quad.point.radius;
                        if (l < r) {
                            l = (l - r) / (l * 2);
                            node.x -= x *= l;
                            node.y -= y *= l;
                            quad.point.x += x;
                            quad.point.y += y;
                        }
                    }
                    return x1 > nx2 || x2 < nx1 || y1 > ny2 || y2 < ny1;
                }
            }
        
            render(np, forceStrength);
        
        </script>