An elementary direct strategy to prove the result is the following. Let us make things simpler, at least initially, referring to the elementary equation in R (the procedure easily extends to Rn and to Hamilton equations in particular)
dxdt=F(x(t)).
Let us assume that two solutions x and x′ exist with the same initial condition x(0)=x0=x′(0). Finally suppose that dF/dx is bounded by some positive real L<+∞ (if dF/dx is continuous, it is automatically true working in a bounded neighbourhood).
I go to prove that it must be x(t)=x′(t) in a sufficiently small right neighbourhood of t=0.
The considered differential equation is equivalent to the integral one (as Ron said)
x(t)=x0+∫t0F(x(u))du,
so that, we also have
x′(t)=x0+∫t0F(x′(u))du,
and thus
x(t)−x′(t)=∫t0F(x(u))−F(x′(u))du.
Therefore, using Lagrange theorem (below, y(u) is some unknown point between x(u) and x′(u)),
|x(t)−x′(t)|≤∫t0|F(x(u))−F(x′(u))|du=∫t0|dFdx|y(u)(x′(u)−x(u))|du
We can now use the bound |dF/dx|≤L, obtaining
|x(t)−x′(t)|≤L∫t0|x(u)−x′(u)|du≤Ltmaxu∈[0,t]|x(u)−x′(u)|
Obviously, if 0≤t≤T
Ltmaxu∈[0,t]|x(u)−x′(u)|≤LTmaxu∈[0,T]|x(u)−x′(u)|
so that
|x(t)−x′(t)|≤LTmaxu∈[0,T]|x(u)−x′(u)|∀t∈[0,T]
which eventually produces, since the bound must hold also for the value t=t0 where the (continuous) function in the left-hand side above attains its maximum,
maxt∈[0,T]|x(t)−x′(t)|≤LTmaxu∈[0,T]|x(u)−x′(u)|,
i.e.
0≤maxt∈[0,T]|x(t)−x′(t)|≤LTmaxt∈[0,T]|x(t)−x′(t)|.
We can take T>0 so small that 0<LT<1. In this case the found inequality can hold only if
maxt∈[0,T]|x(t)−x′(t)|=0,
and thus x(t)=x′(t) for t∈[0,T].
We have found that, at least in a sufficiently small neighbourhood [0,T] of the initial time (however the procedure can be implemented also for t<0), x(t)=x′(t). This reasoning, gluing together several patches, establishes the uniqueness of the solution on its whole domain.
Let us pass to the true Hamilton equations. The proof is essentially identical for a C2 Hamiltonian H defined an open set of R2n, with the following replacements.
R→R2n
x→→x=(→q,→p)
F→→F(→x):=S∇→xH(→x)
where S is the symplectic 2n×2n matrix with −I, I on the anti diagonal and 0, 0 on the diagonal when viewed as 2x2 block matrix.
|⋅|→||⋅||
the Euclidean norm in R2n.