/* ************************************************************************ * runge.cpp: 4th order Runge-Kutta solution for harmonic oscillator * * In this case, d^y2/dx^2 = - y, can be written as dy_1/dx = y_2 and * * dy_2 /dx = - y_1, with y_1 = y, and y_2 = dy/dx. * C.A. Bertulani - 03/20/2000 */ /**************************************************************/ #include #include #include #include #include typedef double Number; #include "butil.h" #include #define N 2 /* number of equations */ #define dist 0.1 /* stepsize */ #define MIN 0.0 /* minimum x */ #define MAX 10.0 /* maximum x */ void main() { void f(Number x, Number y[], Number dy[]); void rkqs(Number [], Number [], int, Number *, Number, Number, Number [], Number *, Number *, void (*f)(Number, Number [], Number [])); void runge(Number y[], Number dydx[], int n, Number x, Number step, Number yout[], void (*f)(Number, Number [], Number [])); void odeint(Number ystart[], int nvar, Number x1, Number x2, Number eps, Number h1, Number hmin, void (*derivs)(Number, Number [], Number []), void (*rkqs)(Number [], Number [], int, Number *, Number, Number, Number [], Number *, Number *, void (*)(Number, Number [], Number []))); Number x, eps, *y, *dy, *yout; y = vector(N); dy = vector(N); yout = vector(N); y[1] = 1.0; /* initial position */ y[2] = 0.0; /* initial velocity */ /* Test with a simple Runge-Kutta - Follow step by step till x=MAX */ ofstream out("runge1.dat"); /* save data in runge1.dat */ for(x = MIN; x <= MAX ; x += dist) { f(x, y, dy); /* get derivatives */ runge(y, dy, N, x, dist, y, f); out << setw(3) << x << setw(15) << setprecision(4) << y[1]; out << setw(15) << setprecision(4) << cos(x+dist) << endl; } cout << "data stored in runge1.dat\n"; /* Runge-Kutta with adaptive step control - Integrate to x=MAX */ ofstream out2("runge2.dat"); /* save data in runge2.dat */ y[1] = 1.0; /* initial position */ y[2] = 0.0; /* initial velocity */ eps=1.e-6; /* Accuracy */ odeint(y, N, MIN, MAX, eps, dist, 0.0, f, rkqs); out2 << setw(3) << MAX << setw(15) << setprecision(4) << y[1]; out2 << setw(15) << setprecision(4) << cos(MAX) << endl; cout << "data stored in runge2.dat\n"; } #undef dist /* stepsize */ #undef MIN /* minimum x */ #undef MAX /* maximum x */ #undef N /* number of equations */ /*-----------------------end of main program--------------------------*/ /************************************************************************ /* definition of equations - this is the harmonic oscillator */ /* derivatives: y_1 and y_2 */ void f(Number x, Number y[], Number dy[]) { dy[1]=y[2]; /* RHS of first equation */ dy[2]=-y[1]; /* RHS of second equation */ } /*******************************************************************************/ /******************************************************************************* * 4th order Runge-Kutta subroutine Given values for the variables y[1..n] and their derivatives dydx[1..n] known at x, advance solution over an interval h and return the incremented variables as yout[1..n], which need not be a distinct array from y. The user supplies the routine derivs(x,y,dydx), which returns derivatives dydx at x. C.A. Bertulani - 03/20/2000 * *******************************************************************************/ void runge(Number y[], Number dydx[], int n, Number x, Number h, Number yout[], void (*derivs)(Number, Number [], Number [])) { int i; Number xh,hh,h6,*dym,*dyt,*yt; dym=vector(n); dyt=vector(n); yt=vector(n); hh=h*0.5; h6=h/6.0; xh=x+hh; for (i=1;i<=n;i++) yt[i]=y[i]+hh*dydx[i]; /* First step */ (*derivs)(xh,yt,dyt); /* Second step */ for (i=1;i<=n;i++) yt[i]=y[i]+hh*dyt[i]; (*derivs)(xh,yt,dym); /* Third step */ for (i=1;i<=n;i++) { yt[i]=y[i]+h*dym[i]; dym[i] += dyt[i]; } (*derivs)(x+h,yt,dyt); /* Fourth step */ for (i=1;i<=n;i++) /* Accumulate increments with proper weights */ yout[i]=y[i]+h6*(dydx[i]+dyt[i]+2.0*dym[i]); delete []yt; delete []dyt; delete []dym; } /******************************************************************************* * routine odeint * 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 by values at the end of the integration interval. derivs is the user-supplied routine for calculating the rigt-hand side derivative, while rkqs is teh name of the stepper routine to be used. * C.A. Bertulani - 03/20/2000 * *******************************************************************************/ #define MAXSTP 10000 #define TINY 1.0e-30 void odeint(Number ystart[], int nvar, Number x1, Number x2, Number eps, Number h1, Number hmin, void (*derivs)(Number, Number [], Number []), void (*rkqs)(Number [], Number [], int, Number *, Number, Number, Number [], Number *, Number *, void (*)(Number, Number [], Number []))) { int nstp,i; Number x,hnext,hdid,h; Number *yscal,*y,*dydx; yscal=vector(nvar); y=vector(nvar); dydx=vector(nvar); x=x1; h=SIGN(h1,x2-x1); for (i=1;i<=nvar;i++) y[i]=ystart[i]; for (nstp=1;nstp<=MAXSTP;nstp++) { /* Take at most MAXSTP steps */ (*derivs)(x,y,dydx); for (i=1;i<=nvar;i++) /* Scaling to monitor accuracy. This general-purpose choice */ /* can be modified if needed. */ yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY; 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 ((x-x2)*(x2-x1) >= 0.0) { /* Are we done? */ for (i=1;i<=nvar;i++) ystart[i]=y[i]; delete []dydx; delete []y; delete []yscal; return; /* normal exit */ } if (fabs(hnext) <= hmin) cout << "Step size too small in odeint"; h=hnext; } cout << "Too many steps in routine odeint"; } #undef MAXSTP #undef TINY /******************************************************************************* * routine rkqs Used in adaptive size Runge-Kutta integration. 5th order Runge-Kutta step with monitoring of local truncation error to ensure acuracy and adjust stepsize. Input are the independent variable vector y[1..n] and its derivative dydx[1..n] at the starting value of the independent variable x. Also input are 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 by their 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 size derivatives. C.A. Bertulani - 03/20/2000 * *******************************************************************************/ #define SAFETY 0.9 #define PGROW -0.2 #define PSHRNK -0.25 #define ERRCON 1.89e-4 /*The value ERRCON equals (5/SAFETY) raised to the power (1/PGROW).*/ void rkqs(Number y[], Number dydx[], int n, Number *x, Number htry, Number eps, Number yscal[], Number *hdid, Number *hnext, void (*derivs)(Number, Number [], Number [])) { void rkck(Number y[], Number dydx[], int n, Number x, Number h, Number yout[], Number yerr[], void (*derivs)(Number, Number [], Number [])); int i; Number errmax,h, htemp, xnew,*yerr,*ytemp; yerr=vector(n); ytemp=vector(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 the required accuracy */ 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) cout << "stepsize underflow in rkqs"; } if (errmax > ERRCON) *hnext=SAFETY*h*pow(errmax,PGROW); else *hnext=5.0*h; /* No more than a factor 5 of increase */ *x += (*hdid=h); for (i=1;i<=n;i++) y[i]=ytemp[i]; delete []ytemp; delete []yerr; } #undef SAFETY #undef PGROW #undef PSHRNK #undef ERRCON /******************************************************************************* * routine rkck Used in adaptive size Runge-Kutta integration. Given values for the variables y[1..n] and their derivatives dydx[1..n] known at x, advance solution over an interval h and return the incremented variables as yout[1..n]. Also returns an estimate of the local truncation error in yout using embedded fourth-order method. The user supplies the routine derivs(x,y,dydx), which returns derivatives dydx at x. C.A. Bertulani - 03/20/2000 * *******************************************************************************/ void rkck(Number y[], Number dydx[], int n, Number x, Number h, Number yout[], Number yerr[], void (*derivs)(Number, Number [], Number [])) { int i; static Number 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.0/14336.0; Number dc1=c1-2825.0/27648.0,dc3=c3-18575.0/48384.0, dc4=c4-13525.0/55296.0,dc6=c6-0.25; Number *ak2,*ak3,*ak4,*ak5,*ak6,*ytemp; ak2=vector(n); ak3=vector(n); ak4=vector(n); ak5=vector(n); ak6=vector(n); ytemp=vector(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++) /* Acumulate 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 differences between fourth and fifth order methods */ delete []ytemp; delete []ak6; delete []ak5; delete []ak4; delete []ak3; delete []ak2; } /*******************************************************************************/