# Why am I getting so much error for my Runge Kutta Fehlberg solver?

My current project is a reprogramming of a protein folding model involving the solution of thousands of ODEs in C++. I've been making some stop and start progress as I'm writing the solver to run totally on the GPU. I finally got it integrating but when I try to solve dC/dt=-C using the fifth order solution from the RKF45 algorithm with a fixed step size of h=.001, I get decision from the calculator value of e^{-t} on the order of 10^{-4} when is expect global error on the order of 10^{-12}.

Is it because I'm not using adaptive step sizes and error control using both the fourth and fifth order solutions that this is occurring? My thought was that since an adaptive algorithm looks at the difference between the fourth and fifth order guesses and basically treats the fifth order as the "correct answer"that I'd just use that answer to try out the integrator.

EDIT :

I thought that since

Global Error = (Number of Points) * (Local Error)

and

Local Error = h^O

where h is the step size and O is the order, I calculate

O = ln(Global Error/Number of Points)/ln h

Doing an order analysis both the way I thought correct and the way suggested in the answer, I get the following :

Number of Points    h       Global Error    My Way      The Answer's Way
10                  0.1     2.89E-06        6.539E+00
20                  0.05    7.09E-08        6.495E+00   5.350E+00
40                  0.025   1.76E-07        5.216E+00   -1.315E+00
50                  0.02    1.128E-07       5.089E+00   2.003E+00
80                  0.0125  4.401E-08       4.866E+00   2.002E+00
100                 0.01    2.816E-08       4.775E+00   2.001E+00
500                 0.002   1.126E-09       4.316E+00   2.000E+00
1000                0.001   2.810E-10       4.184E+00   2.003E+00
10000               0.0001  3.000E-12       3.881E+00   1.972E+00


So, despite my quasi breakthrough last night, where I did indeed have a programming error, the order of my solver is still off. The aforementioned error associated with the .001 stepsize has decreased dramatically, however, from O(10^-4) to O(10^-10).

2nd EDIT

Here's a graph of LN(E)-LN(a) vs. LN(h). I tossed on a linear fit so I'd expect that to reveal the order to me. What's really confusing now is that I have three calculations for the order of my method. The graph is telling me ~3, Godric's method tells me almost exactly 2 and my method shows it changing from 6->4 for increasing step size. I am more inclined to believe my calculations because the Global Error is so low (~10^-10) and closest to what I'd expect (10^-15*10^3 = 10^-12) for that stepsize (10^-3). EDIT 3:

My Code. Pardon the noobishness.

#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <iostream>
#include <iomanip>                      //display 2 decimal places
#include <math.h>
using namespace std;

/*
Things to do log:
1. Decide how adaptive step sizes is handled by the integrator loop.  We can no longer say how many points to integrate if that is to be determined in the kernel
itself, only what the initial step size is.
2. Figure out how I'm going to graph all this stuff.
3. On paper solve for the step size.
a. arrMin kernel to get min of array
b. Decide how to implement tolerances
5. Use multiples of 256 for maxlength

Important things accomplished:
1. Can create 2D array for storing concentrations over time.
2. Can store a 1D array in 2D array.
3. Can calculate all the k-arrays given static step size.
4. Can take a step given static step size.
5. Can step through an integration from t0 to tf
6. Can accomodate fixed time-step
7. Calculated the current order of integrator using logs.
8. Learned how to profile the program.
9. Fixed the GPU timeout problem with WDDM.

Ideas:
1. Just use large fixed size array
*/
__global__ void rkf5(double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*, int*, size_t, double*, double*, double*);
__global__ void calcK(double*, double*, double*);
__global__ void k1(double*, double*, double*);
__global__ void k2(double*, double*, double*);
__global__ void k3(double*, double*, double*);
__global__ void k4(double*, double*, double*);
__global__ void k5(double*, double*, double*);
__global__ void k6(double*, double*, double*);
__global__ void arrSub(double*, double*, double*);
__global__ void arrMult(double*, double*, double*);
__global__ void arrInit(double*, double);
__global__ void arrCopy(double*, double*);
__device__ void setup(double , double*, double*, double*, double*, int*);
__device__ double flux(int, double*) ;
__device__ double knowles_flux(int, double*);
__device__ void calcStepSize(double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, double*, int*);
__global__ void storeConcs(double*, size_t, double*, int);
__global__ void takeFourthOrderStep(double*, double*, double*, double*, double*, double*, double*);
__global__ void takeFifthOrderStep(double*, double*, double*, double*, double*, double*, double*, double*);

