# How to modify adaptive step size Runge-Kutta cash karp algorithm for higher accuracy?

I'm following a paper on making your own cosmic microwave background perturbation solver code by Peter Callin https://arxiv.org/pdf/astro-ph/0606683.pdf In the programming techniques section V, the author says to use Cash Karp adaptive step size Runge Kutta method for numerical integration of the perturbation equations, and he references the algorithm from Numerical Recipies in C book with tolerance eps = 1E-11.

I'm trying to do that and the code below is taken from that book. The problem is I'm not able to solve my system with high accuracy or low eps. I always get errors for eps < 1E-4. I've already changed the types to double everywhere.


///// My system//////////
double z(double x){
return sqrt(omega_dark_energy*pow(10, x*4.0) + omega_m*pow(10, x) + omega_radiation);
}

int index_phi = 1, index_theta_0 = 2, index_theta_1 = 3, index_delta_dm = 4, index_velocity_dm = 5;

void derivs(double x, double y[], double dydx[]){
/* 1 = potential
2 = theta-0
3 = theta-1
4 = delta-dm
5 = velocity-dm */
y = vector(1,nvar);
dydx = vector(1, nvar);
dydx[index_phi] = log(10)*((pow(10,x)*omega_dm*y[index_delta_dm] + 4.0*omega_radiation*y[index_theta_0])/(2.0*z(x)*z(x)) - y[index_phi]*(pow(10,x*2)*pow(wavenumber,2)*cspeed*cspeed/(3.0*pow(Hubble_const_mpc,2)*z(x)*z(x)) + 1));
dydx[index_theta_0] = log(10)*(-wavenumber*y[index_theta_1]*pow(10,x)*cspeed/(Hubble_const_mpc*z(x)) - (pow(10,x)*omega_dm*y[index_delta_dm] + 4*omega_radiation*y[index_theta_0])/(2*z(x)*z(x)) + y[index_phi]*(pow(10,x*2)*pow(wavenumber,2)*cspeed*cspeed/(3*pow(Hubble_const_mpc,2)*z(x)*z(x)) + 1));
dydx[index_theta_1] = log(10)*((y[index_theta_0] - y[index_phi])*wavenumber*pow(10,x)*cspeed/(3.0*Hubble_const_mpc*z(x)));
dydx[index_delta_dm] = log(10)*(-wavenumber*pow(10,x)*y[index_velocity_dm]*cspeed/(Hubble_const_mpc*z(x)) - 3.0*(pow(10,x)*omega_dm*y[index_delta_dm] + 4*omega_radiation*y[index_theta_0])/(2*z(x)*z(x)) + 3.0*y[index_phi]*(pow(10,x*2)*pow(wavenumber,2)*cspeed*cspeed/(3*pow(Hubble_const_mpc,2)*z(x)*z(x)) + 1));
dydx[index_velocity_dm] = log(10)*(-wavenumber*y[index_phi]*pow(10,x)*cspeed/(Hubble_const_mpc*z(x)) - y[index_velocity_dm]);
}
//////////////////////////

