# Simulating n-bodies and functional programming

### The n-bodies problem

There’s been some activity in the #haskell IRC channel recently towards coming up with a better solution to the n-bodies problem.  In this problem, you’re given a set of bodies, and are asked to simulate Newtonian gravity on them.  Most of the solutions I’ve seen go past have involved some sort of crazy state monad based stuff, reading from and writing to STUArrays and other imperative ideas.  So, I thought I’d try and solve the problem using a pure functional style instead.

### Reactive

I’ve recently been playing with Reactive a lot, initially for my work, but also just because it’s so damned cool.  Reactive is a library for Haskell that lets you describe time varying values in a purely functional way.  That is, out goes the concept of state, and in comes the concept of describing exactly what is going on.

### Back to n-bodies

In order to simulate newtonian physics, we need the bodies’ mass,the positions of the bodies, and the velocities of the bodies:

```data Body = Body { mass     :: Double
, position :: Point3  Double
, velocity :: Vector3 Double
}```

That was pretty straight forward, lets get down to the core of the problem – simulating gravity:

```accelBetween :: Body -> Body -> Vector3 Double
accelBetween p1 p2 =
(mass p2 *^ vecBetween) / ((mag vecBetween) ^ 3)
where
vecBetween = position p2 ^-^ position p1

computeAccels :: Body -> [Body] -> Vector3 Double
computeAccels p = sum . map (accelBetween p)```

Well, that was surprisingly easy!  The `accelBetween` function describes the newtonian gravity equation just as the maths does, and summing up the acceleration due to all the other planets in the system is fairly straightforward.

What we’ve not seen yet though, is any of that lovely functional reactive programming I was talking about.  We know now how to compute the acceleration affecting any one body at any one time, but what we want to know is the acceleration affecting a body at *any* time.  To do this, we’re going to need to move from doing computations on `Body`s, to instead doing computation on `Behavior`s of `Body`s.

```bodyAccels :: [Behaviour Body] -> [Behaviour (Vector3 Double)]
bodyAccels ps = withOthers (liftA2 computeAccels) ps

withOthers :: (a -> [a] -> b) -> [a] -> [b]
withOthers f xs =
withOthers' f [] xs
where
withOthers' _ _  []       = []
withOthers' f os (x : xs) = f x (os ++ xs) : withOthers' f (x : os) xs
```

The `withOthers` function here is just like map, but it passes in all the other values in the list in a second argument to the function.

So then, `bodyAccels` computes a continuous acceleration function for all the bodies in the system.  For each body, it runs the `computeAccels` function, giving it all other bodies in the system as its second argument.  Crucially, `liftA2` allows us to do this in the `Behavior` Applicative, so we are no longer computing it on rigid, static bodies, but instead, on all the positions and velocities the bodies in the system will ever have (isn’t lazyness great!).

Finally, we can get from these accelerations down to the velocities, and then positions of the bodies using integration on the acceleration:

```bodyVel :: Body -> Behaviour (Vector3 Double) -> Behaviour (Vector3 Double)
bodyVel p acc = (velocity p ^+^) <\$> integral dt acc```

I’m sure you can imagine what the bodyPos function looks like.  You may wonder what the `dt` in here is talking about.  This is an unfortunate effect of not being able to mathematically integrate arbitrary functions.  Instead, we must use euler integration, and that requires us to provide times at which to take samples.  The `dt` argument is an event which ticks reasonably fast, and progresses our simulation:

```dt :: Event ()
dt = atTimes [0,0.01..]```

So, now we are able to combine all our efforts together, and solve the whole n-bodies problem:

```nbodies :: [Body] -> [Behaviour Body]
nbodies ps = pbs
where
pbs = Body <\$> (mass <\$> ps) <*> pps <*> pvs
pps = zipWith bodyPos ps pvs
pvs = zipWith bodyVel ps pas
pas = bodyAccels pbs```

### Conclusions

