home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / bsstep.c next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  1.7 KB  |  70 lines
  1. /*
  2. ### Bulirsch-Stoer integration one stepper ###
  3.  
  4. Note:    a slightly modified version of bsstep in "Numerical Recipes"
  5.     c convention of [0,n-1] is used for array indices
  6. */
  7.  
  8. #include <math.h>
  9. #define IMAX 11
  10. #define NUSE 7
  11. #define SHRINK 0.95
  12. #define GROW 1.2
  13.  
  14. extern int_option;
  15. double **ddd = 0, *xxx = 0;
  16.  
  17. void bsstep(y,dydx,nv,xx,htry,eps,yscal,hdid,hnext)
  18. double y[],dydx[],*xx,htry,eps,yscal[],*hdid,*hnext;
  19. int nv;
  20. {
  21.     int i,j;
  22.     double xsav,xest,h,errmax,errnew,temp;
  23.     double *ysav,*dysav,*yseq,*yerr,*dvector(),**dmatrix();
  24.     static int nseq[IMAX]={2,4,6,8,12,16,24,32,48,64,96};
  25.     void free_dmatrix(),free_dvector();
  26.  
  27.     ysav = dvector(0,nv-1);
  28.     dysav = dvector(0,nv-1);
  29.     yseq = dvector(0,nv-1);
  30.     yerr = dvector(0,nv-1);
  31.     xxx = dvector(0,IMAX-1);
  32.     ddd = dmatrix(0,nv-1,0,NUSE-1);
  33.     h = htry;
  34.     xsav = (*xx);
  35.     for(i=0;i<nv;i++) {
  36.         ysav[i] = y[i];
  37.         dysav[i] = dydx[i];
  38.     }
  39.     for(;;){
  40.         for(i=0;i<IMAX;i++){
  41.             (void) mmid(ysav,dysav,nv,xsav,h,nseq[i],yseq,0);
  42.             xest = (temp = h / nseq[i] , temp * temp); 
  43.             (void) rzextr(i,xest,yseq,y,yerr,nv,NUSE);
  44.             errmax = 0.0;
  45.             for(j=0;j<nv;j++){
  46.                 if(errmax < (errnew = fabs(yerr[j] / yscal[j])))
  47.                     errmax = errnew;
  48.             }
  49.             if(errmax < eps) {
  50.                 if(int_option)
  51.                     (void) mmid(ysav,dysav,nv,xsav,h,nseq[i],yseq,1);
  52.                 *xx += h;
  53.                 *hdid = h;
  54.                 *hnext = i == NUSE-1 ? h * SHRINK : 
  55.                     i == NUSE-2 ? h * GROW : (h * nseq[NUSE-2]) / nseq[i];
  56.                 free_dmatrix(ddd,0,nv-1,0,NUSE-1);
  57.                 free_dvector(xxx,0,IMAX-1);
  58.                 free_dvector(yerr,0,nv-1);
  59.                 free_dvector(yseq,0,nv-1);
  60.                 free_dvector(dysav,0,nv-1);
  61.                 free_dvector(ysav,0,nv-1);
  62.                 return;
  63.             }
  64.         }
  65.         h *= 0.25;
  66.         for(i=0;i < (IMAX - NUSE)/2;i++) h /= 2.0;
  67.         if((*xx + h) == (*xx)) system_mess_proc(1,"bsstep: step size underflow. Stop!");
  68.     }
  69. }
  70.