If a body falls down for one second (due to gravity g) this leads to a position error of 0.5*g which is quite big (in my opinion).
Given a gravity of -10 and one timestep of 1s the error would be indeed 5m, but this is an unrealistic assumption. The position difference scales with the size of the timestep. So for a more realistic assumption of a timestep of 1/60 the absolute position difference is 83 mm after 1s. Also noting that all particles are subject to quite alot constraints this is really negligible. The only thing you want to make sure is that the solver sees the gravity directly and not with one frame delay like for the normal Euler. Yours and Muller's solver do this. Code:
#include <iostream>const float g = -10.0f;const float dt = 1.0f / 60.0f;int main( void ){float x1 = 10.0f;float v1 = 0.0f;float x2 = 10.0f;float v2 = 0.0f;float x3 = 10.0f;float v3 = 0.0f;for ( int i = 0; i < 60; ++i ){//// Euler//// x(t+dt) = x(t) + v(t) * dt// v(t+dt) = v(t) + a(t) * dt//x1 += v1 * dt;v1 += g * dt;// // Semi-implicit Euler//// v(t+dt) = v(t) + a( t ) * dt// x(t+dt) = x(t) + v(t+dt) * dt// v2 += g * dt;x2 += v2 * dt;//// 2. Order Euler//// x(t+dt) = x(t) + v(t) * dt + 0.5 * a(t) * dt^2// v(t+dt) = v(t) + a(t) * dtx3 += v3 * dt + 0.5f * g * dt * dt;v3 += g * dt;printf( "E1: (%.3f, %.3f), E2: (%.3f, %.3f), E3: (%.3f, %.3f)\n", x1, v1, x2, v2, x3, v3 );}return 0;}
