Formalism to deal with discontinuous potentials in classical mechanics (hard wall, hard spheres) It seems to me that Hamiltonian formalism does not suit well for problems involving instantaneous change of momentum, like particle collisions with hard wall or hard sphere gas model. At least I could not apply it straightforwardly to the simplest possible problem of 1D particle hitting a wall:
            │/ wall
            │/
  particle  │/
───o────────┼────────────> x
            │/
            │/
            │/

My attempt was quite direct. I took the Hamiltonian to be
$$H = \frac{p^2}{2m} + U(x)$$
with potential $U$ defined as
$$U(x) = \cases{ 0, \; x < 0, \\[.5em] K, \; x > 0, \\[.5em] E, \; x= 0.}$$
where $E$ is the particle energy and $K > E$. Hamiltonian equations should read as
$$\cases{\dot x = \frac{p}{m}, \\[.6em] \dot p = - \frac{d U}{d x}.}$$
It is not hard to integrate the first equation, but my attempts to integrate the second one did not lead to any meaningful result (that's why I do not share them here, it was a complete failure).
So I ask whether it is possible to obtain the solution to the problem by directly integrating Hamiltonian equations in the form above, without relying on general mechanical theorems/principles like energy conservation? Or is such an approach completely unsuitable for the task?
If so, what is the general (and elegant) approach to such systems?
There exist a related question on PSE "Hamiltonian function for classical hard-sphere elastic collision", but the setting is more cumbersome.
 A: I think I got how to deal with the problem in a straightforward way, without the passage to a limit.
Let the phase space of the initial problem be the half plane
$$\{ \,(q,p) \; | \; q > 0\},$$
the wall is at $q=0$. In this phase space when the particle reaches the point $(0, p)$ it instantaneously teleports to the point $(0, -p)$. Particle trajectories are thus discontinuous.
The trick now is to glue the half plane into a cone, so that particle trajectories will become continuous. At the same time the particles are considered free on the entire trajectory, yielding the Hamiltonian being just a kinetic energy. Sorry for the lack of appropriate drawings, but I hope it is not so hard to imagine. The operation can be formally achieved with a non-canonical change of coordinates to $(r, \varphi)$:
$$\begin{cases}
  p = 2 r \sin \dfrac{\varphi}{2}, \\[.5em]
  q = 2 r \cos \dfrac{\varphi}{2}.
\end{cases}$$
These are actually polar coordinates in a plane perpendicular to the cone axis. The symplectic form $dp \land dq$ transforms to $2r \, d \varphi \land dr$, or, in other words, the Poisson matrix becomes
$$
\begin{bmatrix}
0 & \dfrac{1}{2r} \\
-\dfrac{1}{2r} & 0
\end{bmatrix}
$$
The Hamiltonian (kinetic energy) is given by
$$H = \frac{2}{m} r^2 \sin^2 \frac{\varphi}{2}.$$
All of the above leads to the Hamiltonian flow of
$$X_H = \frac{1}{2m} r \sin \varphi \frac{\partial}{\partial r} - \frac{1 - \cos \varphi}{m} \frac{\partial}{\partial \varphi}$$
and the equations of motion of
$$\begin{cases}
  \dfrac{d r}{dt} = \dfrac{1}{2m} r \sin \varphi, \\[.5em]
  \dfrac{d \varphi}{dt} = -\dfrac{1}{m} (1 - \cos \varphi).
\end{cases}$$
These are amenable to a relatively simple integration, free of the subtleties of special functions. As result one gets as a solution
$$\begin{gather}
  \cot \dfrac{\varphi(t)}{2} = C_1+ \dfrac{t}{m}, \\[.5em]
  r(t) = C_2 \sqrt{1 + \cot^2 \dfrac{\varphi(t)}{2}},
\end{gather}$$
which also can be obtained directly from the known solution in $(q,p)$ half-plane and coordinate transformation rules.
To sum up, the effect of the wall is accounted for via the change of the phase space topology. Thus, particles are considered free, with a Hamiltonian being a kinetic energy which is preserved. As far as I reckon the phase space is not a cotangent bundle of a configuration space anymore. If this is true along with all the above derivation, this represents probably the most simple case of a phase space that is not a cotangent bundle. 

Further investigation
Although at first I was guided by geometrical reasoning about the phase space, now I have been thinking about the coordinate transformation itself. I came up with another transform, which is much closer to the original coordinates:
$$\begin{cases}
  p = \operatorname{sgn} \! \left(\, \widetilde q \,\right) \, \widetilde p, \\[.5em]
  q = \left|\, \widetilde q \, \right| .
\end{cases}$$
Here it is assumed that $\frac{d}{dx} \operatorname{sgn} x = 0$. Then $dp \land dq = d \, \widetilde p \land d \, \widetilde q$ and the Hamiltonian is the one of a free particle:
$$H = \frac{\; \widetilde p^{\, 2} \!}{2m}.$$
Like with trigonometric coordinate transformation there seems to be a need to choose the right branch for the transform $q \mapsto \widetilde q$, but since $q \geqslant 0$ there is no ambiguity.
A: Perhaps the simplest and most intuitive approach is to regularize the hard wall potential
$$V_0(x)~=~\left\{ \begin{array}{rcl} 0 &\text{for}& x<0 \cr\cr  \infty &\text{for}& x>0\end{array}\right. $$
as 
$$ \lim_{\varepsilon \to 0^+} V_{\varepsilon}(x) ~=~V_0(x).$$
For instance, one could choose the regularized potential as 
$$  V_{\varepsilon}(x)~=~\frac{x}{\varepsilon}\theta(x).$$
This corresponds to constant velocity for $x<0$ and constant acceleration for $x>0$. Next write down a continuous solution for the position $x_{\varepsilon}(t)$ as a function of time $t$, say, for given pertinent initial conditions.
Finally, at the end of the calculation, one should remove the regularization $\varepsilon \to 0^+$ again, and check if the limit $$ \lim_{\varepsilon \to 0^+} x_{\varepsilon}(t)$$ makes physical sense.
A: I will yet add another formalism:
Lets start with the hamiltonian form of Hamilton's Principle. Let  $c: \mathbb R \longrightarrow T^*Q; t\mapsto (q(t),p(t))$ be the trajectory of a particle in the phase space of the configuration space $Q$, we define a subset of $Q$, $C$, where no contac occur between the particles, and $\partial C$ is te set of ponts wer contac has occure, but not penetration. We assume that $c$ will intersect $\partial C$ at time $t_c$, that is $q(t_c)\in \partial C$. $c$ is smooth out of $t_c$. Lets now derive the equations of motin with variations. We start with the action in hamiltoninan formalism with Lagrangian $L = \dot q p -H(q,p)$:
$$
S[c] = \int_0^T L(q,\dot q, p)dt = \int_0^{t^\lambda} L(q,\dot q, p)dt + \int_{t^\lambda} ^T L(q,\dot q, p)dt , 
$$
and with the variation on the path denoted by $c^\lambda$ and the variation on the collision instant by $t^\lambda$, if we set $\delta S =0$, by the Leibniz integration rule: 
$$
\frac{d}{d\lambda}S[c^\lambda] =  \int^{t^\lambda}_0 \left[\frac{\partial L}{\partial c}\delta c +\frac{\partial L}{\partial \dot c}\delta \dot c \right]
dt +  \int^{T}_{t^\lambda} \left[\frac{\partial L}{\partial c}\delta c +\frac{\partial L}{\partial \dot c}\delta \dot c \right]
dt + \left.L\delta t^\lambda\right|_{t^\lambda_-} -L\delta $$
t^\lambda\right|{t^\lambda+}  $$
We integrate by parts as usual and require that $\delta q(0)=\delta q(T)=0$, there is no need to $\delta p$ to vanish as it is not integrated by parts, we get:
$$
  \int^{t^\lambda}_0 \left[\frac{\partial L}{\partial c}\delta c -\frac{d}{dt}\frac{\partial L}{\partial \dot c}\delta c \right]
dt +  \int^{T}_{^\lambda} \left[\frac{\partial L}{\partial c}\delta c -\frac{d}{dt}\frac{\partial L}{\partial \dot c}\delta c \right]
dt - \left[\frac{\partial L}{\partial \dot c}\delta c + L\delta t^\lambda\right]_{t_-^\lambda}^{t^\lambda_+},
$$
as we require the variations to vanish at all time, we can focus only in the third term:
$$
\left[\frac{\partial L}{\partial \dot c}\delta c + L\delta t_c\right]_{t_c^-}^{t_c^+} = \left[p \delta q + (\dot q p - H)\delta t_c\right]_{t_c^-}^{t_c^+}.
$$
We recall that $q(t_c)\in \partial C$,  so the variation $\delta q(t_c) + \dot q(t_c)\delta t_c \in T_{q(t_c)}Q$ and we get:
$$
\left[p \delta q + (\dot q p - H)\delta t^\lambda\right]_{t_-^\lambda}^{t_+^\lambda} = \left[p(\delta q + \dot q p\delta t^\lambda)  - H \delta t^\lambda \right]_{t_-^\lambda}^{t_+^\lambda} =0,
$$
considering the variations arbitrary, we get conservation of energy, $H$ has to  be constant and we get that the change of momentum $p(t_+) -p(t_-)$ has to be normal to $T_{q(t_c)}Q$, and the conservation of energy fixes the magnitude. The major cave beat is that it is necessary to interpret which direction the particle exits the collision, because it is not possible to extract the information from the equations
Sources:
http://thesis.library.caltech.edu/1934/1/01thesis.pdf
http://people.math.sfu.ca/~van/papers/40603.pdf
A: The problem here is not about a correct "physical" formalism. The problem is that we are seeking singular solutions. 
Let me explain (however, using Newtonian formalism, just because I'm more comfortable with it).  
When talking about an initial value problem of an ODE of the form: $\ddot{x}(t) = \frac{d}{dx}E(x(t))$ one should have a uniform Lipschitz continuity of the right hand side to have a unique local classical solution according to the Picard-Lindelöf theorem. Clearly, in your case of a delta potential (i.e. hardwall) this is not the case since we have $\ddot{x}(t)=-\delta^{'} (x(t))$. 
A common practice among mathematicians in such a case is to seek weak solutions (e.g. in the sense of Schwartz-Distribution theory). Instead of looking for a smooth solution $x(t)$, one looks for less regular solutions. This is in perfect accordance with the physical insight since we are expecting particles that have instantaneous change of momentum. To write the weak formulation of the problem we define $\mathcal{D}(\mathbb{R})$ to be the space of smooth functions with compact support and $\mathcal{D}^{'}(\mathbb{R})$ to be its dual space. The weak formulation reads: find $x \in \mathcal{D}^{'}(\mathbb{R})$ such that:
$$  \int \ddot{x} (s) \phi(s) \; ds = - \int \delta^{'} (x(s)) \phi(s) \; ds $$ 
where $\phi \in \mathcal{D}(\mathbb{R})$. Now, using integration by parts, one can relax the conditions on $x(t)$. A solution to the last equation is said to be a weak solution of our initial problem. However, here we have a problem resulting from the linear nature of Schwartz' distribution theory. The composition of a distribution $\delta ^{'}$ with an element $x(t)$ is only defined when $x(t)$ is smooth. In other words, we can not seek a weak solution $x(t)$ and at the same time consider the potential to be a distribution.  
A proper and rigorous way to deal with this problem is using a nonlinear theory of distributions, where the composition of a distribution is defined, even for less regular objects. This is possible in Colombeau's algebra of tempered generalised functions.  
In "Geometric Theory of Generalised Functions with Applications to General Relativity" by Michael Grosser et al. Newton's equation with a delta potential was proven to have a unique solution in Colombeau's algebra. This solution is parametrised by a regularisation parameter $\varepsilon$. Moreover, it's been proven that when taking the limit of the solution as $\varepsilon \to 0$ it converges to the solution one expects in physics (i.e. the trajectory experiences a reflection at the potential).  
Note that this is very similar to the regularisation procedure described in one of the answers. However, in this algebra the regularisation is made more rigorous and consistent with distribution theory.     
A: Try using $U(x) = \Theta(x)$, $\Theta(x) \equiv \begin{cases} 1, & x>0\\ 0, & x<0\end{cases}$ . The force becomes $\dot{p} = - \delta(x)$, as anticipated. Because the force is infinite with direction opposite that of the momentum, the particle cannot cross the plane $x=0$.
Another way to see that the particle won't cross $x=0$ is to integrate the equation of motion over $x$ just before and just after $x=0$:
$$ -1 = m\int_{0_-}^{0_+} \dot{v} dx = \frac{m}{2} (v^2(0_+)- v^2(0_-)),$$
from which we obtain $v(0_+) = 0$, because the energy is $1$.