//Error checking that I don't understand yet.
#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, char *file, int line, bool abort=true)
{
if (code != cudaSuccess)
{
fprintf(stderr,"GPUassert: %s %s %d\n", cudaGetErrorString(code), file, line);
if (abort) exit(code);
}
}

//Main program.
int main(int argc, char** argv)
{
//std::cout << std::fixed;          //display 2 decimal places
//std::cout << std::setprecision(12);   //display 2 decimal places
const int maxlength = 1;            //Number of discrete concentrations we are tracking.
double concs[maxlength];            //Meant to store the current concentrations
double temp1[maxlength];                //Used as a bin to store products of Butcher's tableau and k values.
double temp2[maxlength];                //Used as a bin to store products of Butcher's tableau and k values.
double tempsum[maxlength];          //Used as a bin to store cumulative sum of tableau and k values
double k1s[maxlength];
double k2s[maxlength];
double k3s[maxlength];
double k4s[maxlength];
double k5s[maxlength];
double k6s[maxlength];
const int numpoints = 40;
double to = 0;
double tf = 1;
//double dt = static_cast<double>(.5)/static_cast<double>(64);
double dt = (tf-to)/static_cast<double>(numpoints);
double mo = 1;
double concStorage[maxlength][numpoints];   //Stores concs vs. time

//Initialize all the arrays on the host to ensure arrays of 0's are sent to the device.
//Also, here is where we can seed the system.
std::cout<<dt;
std::cout<<"\n";
concs=mo;
std::cout<<concs;
std::cout<<" ";
for (int i=0; i<maxlength; i++)
{
for (int j=0; j<numpoints; j++)
concStorage[i][j]=0;
concs[i]=0;
temp1[i]=0;
temp2[i]=0;
tempsum[i]=0;
k1s[i]=0;
k2s[i]=0;
k3s[i]=0;
k4s[i]=0;
k5s[i]=0;
k6s[i]=0;
std::cout<<concs[i];
std::cout<<" ";
}
concs=mo;
std::cout<<"\n";

//Define all the pointers to device array memory addresses. These contain the on-GPU
//addresses of all the data we're generating/using.
double *d_concs;
double *d_temp1;
double *d_temp2;
double *d_tempsum;
double *d_k1s;
double *d_k2s;
double *d_k3s;
double *d_k4s;
double *d_k5s;
double *d_k6s;
double *d_dt;
int *d_maxlength;
int *d_numpoints;
double *d_to;
double *d_tf;
double *d_concStorage;

//Calculate all the sizes of the arrays in order to allocate the proper amount of memory on the GPU.
size_t size_concs = sizeof(concs);
size_t size_temp1 = sizeof(temp1);
size_t size_temp2 = sizeof(temp2);
size_t size_tempsum = sizeof(tempsum);
size_t size_ks = sizeof(k1s);
size_t size_maxlength = sizeof(maxlength);
size_t size_numpoints = sizeof(numpoints);
size_t size_dt = sizeof(dt);
size_t size_to = sizeof(to);
size_t size_tf = sizeof(tf);
size_t h_pitch = numpoints*sizeof(double);
size_t d_pitch;

//Calculate the "pitch" of the 2D array.  The pitch is basically the length of a 2D array's row.  IT's larger
//than the actual row full of data due to hadware issues.  We thusly will use the pitch instead of the data
//size to traverse the array.
gpuErrchk(cudaMallocPitch( (void**)&d_concStorage, &d_pitch, numpoints * sizeof(double), maxlength));

//Allocate memory on the GPU for all the arrrays we're going to use in the integrator.
gpuErrchk(cudaMalloc((void**)&d_concs, size_concs));
gpuErrchk(cudaMalloc((void**)&d_temp1, size_temp1));
gpuErrchk(cudaMalloc((void**)&d_temp2, size_temp1));
gpuErrchk(cudaMalloc((void**)&d_tempsum, size_tempsum));
gpuErrchk(cudaMalloc((void**)&d_k1s, size_ks));
gpuErrchk(cudaMalloc((void**)&d_k2s, size_ks));
gpuErrchk(cudaMalloc((void**)&d_k3s, size_ks));
gpuErrchk(cudaMalloc((void**)&d_k4s, size_ks));
gpuErrchk(cudaMalloc((void**)&d_k5s, size_ks));
gpuErrchk(cudaMalloc((void**)&d_k6s, size_ks));
gpuErrchk(cudaMalloc((void**)&d_maxlength, size_maxlength));
gpuErrchk(cudaMalloc((void**)&d_numpoints, size_numpoints));
gpuErrchk(cudaMalloc((void**)&d_dt, size_dt));
gpuErrchk(cudaMalloc((void**)&d_to, size_to));
gpuErrchk(cudaMalloc((void**)&d_tf, size_tf));

//Copy all initial values of arrays to GPU.
gpuErrchk(cudaMemcpy2D(d_concStorage, d_pitch, concStorage, h_pitch, numpoints*sizeof(double), maxlength, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_concs, &concs, size_concs, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_temp1, &temp1, size_temp1, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_temp2, &temp2, size_temp2, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_tempsum, &tempsum, size_tempsum, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_k1s, &k1s, size_ks, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_k2s, &k2s, size_ks, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_k3s, &k3s, size_ks, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_k4s, &k4s, size_ks, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_k5s, &k5s, size_ks, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_k6s, &k6s, size_ks, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_maxlength, &maxlength, size_maxlength, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_numpoints, &numpoints, size_numpoints, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_dt, &dt, size_dt, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_to, &to, size_to, cudaMemcpyHostToDevice));
gpuErrchk(cudaMemcpy(d_tf, &tf, size_tf, cudaMemcpyHostToDevice));

//Run the integrator.
rkf5<<<1,1>>>(d_concs, d_concStorage, d_temp1, d_temp2, d_tempsum, d_k1s, d_k2s, d_k3s, d_k4s, d_k5s, d_k6s, d_maxlength, d_numpoints, d_pitch, d_dt, d_to, d_tf);
gpuErrchk( cudaPeekAtLastError() );

//Copy concentrations from GPU to Host.  Almost defunct now that transferring the 2D array works.
cudaMemcpy(concs, d_concs, size_concs, cudaMemcpyDeviceToHost);
//Copy 2D array of concentrations vs. time from GPU to Host.
gpuErrchk( cudaMemcpy2D(concStorage, h_pitch, d_concStorage, d_pitch, numpoints*sizeof(double), maxlength, cudaMemcpyDeviceToHost) );

//Print concentrations after the integrator kernel runs.  Used to test that data was transferring to and from GPU correctly.
std::cout << "\n";
for (int i=0; i<maxlength; i++)
{
std::cout<<concs[i];
std::cout<<" ";
}

double a;
double b;
double c;
for(int i = 0; i< 10; i++)
{
a[i]=0;
b[i]=0;
c[i]=0;
}

//Print out the concStorage array after the kernel runs.  Used to test that the 2D array transferred correctly from host to GPU and back.
std::cout << "\n\n";
std::cout << "Calculated Array";
std::cout << "\n\n";
for (int i=0; i<maxlength; i++)
{
for(int j=0; j<numpoints; j++)
{
if (j%(numpoints/10)==0)
{
a[j/(numpoints/10)]=concStorage[i][j];
std::cout<<concStorage[i][j];
std::cout<<"   ";
}
}
std::cout << "\n";
}
std::cout << "\n";
std::cout << "Exponential";
std::cout << "\n\n";
for (int i=0; i<10; i++)
{
b[i]=exp(-i*(tf-to)/10);
std::cout<<exp(-i*(tf-to)/10);
std::cout<<"   ";
}
std::cout << "\n\n";
std::cout << "Error Array";
std::cout << "\n\n";
for (int i=0; i<10; i++)
{
c[i]=a[i]-b[i];
std::cout<<c[i];
std::cout<<"   ";
}
std::cout << "\n\n";

return 0;
}
//Main kernel.  This is mean to be run as a master thread that calls all the other functions and thusly "runs" the integrator.
__global__ void rkf5(double* concs, double* concStorage, double* temp1, double* temp2, double* tempsum, double* k1s, double* k2s, double* k3s, double* k4s, double* k5s, double* k6s, int* maxlength, int* numpoints, size_t pitch, double* dt, double* to, double* tf)
{
/*
axy variables represent the coefficients in the Butcher's tableau where x represents the order of the step and the y value corresponds to the ky value
the coefficient gets multiplied by.  Have to cast them all as doubles, or the ratios evaluate as integers.
e.g. a21 -> a21 * k1
e.g. a31 -> a31 * k1 + a32 * k2
*/
double a21 = static_cast<double>(.25);

double a31 = static_cast<double>(3)/static_cast<double>(32);
double a32 = static_cast<double>(9)/static_cast<double>(32);

double a41 = static_cast<double>(1932)/static_cast<double>(2197);
double a42 = static_cast<double>(-7200)/static_cast<double>(2197);
double a43 = static_cast<double>(7296)/static_cast<double>(2197);

double a51 = static_cast<double>(439)/static_cast<double>(216);
double a52 = static_cast<double>(-8);
double a53 = static_cast<double>(3680)/static_cast<double>(513);
double a54 = static_cast<double>(-845)/static_cast<double>(4104);

double a61 = static_cast<double>(-8)/static_cast<double>(27);
double a62 = static_cast<double>(2);
double a63 = static_cast<double>(-3544)/static_cast<double>(2565);
double a64 = static_cast<double>(1859)/static_cast<double>(4104);
double a65 = static_cast<double>(-11)/static_cast<double>(40);

//for loop that integrates over the specified number of points. Actually, might have to make it a do-while loop for adaptive step sizes
for(int k = 0; k < *numpoints; k++)
{
if (k!=0)
{
arrCopy<<< 1, *maxlength >>>(concs, tempsum);
}
arrInit<<< 1, *maxlength >>>(tempsum, 0);
arrInit<<< 1, *maxlength >>>(temp1, 0);
arrInit<<< 1, *maxlength >>>(temp2, 0);

calcK<<< 1, *maxlength >>>(concs, k1s, dt);         //k1 = dt * flux (concs)
cudaDeviceSynchronize(); //Sync here because kernel continues onto next line before k1 finished

setup(a21, temp1, tempsum, k1s, concs, maxlength);  //tempsum = a21*k1
arrAdd<<< 1, *maxlength >>>(concs, tempsum, tempsum);   //tempsum = concs + a21*k1

calcK<<< 1, *maxlength >>>(tempsum, k2s, dt);       //k2 = dt * flux (concs + a21*k1)

arrInit<<< 1, *maxlength >>>(tempsum, 0);
setup(a31, temp1, tempsum, k1s, concs, maxlength);  //temp1sum = a31*k1
setup(a32, temp1, tempsum, k2s, concs, maxlength);  //tempsum = a31*k1 + a32*k2
arrAdd<<< 1, *maxlength >>>(concs, tempsum, tempsum);   //tempsum = concs + a31*k1 + a32*k2

calcK<<< 1, *maxlength >>>(tempsum, k3s, dt);           //k3 = dt * flux (concs + a31*k1 + a32*k2)

arrInit<<< 1, *maxlength >>>(tempsum, 0);
setup(a41, temp1, tempsum, k1s, concs, maxlength);  //tempsum = a41*k1
setup(a42, temp1, tempsum, k2s, concs, maxlength);  //tempsum = a41*k1 + a42*k2
setup(a43, temp1, tempsum, k3s, concs, maxlength);  //tempsum = a41*k1 + a42*k2 + a43*k3
arrAdd<<< 1, *maxlength >>>(concs, tempsum, tempsum);   //tempsum = concs + a41*k1 + a42*k2 + a43*k3

calcK<<< 1, *maxlength >>>(tempsum, k4s, dt);           //k4 = dt * flux (concs + a41*k1 + a42*k2 + a43*k3)

arrInit<<< 1, *maxlength >>>(tempsum, 0);
setup(a51, temp1, tempsum, k1s, concs, maxlength);  //tempsum = a51*k1
setup(a52, temp1, tempsum, k2s, concs, maxlength);  //tempsum = a51*k1 + a52*k2
setup(a53, temp1, tempsum, k3s, concs, maxlength);  //tempsum = a51*k1 + a52*k2 + a53*k3
setup(a54, temp1, tempsum, k4s, concs, maxlength);  //tempsum = a51*k1 + a52*k2 + a53*k3 + a54*k4
arrAdd<<< 1, *maxlength >>>(concs, tempsum, tempsum);   //tempsum = concs + a51*k1 + a52*k2 + a53*k3 + a54*k4

calcK<<< 1, *maxlength >>>(tempsum, k5s, dt);           //k5 = dt * flux (concs + a51*k1 + a52*k2 + a53*k3 + a54*k4)

arrInit<<< 1, *maxlength >>>(tempsum, 0);
setup(a61, temp1, tempsum, k1s, concs, maxlength);  //tempsum = a61*k1
setup(a62, temp1, tempsum, k2s, concs, maxlength);  //tempsum = a61*k1 + a62*k2
setup(a63, temp1, tempsum, k3s, concs, maxlength);  //tempsum = a61*k1 + a62*k2 + a63*k3
setup(a64, temp1, tempsum, k4s, concs, maxlength);  //tempsum = a61*k1 + a62*k2 + a63*k3 + a64*k4
setup(a65, temp1, tempsum, **k4s**, concs, maxlength);  //tempsum = a61*k1 + a62*k2 + a63*k3 + a64*k4 + a65*k5
arrAdd<<< 1, *maxlength >>>(concs, tempsum, tempsum);   //tempsum = concs + a61*k1 + a62*k2 + a63*k3 + a64*k4 + a65*k5

calcK<<< 1, *maxlength >>>(tempsum, k6s, dt);           //k6 = dt * flux (concs + a61*k1 + a62*k2 + a63*k3 + a64*k4 + a65*k5)

//At this point, temp1 and tempsum are maxlength dimension arrays that are able to be used for other things.

//Calculate acceptable step size before storing the concentrations.
calcStepSize(temp1, temp2, tempsum, concs, k1s, k2s, k3s, k4s, k5s, k6s, dt, maxlength);    //temp1 = 4th Order guess, tempsum = 5th Order guess

//Store the initial conditions in the first column of the storage array.
if (k==0)
{
storeConcs<<< 1, *maxlength >>>(concStorage, pitch, concs, k);  //Store this step's concentrations in 2D array
}
//Store future concentration in next column of storage array.
storeConcs<<< 1, *maxlength >>>(concStorage, pitch, tempsum, k+1);  //Store this step's concentrations in 2D array
}
}
//calcStepSize will take in an error tolerance, the current concentrations and the k values and calculate the resulting step size according to the following equation
//e[n+1]=y4[n+1] - y5[n+1]
__device__ void calcStepSize(double* temp1, double*temp2, double* tempsum, double* concs, double* k1s, double* k2s, double* k3s, double* k4s, double* k5s, double* k6s, double* dt, int* maxlength)
{
//do
//{
takeFourthOrderStep<<< 1, *maxlength >>>(temp1, concs, k1s, k2s, k3s, k4s, k5s);            //Store 4th order guess in temp1
takeFifthOrderStep<<< 1, *maxlength >>>(tempsum, concs, k1s, k2s, k3s, k4s, k5s, k6s);  //Store 5th order guess in tempsum
arrSub<<< 1, *maxlength >>>(temp1, tempsum, temp2)
arrMin<<< 1, *maxlength >>>
//arrMult
//}
//while
}
//takeFourthOrderStep is going to overwrite the old temp1 array with the new array of concentrations that result from a 4th order step.  This kernel is meant to be launched
// with as many threads as there are discrete concentrations to be tracked.
__global__ void takeFourthOrderStep(double* y4, double* concs, double* k1s, double* k2s, double* k3s, double* k4s, double* k5s)
{
double b41 = static_cast<double>(25)/static_cast<double>(216);
double b42 = static_cast<double>(0);
double b43 = static_cast<double>(1408)/static_cast<double>(2565);
double b44 = static_cast<double>(2197)/static_cast<double>(4104);
double b45 = static_cast<double>(-1)/static_cast<double>(5);
int idx = blockIdx.x * blockDim.x + threadIdx.x;
y4[idx] = concs[idx] + b41 * k1s[idx] + b42 * k2s[idx] + b43 * k3s[idx] + b44 * k4s[idx] + b45 * k5s[idx];
}
//takeFifthOrderStep is going to overwrite the old array of concentrations with the new array of concentrations.  As of now, this will be the 5th order step.  Another function can be d
//defined that will take a fourth order step if that is interesting for any reason.  This kernel is meant to be launched with as many threads as there are discrete concentrations
//to be tracked.
//Store b values in register? Constants?
__global__ void takeFifthOrderStep(double* y5, double* concs, double* k1s, double* k2s, double* k3s, double* k4s, double* k5s, double* k6s)
{
double b51 = static_cast<double>(16)/static_cast<double>(135);
double b52 = static_cast<double>(0);
double b53 = static_cast<double>(6656)/static_cast<double>(12825);
double b54 = static_cast<double>(28561)/static_cast<double>(56430);
double b55 = static_cast<double>(-9)/static_cast<double>(50);
double b56 = static_cast<double>(2)/static_cast<double>(55);
int idx = blockIdx.x * blockDim.x + threadIdx.x;
y5[idx] = concs[idx] + b51 * k1s[idx] + b52 * k2s[idx] + b53 * k3s[idx] + b54 * k4s[idx] + b55 * k5s[idx] + b56 * k6s[idx];
}
//storeConcs takes the current array of concentrations and stores it in the cId'th column of the 2D concStorage array
//pitch = memory size of a row
__global__ void storeConcs(double* cS, size_t pitch, double* concs, int cId)
{
//cS is basically the memory address of the first element of the flattened (1D) 2D array.
double* row = (double*)((char*)cS + tIdx * pitch);
row[cId] = concs[tIdx];
}
//Perhaps I can optimize by using shared memory to hold conc values.
__global__ void calcK(double* concs, double* ks, double* dt)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
ks[idx]=(*dt)*flux(idx, concs);
}
//Adds two arrays (a + b) element by element and stores the result in array c.
__global__ void arrAdd(double* a, double* b, double* c)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
c[idx]=a[idx]+b[idx];
}
//Subtracts two arrays (a - b) element by element and stores the result in array c.
__global__ void arrSub(double* a, double* b, double* c)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
c[idx]=a[idx]-b[idx];
}
//Multiplies two arrays (a * b) element by element and stores the result in array c.
__global__ void arrMult(double* a, double* b, double* c)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
c[idx]=a[idx]*b[idx];
}
//Will find the min of errors array.
__global__ void arrMin(double* errors)
{
//extern _shared_ double;
}
//Initializes an array a to double value b.
__global__ void arrInit(double* a, double b)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
a[idx]=b;
}
//Copies array b onto array a.
__global__ void arrCopy(double* a, double* b)
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
a[idx]=b[idx];
}
//Placeholder function for the flux calculation.  It will take the size of the oligomer and current concentrations as inputs.
__device__ double flux(int r, double *concs)
{
return -concs[r];
}
//This function multiplies a tableau value by the corresponding k array and adds the result to tempsum.  Used to
//e.g. setup(a21, temp1, tempsum, k1s, concs, maxlength) => tempsum = a21 * k1
__device__ void setup(double tableauValue, double *temp1, double *tempsum, double *ks, double *concs, int *maxlength)
{
//Sets tempsum to tabVal * k
arrInit<<< 1, *maxlength >>>(temp1, tableauValue);      //Set [temp1] to tableau value, temp1 = a
arrMult<<< 1, *maxlength >>>(ks, temp1, temp1);         //Multiply tableau value by appropriate [k], temp1 = a*k
arrAdd<<< 1, *maxlength >>>(tempsum, temp1, tempsum);   //Move tabVal*k to [tempsum], tempsum = tempsum+temp1
//temp1 = tableauValue * kArray
//tempsum = current sum (tableauValue * kArray)
}


