home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / symp_picard.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  1.7 KB  |  74 lines
  1. /*
  2. ### Picard iteration to solve nonlinear implicit equation ###
  3.  
  4. Input:    (double *) pxv (dim=var_dim): current coords
  5.     (int) dim: the dimension of the phase space
  6. Output:    (double *) pvy (dim=var_dim): returns the coordinates of the next
  7.         orbit.
  8. */
  9. symp_picard(pvy,pvx,time,time_step,dim)
  10. double *pvy,*pvx,time,time_step;
  11. int dim;
  12. {
  13.     int i,it,max_iter,iswit;
  14.     double eps,epsf,err,fabs(),*vnew,*vnewt,*dvec,*vnewp,*dvector();
  15.     extern int stop,nok,nbad;
  16.     extern double cutoff,*param;
  17.     extern int (*f_p)();
  18.  
  19.     vnew = dvector(0,dim-1);
  20.     vnewp= dvector(0,dim-1);
  21.     vnewt= dvector(0,dim-1);
  22.     dvec = dvector(0,dim-1);
  23.     eps = 1.e-14;
  24.     epsf = 1.e-16;
  25.     max_iter = 10;
  26.     for(i=0;i<dim;i++) vnew[i]= *(pvx+i);
  27.     for(i=0;i<dim;i++) vnewp[i]= *(pvx+i);
  28.  
  29.     for(it = 0;it <max_iter;it++){
  30.         for(i=0;i<dim;i++){
  31.             if(vnew[i] > cutoff || vnew[i] < - cutoff) {
  32.                 system_mess_proc(1,"Orbits appear to diverge off to an infinity! Stop!");
  33.                 stop=1;
  34.                 goto done;
  35.             }
  36.         }
  37.         for(i=0;i<dim;i++) vnewt[i] = (vnew[i] + *(pvx+i))/2;
  38.         (int) f_p(vnew,0,vnewt,param,time,dim);
  39.         for(i=0;i<dim;i++) vnew[i] = *(pvx+i) + time_step * vnew[i];
  40.  
  41.         if(stop) {
  42.             system_mess_proc(1,"symp_sqze: sigular matrix");
  43.             goto done;
  44.         }
  45.     
  46.         /* testing for adequacy of roots */
  47.         iswit =1;
  48.         for(i=0;i<dim;i++){
  49.             dvec[i] = vnew[i]-vnewp[i];
  50.             if(fabs(dvec[i]) >eps)
  51.                 iswit = 0;
  52.             vnewp[i]=vnew[i];
  53.         }
  54.         if(iswit)
  55.             goto done;
  56.     }
  57.  
  58. done:
  59.     for(i=0;i<dim;i++) *(pvy+i) = vnew[i];
  60.     err = 1.e-16;
  61.     for(i=0;i<dim;i++){
  62.         if(err > dvec[i])
  63.                   err = dvec[i];
  64.     }
  65.     if(it>=10)
  66.         nbad++;
  67.     else
  68.         nok++;
  69.     (void) free_dvector(vnew,0,dim-1);
  70.     (void) free_dvector(vnewt,0,dim-1);
  71.     (void) free_dvector(vnewp,0,dim-1);
  72.     (void) free_dvector(dvec,0,dim-1);
  73. }
  74.