//////Cash Karp algorithm taken from the referenced book
void rkck(double y[], double dydx[], int n, double x, double h, double yout[],
double yerr[], void (*derivs)(double, double [], double [])){
/*Given values for n variables y[1..n] and their derivatives dydx[1..n] known at x, use
the fifth-order Cash-Karp Runge-Kutta method to advance the solution over an interval h
and return the incremented variables as yout[1..n]. Also return an estimate of the local
truncation error in yout using the embedded fourth-order method. The user supplies the routine
derivs(x,y,dydx), which returns derivatives dydx at x.*/
int i;
static double a2=0.2,a3=0.3,a4=0.6,a5=1.0,a6=0.875,b21=0.2,
b31=3.0/40.0,b32=9.0/40.0,b41=0.3,b42 = -0.9,b43=1.2,
b51 = -11.0/54.0, b52=2.5,b53 = -70.0/27.0,b54=35.0/27.0,
b61=1631.0/55296.0,b62=175.0/512.0,b63=575.0/13824.0,
b64=44275.0/110592.0,b65=253.0/4096.0,c1=37.0/378.0,
c3=250.0/621.0,c4=125.0/594.0,c6=512.0/1771.0,
dc5 = -277.00/14336.0;
double dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0,
dc4=c4-13525.0/55296.0,dc6=c6-0.25;
double *ak2,*ak3,*ak4,*ak5,*ak6,*ytemp;
ak2=vector(1,n);
ak3=vector(1,n);
ak4=vector(1,n);
ak5=vector(1,n);
ak6=vector(1,n);
ytemp=vector(1,n);
for (i=1;i<=n;i++)  //First step.
ytemp[i]=y[i]+b21*h*dydx[i];
(*derivs)(x+a2*h,ytemp,ak2); //Second step
for (i=1;i<=n;i++)
ytemp[i]=y[i]+h*(b31*dydx[i]+b32*ak2[i]);
(*derivs)(x+a3*h,ytemp,ak3);  //Third step.
for (i=1;i<=n;i++)
ytemp[i]=y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);
(*derivs)(x+a4*h,ytemp,ak4);   //Fourth step
for (i=1;i<=n;i++)
ytemp[i]=y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);
(*derivs)(x+a5*h,ytemp,ak5);  //Fifth step.
for (i=1;i<=n;i++)
ytemp[i]=y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);
(*derivs)(x+a6*h,ytemp,ak6);  // Sixth step.
for (i=1;i<=n;i++) //Accumulate increments with proper weights.
yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);
for (i=1;i<=n;i++)
yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
//Estimate error as difference between fourth and fifth order methods.
free_vector(ytemp,1,n);
free_vector(ak6,1,n);
free_vector(ak5,1,n);
free_vector(ak4,1,n);
free_vector(ak3,1,n);
free_vector(ak2,1,n);
}

void rkqs(double y[], double dydx[], int n, double *x, double htry, double eps,
double yscal[], double *hdid, double *hnext,
void (*derivs)(double, double [], double [])){
/*Fifth-order Runge-Kutta step with monitoring of local truncation error to ensure accuracyand
adjust stepsize. Input are the dependent variable vector y[1..n] and its derivative dydx[1..n]
at the starting value of the independent variable x. Also input are the stepsize to be attempted
htry, the required accuracy eps, and the vector yscal[1..n] against which the error is
scaled. On output, y and x are replaced bytheir new values, hdid is the stepsize that was
actually accomplished, and hnext is the estimated next stepsize. derivs is the user-supplied
routine that computes the right-hand side derivatives.*/
void rkck(double y[], double dydx[], int n, double x, double h,
double yout[], double yerr[], void (*derivs)(double, double [], double []));
int i;
double errmax,h,htemp,xnew,*yerr,*ytemp;
yerr=vector(1,n);
ytemp=vector(1,n);
h=htry;  //Set stepsize to the initial trial value.
for (;;) {
rkck(y,dydx,n,*x,h,ytemp,yerr,derivs); //Take a step
errmax=0.0;  //Evaluate accuracy.
for (i=1;i<=n;i++) errmax=FMAX(errmax,fabs(yerr[i]/yscal[i]));
errmax /= eps; //Scale relative to required tolerance
if (errmax <= 1.0) break; //Step succeeded. Compute size of next step.
htemp=SAFETY*h*pow(errmax,PSHRNK);
// Truncation error too large, reduce stepsize.
h=(h >= 0.0 ? FMAX(htemp,0.1*h) : FMIN(htemp,0.1*h));
// No more than a factor of 10.
xnew=(*x)+h;
// if (xnew == *x) nrerror("stepsize underflow in rkqs");
}
if (errmax > ERRCON) *hnext=SAFETY*h*pow(errmax,PGROW);
else *hnext=10.0*h; //No more than a factor of 5 increase. //Changed to 10
*x += (*hdid=h);
for (i=1;i<=n;i++) y[i]=ytemp[i];
free_vector(ytemp,1,n);
free_vector(yerr,1,n);
}

