I want to simulate a N-body system in the most functional way possible. The system I'm modeling are metallic colloidal particles in a solution of peroxide water, as here http://arxiv.org/pdf/1502.07144.pdf. Due to the small size of the particles, the system is in a low Reynold's regime, so there is no inertia. Furthermore, the surface of the particles react with the fluid and there is creation and consumption of hydrogen (H) on their surface. Gradients on the concentration of H will produce forces on the particles, which depend on their material parameters. The particles can show self-organization and motion.

In general, the long range force between the particles is given by the combined action of the hydrodynamic forces and the gradient of the H concentration field, which is yet to be experimentally characterised. We assume that the force between the particles is a radially symmetric, pair-wise function that decays as $1/r^2$. In principle this force can be modified by changing the chemistry of the system and/or the composition of the particles, say, by adding salts and changing the potential to a Yukawa potential instead of $1/r^2$. I'm interested on finding the structures formed by these colloidal particles that show dynamic function, e.g., oscillators, by modifying the forces or the material parameters.

The systems I study have generally only a few particles, so I'm not so concerned with an efficient implementation of the force computation. Rather I want a very abstract code that allows me to eventually run several thousands of simulations of small systems with different parameters, analyse the structures found, and re-run with slightly different parameters —- think of genetic algorithms. A function that runs fast would be nice though.

This is one example of the structures found for two kind of particles interacting:

enter image description here

I already have a C code that integrates the system and works well, but I find the problem of solving a general N-body problem in Mathematica using its built-in integrators interesting in own right and as an example of how straight-forward such code might be in a functions style. (I took some ideas and code from 38687). Further, the Mathematica code is 100 times shorter (albeit slower) and could be easily expanded to simulate molecules, granular matter or the like by simply changing the intereaction force and leaving all the technical details to a fast implementation to compute the force between the particles. I think this would be great for teaching molecular dynamic simmulations for example.

I've been using Mathematica for some years for data analysis but recently started to study more in depth the functional way of programming, so excuse me if the code is not elegant, it's working though. This is where I am at the moment.

Define some particles, define a certain functional form for the force and a interaction matrix that is not necessarily symmetric. (If you simulate gravity, for example, you have a symmetric matrix, in diffusophoretic systems you generally don't.)

Then integrate with NDSolve.

Clear[force, r, type]
numbodies = 9
pos0 = RandomReal[{-10, 10}, {numbodies, 2}]

type = 1 + Boole[OddQ[#]] & /@ 
  Range[numbodies] (*type of particle, we start with 2 but they could \
be more, keepp it balanced to form structures*)

matrixB = 
 2.51 {{1, -1}, {-.55, 
    1.1}} (*interaction matrix between the particles*)

force[j_] := 
     r[j][t] - r[i][t]]) (matrixB[[type[[i]], type[[j]]]]/
      EuclideanDistance[r[j][t], r[i][t]]^2 +
     1/EuclideanDistance[r[j][t], r[i][t]]^12), {i, 
    j]}];(*diffusophoretic force and leonard jones repulsion*)