FINAL EDIT I found my error. It was a typo calculating the K6's. I put asterisks around the offending "k4s" term which should be a "k5s". Anyone who wants a CUDA RKF45 integrator is welcome to use it as it works now and is of the correct order.

Number of Points    h       Global Error    Calculated Order
10          0.1         3.59E-09
20          0.05        1.08E-10    5.06E+00
40          0.025       3.31E-12    5.03E+00
80          0.0125      1.03E-13    5.01E+00
160         0.00625     2.20E-15    5.55E+00
320         0.003125    1.20E-15    8.74E-01
640         0.0015625   1.40E-15    -2.22E-01
1280        0.00078125  2.00E-16    2.81E+00

• Is there any reason you're not using a library that implements this method? – Geoff Oxberry Nov 8 '13 at 19:28
• If you're new to C++, CUDA and wirting ODE solvers this task seems quite challenging. If you have access to MatLab with the Parallel Programming Tool box and an NVidia Graphics card, it might be worth wile to prototype this code in MatLab. It has support for CUDA GPUs (http://www.mathworks.com/discovery/matlab-gpu.html) (and OpenMP and MPI for that matter) these days. Once the code is running in MatLab, you could port it over to C++ or use MeX. Just a thought. – seb Nov 9 '13 at 6:06
• You could use the odeint library (part of Boost) with the Thrust interface. If you're going to roll your own code, you should test your implementation with the method of manufactured solutions. – Geoff Oxberry Nov 9 '13 at 6:14
• @Karsten: It might be worth redoing your analysis assuming that the global error is of the form $E(h) = a * h^{n}$. Then $\log{E} = \log{a} + n \log{h}$, so a log-log plot (or a linear regression on the log-transformed data) will tell you what the error is. A plot would also be more informative to readers of your post; to make sense of your data, I'd have to either plot it myself to see the order of your method, or I'd have to do the regression myself. – Geoff Oxberry Nov 10 '13 at 3:13
• @KarstenChu point taken. I prefer the Python/C combo myself; I just like to keep an open mind. To see what you can do with MatLab these days check out this link from the PDC Cluster in Stockholm, Sweden. They have slides and a lot of example matlab scripts. They include openmp, MPI, cuda. – seb Nov 10 '13 at 9:58

