Haskell Particle simulation - calculating velocities of particles

Question

I am working on a particle simulation program using Haskell. For one of the functions I am trying to determine the new velocities of all the particles in the simulation based on the mass and velocities of all the surrounding particles.

The function is of this form:

accelerate :: Float -> [Particle] -> [Particle]

Particle is a data type that contains the mass, position vector and velocity vector, the 'Float' argument represents the delta time of the respective time step in the simulation

I would like some suggestions on possible functions I can use to traverse the list while calculating the velocities of each of the particles with respect to the other particles in the list.

One possible approach I can think of:

1. assume there is another function 'velocityCalculator' which has the following definition:

   velocityCalculator :: Particle -> Particle -> (Float,Float)
   

This takes two particles and returns the updated velocity vector for the first particle.

2. apply foldl; using the above function as the binary operator, a particle and the list of particles as the arguments, i.e.

   foldl velocityCalculator particle particleList
   

3. iterate through the list of particle, applying foldl to each element and building the new list containing the particles with the updated velocities

I am not sure if this is the most efficient method so any suggestions and improvements are very much appreciated.

PLEASE NOTE -> as I have said I am only looking for suggestions not an answer!

Thanks!

Answer 1

Seems like you are pretty set on using foldl. For example

3. iterate through the list of particle, applying foldl to each element and building the new list containing the particles with the updated velocities

doesn't really make sense. You apply foldl to a list to reduce it to a "summary value", according to some binary summarizing function. It doesn't really make sense to apply it to a single particle.

I am answering this question assuming you are having trouble writing the program in the first place -- it's usually best to do this before worrying about efficiency. Let me know if I have assumed wrong.

I am not sure what rule you want to use to update the velocities, but I assume it's some sort of pairwise force simulation, for example gravity or electromagnetism. If this is so, here are some hints that will guide you to a solution.

type Vector = (Float, Float)

-- Finds the force exerted on one particle by the other.
-- Your code will be simplified if this returns (0,0) when the two
-- particles are the same. 
findForce :: Particle -> Particle -> Vector

-- Find the sum of all forces exerted on the particle
-- by every particle in the list. 
totalForce :: [Particle] -> Particle -> Vector

-- Takes a force and a particle to update, returns the same particle with
-- updated velocity. 
updateVelocity :: Vector -> Particle -> Particle

-- Calculate mutual forces and update all particles' velocities
updateParticles :: [Particle] -> [Particle]

Each of these functions will be quite short, one or two lines. If you need further hints for which higher-order functions to use, pay attention to the type of the function you are trying to write, and note

map :: (a -> b) -> [a] -> [b]           -- takes a list, returns a list
filter :: (a -> Bool) -> [a] -> [a]     -- takes a list, returns a list
foldl :: (a -> b -> a) -> a -> [b] -> a -- takes a list, returns something else
foldr :: (a -> b -> b) -> b -> [a] -> b -- takes a list, returns something else

Answer 2

You might achieve a speed up by a factor of 2 by memoization, if you give every Particle a particle_id :: Int ID and then define that:

forceOf a b | particle_id a > particle_id b = -(forceOf b a)
            | otherwise = (pos a - pos b) *:. charge a * charge b / norm (pos a - pos b) ^ 3

where (*:.) :: Vector -> Double -> Vector is vector-scalar multiplication so the above is a 1/r^2 force law. Notice that here we memoize pos a - pos b and then we also memoize forceOf a b for use as forceOf b a.

Then you want to use dvs = [dt * sum (map (forceOf a) particles) / mass a | a <- particles] to get a list of changes in velocity, then zipWith (+) (map velocity particles) dvs

One problem is that this approach doesn't do so well with numerical uncertainty: everything for time $t+1$ is based on things that were true at time $t$. You can start to solve this problem by solving a matrix equation; instead of v₊ = v + dt * M v (where v = v(t) and v₊ = v(t + 1)), you can write v₊ = v + dt * M v₊, so that you have v₊ = (1 − dt * M)^-1 v. That can often be more numerically stable. It is potentially even better to mix the two solutions 50/50. v₊ = v + dt * M (v + v₊) / 2.There are lots of options here.