odesys = Table[{r[j]'[t] == force[j], r[j][0] == pos0[[j]]}, {j, 
depvars = Flatten[Table[{r[j]}, {j, numbodies}]];
tfin = 500;
sol = NDSolve[odesys, depvars, {t, 0, tfin}][[1]];

pos = Array[r, {numbodies}] /. sol;
  ParametricPlot[#[t] & /@ pos, {t, Max[0, tfin1 - 50], tfin1}, 
   ImageSize -> 400, Frame -> True, PlotRange -> {-20, 20}], 
   MapIndexed[{Hue[.35 type[[First[#2]]] ], Disk[#1[tfin1], 0.5]} &, 
    pos]]], {tfin1, 0.1, tfin, 1}]

So, my questions are:

Is this a proper way to functionally code an N-body integration? Are there better/more-elegant ways of doing it? Can the performance be improved? Is the Animate of the plot re-computing the solution for each final time it shows? Can I pass parameters to the interaction matrix without Mathematica complaining? My objective is to be able to answer questions concerning the interaction matrix and the structures formed, so this code is just a little part of big multi-dimensional minimization, say, to find the matrix/force law/initial-condition that maximizes the rotational energy in the system.

(if you are wondering why the values of the interaction matrix, they allow for the formation of a pure oscillator in some cases.)

Feel free to edit the title of the question if it's not well chosen. I learnt a lot while writing this question so maybe it's not the best phrasing now that I have NDSolve working.


I think I haven't been clear enough with the generality of the problem. This is basically the same code to solve a 3-body gravitational problem.

Clear[force, r, type, v]
numbodies = 3
pos0 = RandomReal[{-10, 10}, {numbodies, 2}]
vel0 = RandomReal[{-.1, .1}, {numbodies, 2}]
vel0 = (# - Mean[vel0]) & /@ vel0
type = 1;
matrix = {{1, 1}, {1, 1}} (*interaction matrix between the particles*)
force[j_] := 
 Sum[-(Normalize[r[j][t] - r[i][t]]) matrix[[type, type]]/(.1 + 
      EuclideanDistance[r[j][t], r[i][t]]^2), {i, 
    j]}];(*gravitational force with a kernel*)

odesys = Table[{v[j]'[t] == force[j],
   r[j]'[t] == v[j][t],
   r[j][0] == pos0[[j]],
   v[j][0] == vel0[[j]]}, {j, numbodies}];

depvars = Flatten[Table[{r[j], v[j]}, {j, numbodies}]];

tfin = 200;
sol = NDSolve[odesys, depvars, {t, 0, tfin}][[1]];

pos = Array[r, {numbodies}] /. sol;

  ParametricPlot[#[t] & /@ pos, {t, Max[0, tfin1 - 50], tfin1}, 
   ImageSize -> 400, Frame -> True, PlotRange -> {-20, 20}], 
  Graphics[MapIndexed[{Hue[.35], Disk[#1[tfin1], 0.5]} &, 
    pos]]], {tfin1, 0.1, tfin, 1}]

And this is how it looks (don't know why the flinchs appear when I export to gif...):

enter image description here

In this case the interaction matrix is just a scalar, and it's where the Gm1m2 should go. Also the particles have only one type. The equation of motion in this case is a bit more complicated since you need equations for both position and velocity, but the way of solving the problem is exactly the same. I hope this clarify a bit my question.

  • $\begingroup$ I have edited the title, but am not certain my edit is a sufficiently better description of your problem that it should stand. I you think it is not, please go back to your original title. I also tried to improve your sentence structure, but in some cases the wording was ambiguous enough that I may have made some bad changes. Please review my 'corrections' and correct them as you see fit. $\endgroup$
    – m_goldberg
    Commented Feb 1, 2016 at 7:22
  • $\begingroup$ @m_goldberg thanks for the edits. I added another example showing how the same approach works in different systems, not only fluid dynamics. $\endgroup$ Commented Feb 1, 2016 at 15:25
  • $\begingroup$ If you rewrite your matrixB so that it's a function of parameter type (e.g., matrixB[1] := ..., matrixB[2] := ...) instead of a list, you can easily modify it later to be dependent on an additional variable, like matrixB[t_, 1] := ... or matrixB[t_][1] := .... $\endgroup$
    – shrx
    Commented Feb 4, 2016 at 9:48
  • $\begingroup$ As per performance improvements, it depends on the size of your system. For many particles you could gain performance by utilising GPU cores with CUDA or OpenCL. see: wolfram.com/mathematica/new-in-8/cuda-and-opencl-support/… $\endgroup$
    – shrx
    Commented Feb 4, 2016 at 9:52
  • $\begingroup$ thanks for the input @shrx but if you look at the code of Nbody.cl is a typical C code for the GPU, which yes is faster but not beautiful. What I'm trying to do now is a tree-code force computation in a functional way, that goes more along the lines of my question. $\endgroup$ Commented Feb 5, 2016 at 4:13

1 Answer 1


In version 12.0 you can use NBodySimulation for that:

data  = NBodySimulation[
   "InverseSquare", {<|"Mass" -> 1, "Position" -> pos0[[1]], 
     "Velocity" -> vel0[[1]]|>,
    <|"Mass" -> 1, "Position" -> pos0[[2]], 
     "Velocity" -> vel0[[2]]|>,
    <|"Mass" -> 1, "Position" -> pos0[[3]], 
     "Velocity" -> vel0[[3]]|>}, tfin];

Visualize the data:

ParametricPlot[Evaluate[data[All, "Position", t]], {t, 0, tfin}]

enter image description here

  • 1
    $\begingroup$ is it possible to write your own force? the example in the documentation only shows the use of "Harmonic" or "InverseSquare", not a user defined function f[statei, statej]. To make it more concrete, say how would one write Lenard Jones? (maybe this should be a separate questionn) $\endgroup$ Commented Apr 27, 2019 at 20:12
  • 7
    $\begingroup$ I have the same questions: 1) how to define your own "PairwiseForce" or "PairwisePotential" function for NBodySimulation 2) How to add "ExternalForce", or friction in such systems? The manual (help) gives no explanation or examples. $\endgroup$
    – DavidG
    Commented Jun 14, 2019 at 18:58

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.