/home/peifer/diffit/modules/integrateRK.cc

Go to the documentation of this file.
00001 #include <math.h>
00002 #include <iostream>
00003 #include "../nr.h"
00004 #include "../def.h"
00005 
00006 using namespace std;
00007 
00008 #define TINY 1.0e-30
00009 
00010 /* extension by w.h. 1/98: f means full
00011 if ode_fxp!=NULL, ALL intermediate x-values will be stored there.
00012 Then we can calculate derivatives w.resp. to parameters or initial values
00013 by integrating the curve with the same steps and a light routine like rk4.
00014 We need no intermediate y-values for this purpose => We don't want to set
00015 ode_dxsav=0 and use ode_xp. 
00016 If desired, allocate ode_fxp[1..ode_fkmax] in the calling program. */
00017 
00018 
00019 void integrateRK(Glob *globs,double ystart[], double p[], long nvar, double x1, double x2,double eps,
00020                  double h1, double hmin, long *nok, long *nbad,
00021                  void (*derivs) (double, double[], double[], double[]))
00022 {
00023   long ode_fkmax = 0, ode_fkount = 0, ode_kmax = 0, ode_kount = 0;
00024   double ode_dxsav = 0, *ode_fxp = 0, *ode_xp = 0, **ode_yp = 0;
00025   long rkqs_ign = 0;
00026   long nstp,i;
00027   double xsav,x,hnext,hdid,h;
00028   double *yscal,*y,*dydx;
00029   long odeint_MAXSTP=globs->maxstp;
00030 
00031   yscal=dvector(1,nvar);
00032   y=dvector(1,nvar);
00033   dydx=dvector(1,nvar);
00034   x=x1;
00035   h=SIGN(h1,x2-x1);
00036   *nok = (*nbad) = ode_kount = 0;
00037   for (i=1;i<=nvar;i++) y[i]=ystart[i];
00038   if (ode_dxsav<0) nrerror("odeint:dxsav<0");
00039   if (ode_kmax>0 && ode_kount<ode_kmax-1) {
00040      ode_xp[++ode_kount]=x;
00041      for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
00042      xsav=x;
00043   } /* save first point, w.h. 19.03.99 */
00044   for (nstp=1;nstp<=odeint_MAXSTP;nstp++) {
00045     (*derivs)(x,y,dydx,p);
00046     for (i=1;i<=nvar;i++)
00047       yscal[i]=fabs(y[i])+fabs(dydx[i]*h)+TINY;
00048     if (ode_fkmax) {
00049       if (ode_fkount>=ode_fkmax) nrerror("ode_fkount exceeds ode_fkmax");
00050       ode_fxp[++ode_fkount]=x;
00051     }
00052     if (ode_kmax>0 && ode_kount<ode_kmax-1 && fabs(x-xsav)>ode_dxsav) {
00053       ode_xp[++ode_kount]=x;
00054       for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
00055       xsav=x;
00056     }
00057     if ((x+h-x2)*(x+h-x1) > 0.0) h=x2-x;
00058     rkqs(y,dydx,p,nvar,&x,h,eps,yscal,&hdid,&hnext,globs->rkqs_ign,derivs);
00059     if (hdid == h) ++(*nok); 
00060     else ++(*nbad);
00061     if ((x-x2)*(x2-x1) >= 0.0) {
00062       for (i=1;i<=nvar;i++) ystart[i]=y[i];
00063       if (ode_kmax>0 && ode_kount<ode_kmax-1) {
00064         ode_xp[++ode_kount]=x;
00065         for (i=1;i<=nvar;i++) ode_yp[i][ode_kount]=y[i];
00066       }
00067       free_dvector(dydx,1,nvar);
00068       free_dvector(y,1,nvar);
00069       free_dvector(yscal,1,nvar);
00070       return;
00071     }
00072     if (fabs(hnext) <= hmin)
00073       {
00074         cerr << "Step size too small in odeint\n";
00075         throw 1;
00076       }
00077     h=hnext;
00078   }
00079   cerr << "Too many steps in routine odeint; try -maxstp <n> greater n\n";
00080   throw 1;
00081 }
00082 #undef TINY

Generated on Mon Jan 29 17:09:12 2007 for Diffit by  doxygen 1.4.6