[ODE] ODE spatial inconsistency

[pilat@cpsc.ucalgary.ca]

I've run into a problem with ODE.  I simulate simple shapes joined with
universal joints.  When i create my shapes with the same position, they
have exactly the same behaviour.  Problem is, when I create them in a
different position, the behaviour can very greatly - even if the
difference in position is only 0.01 instead of 0.0 on one of the axes. 
I've thought that my simulation might be broken, but after many hours
spent on testing various things, I isolated the problem to ODE having
inconsistent behaviour.

I've created a simple sample that demonstrates the problem. I'll paste the
sample code below if anybody is interested. Using startx=0.0 the result is
consistently (-0.204035, -0.00123104, 0.463304), however setting it to
0.01 produces (-0.19469, 0.000128399, 0.4633) consistently.  I understand
the change in the x coordinate (although it does not vary by 0.01 in the
result) but why the change in the y & z?  The differences in the sample
might not seem large but they are quite large in a more complex setting
and/or with different parameter values.  Could the problem be due to the
contact settings?

Anybody have any idea why this is the case and perhaps how to fix it? 
I've tried changing step size, ERP, CFM, and they all produce a different
result but there are still inconsistencies between the two results (with
x=0.0 and x=0.01).

Thanks,
Marcin

---
#include <iostream>
#include <ode/ode.h>

dWorldID	_world;			// ODE world object
dSpaceID	_space;			// default ODE space
dJointGroupID _contacts;	// contact group
dGeomID	_ground;										// ground geom

void collisionCallback(void *data, dGeomID o1, dGeomID o2)
{
    // maximum 20 contacts
	dContact contacts[20];
    int num = dCollide (o1,o2,20,&contacts[0].geom,sizeof(dContact));

    for (int i=0; i<num; i++)
		{
			contacts[i].surface.mode = dContactBounce;

			contacts[i].surface.mu = dInfinity;
			contacts[i].surface.bounce = 0.1;
			contacts[i].surface.bounce_vel = 0.1;

			dJointID c = dJointCreateContact (_world, _contacts, &contacts[i]);
			dJointAttach (c, dGeomGetBody(contacts[i].geom.g1),
dGeomGetBody(contacts[i].geom.g2));
    }
}

int main (int argc, char * const argv[])
{
    // initialize ode
	_world = dWorldCreate();
	_contacts = dJointGroupCreate(0);
	_space = dHashSpaceCreate(0);
	_ground = dCreatePlane(_space, 0,0,1,0);

    double _stepSize = 0.01;

	// set some rules
	dWorldSetERP (_world, 0.015);
	dWorldSetCFM (_world, 0.01);
	dWorldSetGravity (_world,0,0, -3.0);

    double startx = 0.0;

    // create bodies
	dBodyID b1 = dBodyCreate(_world);
    dMass m1;
	dMassSetBox(&m1, 1.0, 1.0, 1.0, 1.0);
    dMassAdjust(&m1, 1.0);
	dBodySetMass(b1, &m1);
    dBodySetPosition(b1, startx, 0.0, 50.0);

    dGeomID g1 = dCreateBox(_space, 1.0, 1.0, 1.0);
	dGeomSetBody(g1, b1);

    dBodyID b2 = dBodyCreate(_world);
    dMass m2;
	dMassSetBox(&m2, 1.0, 1.0, 1.0, 1.0);
    dMassAdjust(&m2, 1.0);
	dBodySetMass(b2, &m2);
    dBodySetPosition(b2, startx+0.5, 0.0, 50.0);

    dGeomID g2 = dCreateBox(_space, 1.0, 1.0, 1.0);
	dGeomSetBody(g2, b2);

    // create joint
	dJointID j = dJointCreateUniversal(_world, 0);
    dJointAttach(j, b1, b2);
	dJointSetUniversalAnchor(j, startx+0.5, 0.0, 50.0);
	dJointSetUniversalAxis1(j, 0.0, 1.0, 0.0);
	dJointSetUniversalAxis1(j, 0.0, 0.0, 1.0);

//    dJointSetUniversalParam(j, dParamVel, 4.0);
//    dJointSetUniversalParam(j, dParamFMax,4.0);

    for (int i=0; i<10000; i++)
    {
        dSpaceCollide(_space, 0, &collisionCallback);
        dWorldStep(_world,(dReal)_stepSize);
        dJointGroupEmpty(_contacts);
    }

    const dReal * const pos1 = dBodyGetPosition(b1);
    std::cout << "b1 end pos: " << pos1[0] << ", " << pos1[1] << ", " <<
pos1[2] << std::endl;

    return 0;
}