#include <math.h>
#include "randomize.h"
//"thread safe" version
namespace libpetey {
//derivative functions are compatible with the GSL
//(return value is ignored)
template <class ind_type, class dep_type>
dep_type **rk_dumb(ind_type t0, dep_type *xinit, long nx, ind_type dt, long nt, void *param,
int (* derivs) (ind_type, dep_type *, dep_type *, void *) ) {
dep_type **xvec;
ind_type t, tph;
dep_type dxall[nx], dx[nx];
dep_type xint[nx];
long i, j;
xvec=new dep_type * [nt+1];
xvec[0]=new dep_type [nx];
for (j=0;j<nx;j++) xvec[0][j]=xinit[j];
t=t0;
for (i=1;i<=nt;i++) {
//allocate space for the next vector of dep. values:
xvec[i]=new dep_type[nx];
//take the first set of derivatives:
(* derivs) (t, xvec[i-1], dx, param);
for (j=0;j<nx;j++) {
xint[j]=xvec[i-1][j]+dx[j]*dt/2; //calculate intermediate dep. values
dxall[j]=dx[j]/6; //start calculating the total deriv
}
//take the second set of derivatives:
tph=dt/2.;
tph=t+tph;
(* derivs) (tph, xint, dx, param);
for (j=0;j<nx;j++) {
xint[j]=xvec[i-1][j]+dx[j]*dt/2;
dxall[j]=dxall[j]+dx[j]/3;
}
//take the third set of derivatives:
(* derivs) (tph, xint, dx, param);
for (j=0;j<nx;j++) {
xint[j]=xvec[i-1][j]+dx[j]*dt;
dxall[j]=dxall[j]+dx[j]/3;
}
//take the fourth set of derivatives:
t=t+dt;
(* derivs) (t, xint, dx, param);
for (j=0;j<nx;j++) {
xvec[i][j]=xvec[i-1][j]+(dxall[j]+dx[j]/6)*dt; //calculate next vector
}
}
return xvec;
}
template <class ind_type, class dep_type>
void rk_dumb(ind_type t0, dep_type *xinit, long nx, ind_type dt, long nt,
dep_type ** xvec, void *param,
int (* derivs) (ind_type, dep_type *, dep_type *, void *) ) {
ind_type t, tph;
dep_type dxall[nx], dx[nx];
dep_type xint[nx];
long i, j;
for (j=0;j<nx;j++) xvec[0][j]=xinit[j];
t=t0;
for (i=1;i<=nt;i++) {
//take the first set of derivatives:
(* derivs) (t, xvec[i-1], dx, param);
for (j=0;j<nx;j++) {
xint[j]=xvec[i-1][j]+dx[j]*dt/2; //calculate intermediate dep. values
dxall[j]=dx[j]/6; //start calculating the total deriv
}
//take the second set of derivatives:
tph=dt/2.;
tph=t+tph;
(* derivs) (tph, xint, dx, param);
for (j=0;j<nx;j++) {
xint[j]=xvec[i-1][j]+dx[j]*dt/2;
dxall[j]=dxall[j]+dx[j]/3;
}
//take the third set of derivatives:
(* derivs) (tph, xint, dx, param);
for (j=0;j<nx;j++) {
xint[j]=xvec[i-1][j]+dx[j]*dt;
dxall[j]=dxall[j]+dx[j]/3;
}
//take the fourth set of derivatives:
t=t+dt;
(* derivs) (t, xint, dx, param);
for (j=0;j<nx;j++) {
xvec[i][j]=xvec[i-1][j]+(dxall[j]+dx[j]/6)*dt; //calculate next vector
}
}
}
//generalized power function:
template <class real>
int test_rk_pfunc(real t, real *x, real *f, void *param) {
int err=0;
void **p2=(void **) param;
int nterm=*(int *) p2[0]; //number of terms
real *c=(real *) p2[1]; //coefficients
real t_n; //t^n
real fact; //n!
t_n=1;
fact=1;
*f=c[0];
for (int i=1; i<nterm; i++) {
t_n*=t;
fact*=i;
*f+=c[i]*t_n/fact;
}
return err;
}
//indefinite integral of generalized power function:
template <class real>
real test_rk_ipfunc(real t, void *param) {
void **p2=(void **) param;
int nterm=*(int *) p2[0]; //number of terms
real *c=(real *) p2[1]; //coefficients
real t_n; //t^n
real fact; //n!
real f;
t_n=1;
fact=1;
f=0;
for (int i=1; i<=nterm; i++) {
t_n*=t;
fact*=i;
f+=c[i-1]*t_n/fact;
}
return f;
}
template <class real>
int test_rk_ts(int nterm, //number of terms
real t1, //integration limits
real t2,
real h) { //step size
void *param[2];
real c[nterm];
real x0=0;
long nt=(t2-t1)/h;
real **x;
real x1, x2;
real xanal;
real k;
real errp, erra; //predicted and actual error
int exit_code;
x=new real *[nt+2];
x[0]=new real[nt+2];
for (int i=0; i<nt+2; i++) x[i]=x[0]+i;
//random coefficients:
for (int i=0; i<nterm; i++) c[i]=ranu();
//initialize parameters:
param[0]=&nterm;
param[1]=c;
//given step size:
rk_dumb(t1, &x0, 1, h, nt, x, param, &test_rk_pfunc<real>);
//initial condition gets copied to itself (***):
rk_dumb(h*nt, x[nt], 1, t2-h*nt, 1, x+nt, param, &test_rk_pfunc<real>);
x1=x[nt+1][0];
//double step size:
nt=(t2-t1)/h/2;
rk_dumb(t1, &x0, 1, 2*h, nt, x, param, &test_rk_pfunc<real>);
//initial condition gets copied to itself (***):
rk_dumb(2*nt*h, x[nt], 1, t2-2*h*nt, 1, x+nt, param, &test_rk_pfunc<real>);
x2=x[nt+1][0];
//analytic result:
xanal=test_rk_ipfunc(t2, param)-test_rk_ipfunc(t1, param);
//predicted error:
k=(x1-x2)/pow(h, 5)/15;
errp=(x1-x2)/15;
erra=xanal-x1;
printf("R-K integrator:\n");
printf("numerical = %g + %g\n", x1, (x1-x2)/15);
printf("analytic = %g (dif=%g)\n", xanal, xanal-x1);
//if error is the right order-of magnitude, then we say it passed:
if (fabs(2*(errp-erra))/(errp+erra)<1) exit_code=0; else exit_code=1;
delete [] x[0];
delete [] x;
return exit_code;
}
template int test_rk_ts<float>(int, float, float, float);
template int test_rk_ts<double>(int, double, double, double);
} //end namespace libpetey