/*User storage for intermediate results. Preset kmax and dxsav in the calling program. If kmax =
0 results are stored at approximate intervals dxsav in the arrays xp[1..kount], yp[1..nvar]
[1..kount], where kount is output by odeint. Defining declarations for these variables, with
memoryallocations xp[1..kmax] and yp[1..nvar][1..kmax] for the arrays, should be in
the calling program.*/
void odeint(double ystart[], int nvar, double x1, double x2, double eps, double h1,
double hmin, int nok, int nbad,
void (*derivs)(double, double [], double []),
void (*rkqs)(double [], double [], int, double *, double, double, double [],
double *, double *, void (*)(double, double [], double []))){
/*Runge-Kutta driver with adaptive stepsize control. Integrate starting values ystart[1..nvar]
from x1 to x2 with accuracy eps, storing intermediate results in global variables. h1 should
be set as a guessed first stepsize, hmin as the minimum allowed stepsize (can be zero). On
output nok and nbad are the number of good and bad (but retried and fixed) steps taken, and
ystart is replaced byvalues at the end of the integration interval. derivs is the user-supplied
routine for calculating the right-hand side derivative, while rkqs is the name of the stepper
routine to be used.*/
int nstp,i;
double xsav,x,hnext,hdid,h;
double *yscal,*y,*dydx;
yscal=vector(1,nvar);
y=vector(1,nvar);
dydx=vector(1,nvar);
// INITIAL CONDIIONS
x=x1;
h=SIGN(h1,x2-x1);
nok = nbad = kount = 0;
for (i=1;i<=nvar;i++) y[i]=ystart[i];
if (kmax > 0) xsav=x-dxsav*2.0; //Assures storage of first step.
for (nstp=1;nstp<=MAXSTP;nstp++) { //Take at most MAXSTP steps.
(*derivs)(x,y,dydx);
for (i=1;i<=nvar;i++)
//Scaling used to monitor accuracy. This general-purpose choice can be modified if need be.
yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
if (kmax > 0 && kount < kmax-1 && fabs(x-xsav) > fabs(dxsav)) {
xp[++kount]=x; //Store intermediate results.
for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
xsav=x;
}
if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x; //If stepsize can overshoot, decrease.
(*rkqs)(y,dydx,nvar,&x,h,eps,yscal,&hdid,&hnext,derivs);
if (hdid == h) ++nok; else ++nbad;
if ((x-x2)*(x2-x1) >= 0.0) { //Are we done?
for (i=1;i<=nvar;i++) ystart[i]=y[i];
if (kmax) {
xp[++kount]=x; //Save final step.
for (i=1;i<=nvar;i++) yp[i][kount]=y[i];
}
free_vector(dydx,1,nvar);
free_vector(y,1,nvar);
free_vector(yscal,1,nvar);
return; //Normal exit.
}
if (fabs(hnext) <= hmin) nrerror("Step size too small in odeint");
h=hnext;
}
nrerror("Too many steps in routine odeint");
}



I always end up getting some sort of error like step size underflow or max number of steps when i try to aim for high accuracy. I cannot seem to understand how did the author of the paper modify this algorithm to solve such a complicated system.

• A good starting point to debug problems is to carefully read the error message and try to understand (i) what the error message actually says, (ii) in a debugger get to that point and understand why the error appears. Right now you seem to not care very much what the error says "I always end up getting some sort of error like...", but it's quite useful to spend time understand the error when trying to understand its cause. Apr 3 at 22:14
• In derivs you define a new vector y filled with zeros or garbage. This supersedes the argument y so that the derivatives are constructed for the wrong state. // When rolling your own solver it is always a good practice to test it on some equations where the solution is known, such as the oscillator, or some MMS construct. Apr 4 at 7:38
• @LutzLehmann Thank you so so much, it's working for a part of my system now, and I'll try to solve the full system now. Thank you very much!!!! Edit: How do I mark the question as answered? Apr 5 at 13:51
• I have made my comment into an answer. With some delay, you should be able to mark this as the (best) answer. Apr 5 at 14:07

At the state of the code in OPost, there is no hope to get any solution, as the passing of arguments in the ODE function derivs is interrupted by creating a local variable y that supersedes the passed argument y.