I'm working with an RKF45 integrator that I have programmed using CUDA C++ on my GPU and am pondering a few questions as I'm trying to track down some issues with my code.  I'm using double precision values in the calculation.  Is it feasible of me to demand of my integrator that the difference between the 4th and 5th order estimates be < 1 x 10 ^-16 if machine precision is that of a double?If I do impose local error to be 10E-16 and I integrate over, for example, 20,000 time steps, should I then be expecting a maximal global error of the product of these two numbers : 2E-12?Indeed if that is the case, that's not too bad, however I'm trying to track the source of a leak where my numbers become unrealistic and I'm trying to rule out machine precision as the source of that error.EDITMy RHS of my equations look like this.double frag_term = 0;    double flux = 0;    if (r == ((maxlength)-1))        flux = -km*(r)*conc[r]+2*(ka)*conc[r-1]*conc[0];    else if ( (r > ((nc)-1)) && (r != ((maxlength)-1)) )    {        frag_term = conc2[maxlength-1] - conc2[r];          flux = -(km)*(r)*conc[r] + 2*(km)*frag_term - 2*(ka)*conc[r]*conc[0] + 2*(ka)*conc[r-1]*conc[0];    }    else if (r == ((nc)-1))    {        frag_term = conc2[maxlength-1] - conc2[r];        flux = (kn)*pow(conc[0],(nc)) + 2*(km)*frag_term - 2*(ka)*conc[r]*conc[0];    }    else if (r < ((nc)-1))        flux = 0;where conc is an array of concentrations whose values I am integrating over time, conc2 is an array of partial sums such that conc2[r] = sum of all conc array items less than and including the rth term.
Machine precision and local error
c++;error estimation;time integration;cuda;runge kutta