I am assuming that you are setting the error tolerance at 1e-12. You are correct that when an adaptive scheme accepts the current step size, it assumes the 5th order scheme was, for all intents and purposes, the "correct" answer. However this is only when it accepts the current step. If the difference between the 4th and 5th order steps are too large, it doesn't accept either as good enough. For testing a single algorithm like this you need to run it with multiple time steps and see the order of convergence to verify it.

Now on the other hand, a step of 0.001 with a 5th order scheme on your test problem (I am assuming was 1D) should be more than sufficient to drive the error to nearly machine epsilon. My guess is that you have a typo somewhere that reduces your integrator to 1st order. Like I said in paragraph 1, run it with multiple time steps and verify its convergence order. This will tell you whether the error you see is reasonable, or if you have an error somewhere in your integrator.

EDIT

To calculate the convergence order of an ode solver, run it on a problem with a known solution for two different time steps $h$ and $\frac{h}{2}$. Your global error should then be:

$\frac{E(h)}{E(\frac{h}{2})} = 2^n$

Note that the 2 comes from the ratio between time steps. You should get $n=5$ for your 5th order scheme.

• As of right now, I don't have an adaptive step algorithm programmed, so that is probably contributing to the error. For the fixed sizes, however, I did take a fun little jaunt into logarithms to calculate the approximate order of my solver as is and I'm getting ~2.1. Using the difference between the 5th order guess and the accepted value as my "global error," I did O = ln(GlobalError/NumberOfPoints) / ln(StepSize) for step sizes of .001, .01 and .1. Does that suggest that I have an error in my code? – Hair of Slytherin Nov 9 '13 at 3:07
• See my edit for how to calculate the order. – Godric Seer Nov 9 '13 at 15:38
• I added a chart where I took several runs of my data and associated errors and plugged them into the equation O = LN( E(h1)/E(h2) ) / LN (h1/h2) for arbitrary stepsizes. That's what I put in the chart as "Your Order" – Hair of Slytherin Nov 10 '13 at 17:01
• Look at the last two data points in your table. You reduced the step size by a factor of 10 and your error dropped by a factor of about 100. Your method is definitely 2nd order. Without seeing your code I don't know what to suggest other than rechecking your coefficients. – Godric Seer Nov 10 '13 at 20:28
• All of your coefficients are correct. The one thing I didn't see was where you were modifying dt for calculating the different K's (i.e. f(t+adt, x+ aK)). Since your equations (test and real) appear to be autonomous, that shouldn't matter though. – Godric Seer Nov 11 '13 at 14:20