home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / fixptlib / mnewt.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-27  |  2.5 KB  |  97 lines
  1. /*
  2. Compute periodic orbit of a period r by newton's method
  3. (Solve F(x)=f^r(x)-x-xoffset= 0)
  4. ------------------------------------------------------
  5. Note:    x[0...n-1] convention 
  6.     
  7. Input:    x[n]= current guess, n= phase space dimension, r= period
  8.     tolx= tolerated error in x (sum), tolf= tolerated error in F(x) (sum)
  9.     ntrial = # of maximum secant step
  10. Output:    x[n]=new guessed solution, *ndid= acutal # of secant step done,
  11.     *err= actual error in x (sum)
  12. */
  13.  
  14. #define DELTAX 1.e-8
  15. #define FREERETURN {*ndid = k;free_dmatrix(alpha,0,n-1,0,n-1);free_dvector(x2,0,n-1);free_dvector(xt,0,n-1);free_dvector(beta,0,n-1);\
  16.     free_ivector(indx,0,n-1);return;}
  17.  
  18. void mnewt(x,ntrial,n,r,tolx,tolf,err,ndid)
  19. int ntrial,n,r,*ndid;
  20. double x[],tolx,tolf,*err;
  21. {
  22.     int i,j,k,*indx,*ivector();
  23.     double tolx10,fabs(),errf,d,*x2,*xt,*beta,**alpha,*dvector(),**dmatrix();
  24.     void usrfun(),usrfun2(),usrfun3(),ludcmp(),lubksb(),free_dmatrix(),free_ivector(),free_dvector();
  25.     extern int stop,enable_period,cur_color,model,mapping_on,fderiv_on;
  26.  
  27.     tolx10 = tolx /10;
  28.     indx = ivector(0,n-1);
  29.     x2 = dvector(0,n-1);
  30.     xt = dvector(0,n-1);
  31.     beta = dvector(0,n-1);
  32.     alpha = dmatrix(0,n-1,0,n-1);
  33.     for(k=0;k<ntrial;k++){
  34.         for(i=0;i<n;i++) x2[i] = x[i] + DELTAX;
  35.         if(mapping_on){
  36.             if(fderiv_on)
  37.                 usrfun(x,alpha,beta,r,n);
  38.             else
  39.                 usrfun2(x,x2,alpha,beta,r,n);
  40.         }
  41.         else {
  42.             usrfun3(x,x2,alpha,beta,r,n);
  43.         }
  44.         if(enable_period==1){
  45.             dist_periodic(xt,beta,n);
  46.             for(i=0;i<n;i++)beta[i]=xt[i];
  47.         }
  48.         /*
  49.         printf("NTRIAL=%d\n",k);
  50.         printf("usrfun x: ");
  51.         for(i=0;i<n;i++) printf("%f ",x[i]);    
  52.         printf("\n");
  53.         printf("usrfun f: ");
  54.         for(i=0;i<n;i++) printf("%f ",beta[i]);    
  55.         printf("\n");
  56.         printf("usrfun alpha:");
  57.         for(j=0;j<n;j++){
  58.             for(i=0;i<n;i++) printf("%f ",alpha[j][i]);    
  59.             printf("\n");
  60.         }
  61.         */
  62.  
  63.         errf = 0.0;
  64.         for(i=0;i<n;i++) errf += fabs(beta[i]);
  65.         if(errf <= tolf) FREERETURN
  66.  
  67.         ludcmp(alpha,n,indx,&d);
  68.         /* should be placed here not to proceed further
  69.         in case of singular matrix (singular matrix often
  70.         arise when all the machine accuracy was lost in
  71.         computing alpha. This happens in case of good
  72.         convergence,too) */
  73.         if(stop){
  74.             FREERETURN
  75.         }
  76.         lubksb(alpha,n,indx,beta);
  77.  
  78.         /*
  79.         printf("new beta:");
  80.         for(i=0;i<n;i++) printf("%f ",beta[i]);    
  81.         printf("\n");
  82.         */
  83.         *err = 0.0;
  84.         for(i=0;i<n;i++){
  85.             *err += fabs(beta[i]);
  86.             x[i] -= beta[i];
  87.         }
  88.         if(*err <= tolx) FREERETURN
  89.         for(i=0;i<n;i++) if(fabs(beta[i]) < tolx10) x[i] -= tolx10;
  90.             
  91.         /* draw intermediate steps */
  92.         (void) draw_record_pwf(x,cur_color,1,3,1,0);
  93.         (void) draw_record_pwf(x,cur_color,1,3,1,0);
  94.     }
  95.     FREERETURN
  96. }
  97.