Possiblity of predicting behavior of system when properties described as functions over space? Suppose I am given a system that consists of a distribution of charged particles(which are all over space and are point-charges). They are described by a set of functions instead of variables. These functions are:


*

*$C(x)$ - a function describing the amount of charge at each point in space. For example, take the one dimensional example of $C(x) = \sin(x)$.

*$M(x)$ - a function describing the amount of mass at each point in space.

*$V(x)$ - the function describing what the sum of all the velocities of the particles at a point is.


Given these are the functions to describe the system at time $t = 0$, is it possible to predict the state of this system (the charge distribution, velocities and mass distributions) at a later time, say, $t = 2$, using Newtonian dynamics? 
My system can be thought of as a fluid, where distribution of mass and charge vary. I ignore charges if that simplifies things. I want to take into account gravity, though, at least. How do I exactly proceed to find the distribution functions at $t = 2$?
 A: If you take a particle approach, then you might want to investigate Particle-in-cell methods. In this particular approach, each point-charge is given a position, velocity, and charge in any dimension (1d, 2d, or 3d) and you evolve the particles using the Lorentz force law (1) and the position-derivative (2):
\begin{align}
m\frac{d\mathbf v}{dt}&=q\mathbf E+q\mathbf v\times\mathbf B \\
\frac{d\mathbf x}{dt}&=\mathbf v
\end{align}
What complicates this is that, in order to compute $\mathbf E$ and $\mathbf B$, you need to use the fact that you have a discrete collection of particles, with charge density $\rho(\mathbf x)$, and, from this, you need to determine the potential, $\phi(\mathbf x)$, and then the associated $\mathbf E$ & $\mathbf B$ fields. State-of-the-art simulations employing this are able to store information of about $10^9$ particles (note that this is on a supercomputing cluster, not on your average desktop).
If you want the fluid approach, then ACuriousMind is correct and you want to investigate Magnetohydrodynamics (MHD). These modify the Euler hydrodynamic equations to incorporate the dynamic effects of a magnetic field:
\begin{align}
\frac{\partial \rho}{\partial t}&=-\nabla\cdot\rho\mathbf v \\
\frac{\partial\rho\mathbf v}{\partial t}&=-\nabla\cdot\left(\rho\mathbf v\mathbf v+(p+\frac12B^2)\mathbb I\right) -\mathbf B\mathbf B\\ 
\frac{\partial e}{\partial t}&=-\nabla\cdot\left(\left[e+p+B^2+\mathbf B\mathbf B\right]\cdot\mathbf v\right)\\
\frac{\partial\mathbf B}{\partial t}&=-\nabla\cdot\left(\mathbf v\mathbf B-\mathbf B\mathbf v\right)
\end{align}
where the $\mathbf a\mathbf b$ is a dyadic, $\rho e=p/(\gamma-1)$ with $\gamma$ the adiabatic index is the equation of state, $\mathbb I$ the identity matrix, and $\rho$ the mass density, all other variables take their normal meaning. We also have the further condition that $\nabla\cdot\mathbf B=0$
In order to solve the MHD equations, you need to discretize the domain (break it up into discrete chunks and apply averaged information at "cell-centers"). The divergenceless magnetic field condition is also tough to enforce on a discrete grid, many papers have been devoted to maintaining this to machine precision!  If you want to add in gravity, you need to add an extra term, $\rho\mathbf g$, to the momentum equation ($\partial_t\rho\mathbf v$ term)--this is often easier using the potential $\nabla^2\phi=\rho\mathbf g$, but also requires more work!.


If you actually do want to model this, note that you will likely be spending the next several months writing and testing the program to do whichever of the two above you've chosen. 
The alternative is to go onto the github and find an MHD or PIC program that someone else spent months/years writing & testing and use that. I cannot stress enough that you will also need to spend a long time (probably on the order of weeks) studying how the code works (because it is not a black box) so that you can maybe trust the results.
A: An approach can be: find the force on position $x$ due to influence of mass on $x'$ all other:
$$f_g(x,x',t) = \frac{G M(x,t)M(x',t)}{\left( x-x' \right)^2} $$
The integral of this equation, which is the resulting force on $x$ from the influence of all other massess, will be equal to momentum change in this point:
$$ M(x,t) \int \frac{G M(x',t)}{\left( x-x' \right)^2}dx' = \frac{d(M(x,t)*V(x,t))}{dt} $$
$$ M(x,t) \int \frac{G M(x',t)}{\left( x-x' \right)^2}dx' = V(x,t)\frac{dM(x,t)}{dt} + M(x,t) \frac{dV(x,t)}{dt} $$
This is the equation you would need to solve, in this approach. My suggestion in general would be to solve it numerically :).
However, in your case where your initial conditions are periodical and assuming all initial velocities are null, I can intuitively say some features of the solution:


*

*Every point of a maximum of M should not move since the resulting force on it is null due to symmetry.

*All mass around this max, and up to the null-mass points should be attracted to the max.

*Your final states should lead to delta functions located where the max were initially.
