Verlet as given by Wikipedia:
- set $\vec x_1=\vec x_0+\vec v_0\,\Delta t+\frac12 A(\vec x_0)\,\Delta t^2$
- for ''n=1,2,...'' iterate $\vec x_{n+1}=2 \vec x_n- \vec x_{n-1}+ A(\vec x_n)\,\Delta t^2.$
S.I. Euler as given by Wikipedia:
$v_{n+1} = v_n + g(t_n, x_n) \, \Delta t\\ x_{n+1} = x_n + f(t_n, v_{n+1}) \, \Delta t$
Starting with S.I. Euler (and simplifying the notation):
$$ x_{n+1} = x_n + hv_{n+1} \implies v_n = \frac{x_n - x_{n-1}}{h}\\ x_{n+1} = x_n + hv_{n+1} = x_n + h(v_n + ha_n) = x_n + h\left(\frac{x_n - x_{n-1}}{h} + ha_n\right) = 2x_n - x_{n-1} + h^2a_n $$
OK so they're the same thing, cool. That matches my tests (some simple projectile stuff + wind resistance + weird acceleration functions) where Verlet and SI Euler perform to within floating point error. However the Verlet page says, "The global error of all Euler methods is of order one, whereas the global error of [Verlet] is, similar to the midpoint method, of order two." It also says, "The global error in position, in contrast, is $O(\Delta t^{2})$ and the global error in velocity is $O(\Delta t^{2})$." The SI Euler page says "The semi-implicit Euler is a first-order integrator, just as the standard Euler method." How can they have different orders when they are the same method and seem to produce identical results?
EDIT: TL:DR Jorge's question was "I thought the symplectic Euler was O[1], how is it producing an O[2]?" and my answer is "well because you've combine a symplectic Euler update with its time shifted (aka adjoint) counterpart. So you're not really looking at the symplectic Euler anymore."
#---------- Extended Answer -----------#
This was a very interesting question. I think in your derivation you have done a thing that makes your "symplectic Euler" NOT a symplectic Euler. The Velocity Verlet, is known to be the composition of 2 symplectic Euler methods. Knowing this helped me recognize what had happened in your proof. (See "Geometric Numerical Integration" by Hairer - ChVI theorem 3.4)
Let's start with your definition of symplectic Euler method $$v_{n+1} = v_n + ha_n$$ $$x_{n+1} = x_n + hv_{n+1}$$
In your derivation you do the following \begin{align} \text{line 1: } && x_{n+1} &= x_n + hv_{n+1} \\ \text{line 2: }&& &= x_n + h(v_n + ha_n) \\ \text{line 3: }&& &= x_n + h\left(\frac{x_n - x_{n-1}}{h} + ha_n\right) \\ \text{line 4: }&& &= 2x_n - x_{n-1} + h^2a_n \end{align}
Moving from Line 2 to line 3 you insert a $v_n$. The velocity has an explicit update $v_{n+1} = v_n + h a_n$. However to arrive at $v_n$ you use do NOT use $v_n = v_{n-1} + h a_{n-1}$. Rather you reconfigure the implicit position $x_{n+1}=x_n+hv_{n+1}$. So you aren't using the same mapping $\phi$ from $v_{n-1} \rightarrow v_n$ as you are from $v_{n} \rightarrow v_{n+1}$. In equations this looks like, $\phi_a(v_{n-1}) = v_n$ and $\phi_b(v_n) = v_{n+1}$. Where, the mappings are related as $\phi^*_a = \phi_b$ so you (unknowingly) constructed a Verlet update $\phi^*_a \circ \phi_a$.