I was working through some formulas with this recently, and it turns out that using this approximation still gives a physically-correct trajectory, for a slightly different initial velocity.
ie. We'd expect a body starting at position \$\vec p_0\$ with initial velocity \$\vec v_0\$ under constant acceleration \$\vec a\$ (eg. due to gravity) to follow the parabola...
$$\vec p(t) = \vec p_0 + \vec v_0 \cdot t + \frac {\vec a} 2 \cdot t^2$$
But using the symplectic Euler integration with a fixed timestep \$\Delta t\$ in eg. Box2D, we instead see the parabola...
$$\vec p(t) = \vec p_0 + (\vec v_0 + \frac {\vec a} 2 \Delta t) \cdot t + \frac {\vec a} 2 \cdot t^2$$
This is still a physically-plausible trajectory (ie. given the motion observed on frame 1 & 2, all subsequent frames correctly match that observation), just with a slightly different boundary condition. And the error shrinks with smaller timesteps. The trick is that only someone with access to the "true" under-the-hood value of \$v_0\$ could tell the difference. Any other observer would simply conclude that the velocity at the start was different, and the rest of the motion proceeds correctly from there.
Given that the initial velocities are usually either:
fully controlled by the game developer (eg. specifying a jump impulse or launch velocity of a projectile), it's trivial for the developer to subtract a half-timestep of acceleration from the input up front to get back on the intended curve, without complicating the integrator.
extremely approximate anyway (eg. the bounce resulting from a collision resolved with discrete overlap detection & restitution along the shortest penetration, rather than true continuous material simulation)
In practice, a realtime physics engine sweeps all kinds of things under the rug, like the precise moment within the frame when an impulse was applied (if it came near the start of the interval, then using the average velocity might make your simulation less accurate, not more) so it's not necessarily in the business of "true reproduction," just "plausible simulation"
...so the designers of these physics engines may have felt that this small approximation of the initial conditions (only half a timestep's worth of error) was reasonably harmless, in the contexts in which these engines tend to do their work.
Granted, this decision may have been made back in a day when processors were less powerful, and the simplification of an extra vector multiply & add out of the inner loop of the physics engine was significant. That might well not be the bottleneck today, but established physics engines could be wary to change the core integrator, for fear of breaking software, shared assets, and libraries that depend on the old behaviour.
If you're interested in experimenting with using the average velocity over the frame, instead of treating all acceleration as occurring at the beginning, it's worth a try. It's not physically wrong, but you might or might not find much benefit in terms of the visual plausibility of the simulation that results.
s = (v + 0.5*a*dt) * dt
, to get the average velocity over the interval? If we useds = 0.5 * v * dt
we'd only ever move at half the velocity, even when our motion is unchanging. \$\endgroup\$