the Lagrangian of a particle system: , where is the particle specific internal energy. The Euler–Lagrange equation of variational mechanics reads, for each particle: When applied to the above Lagrangian, it gives the following momentum equation: where the chain rule has been used, since depends on , and the latter, on the position of the particles. Using the thermodynamic property we may write Pluging the SPH density interpolation and differentiating explicitly leads to which is the SPH momentum equation already mentioned, where we recognize the operator. This explains why linear momentum is conserved, and allows conservation of angular momentum and energy to be conserved as well. Time integration From the work done in the 80's and 90's on numerical integration of point-like particles in large accelerators, appropriate time integrators have been developed with accurate conservation properties on the long term; they are called symplectic integrators. The most popular in the SPH literature is the leapfrog scheme, which reads for each particle : where is the time step, superscripts stand for time iterations while is the particle acceleration, given by the right-hand side of the momentum equation. Other symplectic integrators exist (see the reference textbook). It is recommended to use a symplectic (even low-order) scheme instead of a high order non-symplectic scheme, to avoid error accumulation after many iterations. Integration of density has not been studied extensively (see below for more details). Symplectic schemes are conservative but explicit, thus their numerical stability requires stability conditions, analogous to the Courant-Friedrichs-Lewy condition (see below). Boundary techniques In case the SPH convolution shall be practiced close to a boundary, i.e. closer than , then the integral support is truncated. Indeed, when the convolution is affected by a boundary, the convolution shall be split in 2 integrals, where is the compact support ball centered at , with radius , and denotes the part of the compact support inside the computational domain, . Hence, imposing boundary conditions in SPH is completely based on approximating the second integral on the right hand side. The same can be of course applied to the differential operators computation, Several techniques has been introduced in the past to model boundaries in SPH. Integral neglect The most straightforward boundary model is neglecting the integral, such that just the bulk interactions are taken into account, This is a popular approach when free-surface is considered in monophase simulations. The main benefit of this boundary