scientific computing - Solving floating-point rounding issues C++ -

i develop scientific application (simulation of chromosomes moving in cell nucleus). chromosomes divided in small fragments rotate around random axis using 4x4 rotation matrices.

the problem simulation performs hundreds of billions of rotations, therefore floating-point rounding errors stack , grow exponentially, fragments tend "float away" , detach rest of chromosome time passes.

i use double precision c++. soft runs on cpu moment ported cuda, , simulations can last 1 month @ most.

i have no idea how somehow renormalize chromosome, because fragments chained (you can see doubly linked-list), think best idea, if possible.

do have suggestions ? feel bit lost.

thank much,

h.

edit: added simplified sample code. can assume matrix math classical implementations.

// rotate 1000000 times (int = 0; < 1000000; ++i) {     // pick random section start     int istart = rand() % chromosome->length;      // pick end 20 segments further (cyclic)     int iend = (istart + 20) % chromosome->length;      // build rotation axis     vector4 axis = chromosome->segments[istart].position - chromosome->segments[iend].position;     axis.normalize();      // build rotation matrix , translation vector     matrix4 rotm(axis, rand() / float(rand_max));     vector4 oldpos = chromosome->segments[istart].position;      // rotate each segment between istart , iend using rotm     (int j = (istart + 1) % chromosome->length; j != iend; ++j, j %= chromosome->length)     {         chromosome->segments[j].position -= oldpos;         chromosome->segments[j].position.transform(rotm);         chromosome->segments[j].position += oldpos;     } } 

you need find constraint system , work keep within reasonable bounds. i've done bunch of molecular collision simulations , in systems total energy conserved, every step double check total energy of system , if varies threshold, know time step poorly chosen (too big or small) , pick new time step , rerun it. way can keep track of what's happening system in real time.

for simulation, don't know conserved quantity have, if have one, can try keep constant. remember, making time step smaller doesn't increase accuracy, need optimize step size amount of precision have. i've had numerical simulations run weeks of cpu time , conserved quantities within 1 part in 10^8, possible, need play around some.

also, tomalak said, maybe try reference system start time rather previous step. rather moving chromosomes keep chromosomes @ start place , store them transformation matrix gets current location. when compute new rotation, modify transformation matrix. may seem silly, works because errors average out 0.

for example, lets have particle sits @ (x,y) , every step calculate (dx, dy) , move particle. step-wise way this

t0 (x0,y0) t1 (x0,y0) + (dx,dy) -> (x1, y1) t2 (x1,y1) + (dx,dy) -> (x2, y2) t3 (x2,y2) + (dx,dy) -> (x3, y3) t4 (x3,30) + (dx,dy) -> (x4, y4) ... 

if reference t0, this

t0 (x0, y0) (0, 0) t1 (x0, y0) (0, 0) + (dx, dy) -> (x0, y0) (dx1, dy1) t2 (x0, y0) (dx1, dy1) + (dx, dy) -> (x0, y0) (dx2, dy2) t3 (x0, y0) (dx2, dy2) + (dx, dy) -> (x0, y0) (dx3, dy3) 

so @ time, tn, real position have (x0, y0) + (dxn, dyn)

now simple translation example, you're not going win much. rotation, can life saver. keep matrix euler angles associated each chromosome , update rather actual position of chromosome. @ least way won't float away.