A lot of assumptions are made about how we must write programs.  Often, even beautiful mathematical problems end up described as horrible state-full masses of code that obscure what it is we’re trying to compute.  I’ve presented a solution to the  n-bodies problem using the Reactive library to get a handle on time in a purely functional setting, it turned out that this was rather beautiful!

Full code for my solution can be found below.

```module NBodies (Body(..), nbodies) where

import Control.Applicative

import Data.VectorSpace
import Data.AffineSpace

import FRP.Reactive

import Graphics.FieldTrip

data Body = Body { mass     :: Double
, position :: Point3  Double
, velocity :: Vector3 Double
}

nbodies :: [Body] -> [Behaviour Body]
nbodies ps = pbs
where
pbs = Body <\$> (mass <\$> ps) <*> pps <*> pvs
pps = zipWith bodyPos ps pvs
pvs = zipWith bodyVel ps pas
pas = bodyAccels pbs

accelBetween :: Body -> Body -> Vector3 Double
accelBetween p1 p2 =
(mass p2 *^ vecBetween) / ((mag vecBetween) ^ 3)
where
vecBetween = position p2 ^-^ position p1

computeAccels :: Body -> [Body] -> Vector3 Double
computeAccels p = sum . map (accelBetween p)

bodyAccels :: [Behaviour Body] -> [Behaviour (Vector3 Double)]
bodyAccels ps = withOthers (liftA2 computeAccels) ps

withOthers :: (a -> [a] -> b) -> [a] -> [b]
withOthers f xs =
withOthers' f [] xs
where
withOthers' _ _  []       = []
withOthers' f os (x : xs) = f x (os ++ xs) : withOthers' f (x : os) xs

bodyVel :: Body -> Behaviour (Vector3 Double) -> Behaviour (Vector3 Double)
bodyVel p acc = (velocity p ^+^) <\$> integral dt acc

bodyPos :: Body -> Behaviour (Vector3 Double) -> Behaviour (Point3  Double)
bodyPos p vel = (position p .+^) <\$> integral dt vel

dt :: Event ()
dt = atTimes [0,0.01..]
```

## 10 comments on “Simulating n-bodies and functional programming”

1. axman6 says:

I think it should be mentioned that the original problem being discussed was the shootout  nbodies problem, who’s soul purpose is to iterate n bidoes 50,000,000 times in the shortest possible time. But good work man, I’ll have to play with this sometime.

• beelsebob says:

The reason I didn’t mention that is that the n-bodies problem is a pretty general problem, it’s not constrained solely to the great language shootout. I intend to take this further, and use it within an app for simulating planetary physics.

2. Shane says:

It’s very neat to see the n-bodies problem done in reactive form. Thank you for the article. I’ve ported a simple physics engine to Haskell, and one thing I was wondering about is how one might optimize for anti-symmetric forces like gravity? For example, Suppose the force on body 2 wrt 1 is the negative of the force on body 1 wrt 2, F_{12} = -F_{21}. Is there a nice way to functionally avoid having two function evaluations?

• beelsebob says:

That’s a really good question, and was one of the first things I was thinking of doing with this when I started optimising it – I’ll get back to you.

3. Chris Dew says:

I think I may have just started to (needlessly) re-invent reactive programming, as I hadn’t heard of it before I read this post. Thanks for the enlightenment.

Is what I’m attempting (in the link above) reactive programming, or have I misunderstood?

• beelsebob says:

Yes, that looks very much like FRP in the making 🙂

4. wagstaff says:

This code appears to be a long way from working. So fas as I can see, almost no part of it typechecks.

• beelsebob says:

Yep, the code was written when I didn’t even have a working Reactive version, so it’s not bang on right. It only takes a few tweaks to get it there though. Having said that, if you do tweak it, unfortunately, Reactive is still a little too buggy to run it. Close, but no cigar.

5. wagstaff says:

“A few tweaks” – perhaps depends on your level of expertise!

I’ve got most of this compiling, but I don’t suppose you’d care to provide hints for corrected versions of bodyAccels and nbodies, would you?

6. agrigioli says:

What Word Press theme do you use?