home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / complib / symp_sqze.c < prev    next >
Encoding:
C/C++ Source or Header  |  1989-11-18  |  3.9 KB  |  162 lines
  1. /*
  2. ### solve implicit nonlinear equation by secant's method ###
  3.  
  4. Input:    (double *) pxv (dim=var_dim): current coords
  5.     (double *) pvx1, *pvx2 (dim=var_dim/2): two guesses
  6.         (only coords parts of the phase space variables are used here).
  7.     (int) dim: the dimension of the phase space
  8. Output:    (double *) pvx1 (dim=var_dim): returns the coordinates of the next
  9.         orbit.
  10. */
  11. symp_sqze(pvx1,pvx2,pvx,time,time_step,dim)
  12. double *pvx1,*pvx2,*pvx,time,time_step;
  13. int dim;
  14. {
  15.     int i,j,it,max_secant,index,iswit,ii,i2,dim2,identity;
  16.     double eps,epsf,err;
  17.     double **vnew,**vjac,*dvec,*dvecn,**fiter,*vfull,*vfull2,*dvector(),**dmatrix();
  18.     double fabs(),log10();
  19.     extern int stop;
  20.     extern double cutoff,*param;
  21.     extern int (*f_p)();
  22.  
  23.     
  24.     dim2 = dim/2;
  25.     vnew = dmatrix(0,dim2-1,0,1);
  26.     vjac = dmatrix(0,dim2-1,0,dim2-1);
  27.     fiter = dmatrix(0,dim2-1,0,dim2);
  28.     dvec = dvector(0,dim2-1);
  29.     dvecn = dvector(0,dim2-1);
  30.     vfull = dvector(0,dim-1);
  31.     vfull2 = dvector(0,dim-1);
  32.     eps = 1.e-14;
  33.     epsf = 1.e-16;
  34.     max_secant = 21;
  35.     for(i=0;i<dim2;i++){
  36.         i2 = i*2;
  37.         vnew[i][0]= *(pvx1+i2);
  38.         vnew[i][1]= *(pvx2+i2);
  39.     }
  40.     /*
  41.     it=0;
  42.     printf("%%v%d x1,x2=%18.10e %18.10e z1,z2=%18.10e %18.10e\n",it,vnew[0][0],vnew[0][1],vnew[1][0],vnew[1][1]);
  43.     */
  44.     for(it = 0;it <max_secant;it++){
  45.         for(i=0;i<dim2;i++){
  46.             if(vnew[i][1] > cutoff || vnew[i][1] < - cutoff) {
  47.                 system_mess_proc(1,"Orbits appear to diverge off to an infinity! Stop!");
  48.                 stop=1;
  49.                 goto done;
  50.             }
  51.         }
  52.         for(index=0;index<dim2+1;index++){
  53.             for(i=0;i<dim2;i++){
  54.                 i2=i*2;
  55.                 vfull[i2]=vnew[i][0];
  56.                 vfull[i2+1] = *(pvx+i2+1);
  57.             }
  58.             if(index != dim2)
  59.                 vfull[2*index] = vnew[index][1];
  60.             /*
  61.             printf("%%vfull[%d] %18.12f %18.12f %18.12f %18.12f\n",index,vfull[0],vfull[1],vfull[2],vfull[3]);
  62.             */
  63.             (int) f_p(vfull2,1,vfull,param,time,dim);
  64.             /*
  65.             printf("%%vfull2[%d] %18.12f %18.12f %18.12f %18.12f\n",index,vfull2[0],vfull2[1],vfull2[2],vfull2[3]);
  66.             */
  67.             for(i=0;i<dim2;i++){
  68.                 i2=i*2;
  69.                 fiter[i][index]= vfull[i2] - *(pvx+i2) - time_step * vfull2[i2]; 
  70.             }
  71.             /*
  72.             printf("%%fiter[][%d] %18.12f %18.12f\n",index,fiter[0][index],fiter[1][index]);
  73.             */
  74.         }
  75.         for(i=0;i<dim2;i++){
  76.             for(j=0;j<dim2;j++)
  77.                 vjac[i][j]=(fiter[i][j]-fiter[i][dim2])/(vnew[j][1]-vnew[j][0]);
  78.             dvec[i]= -fiter[i][dim2];
  79.         }
  80.         printf("%%vjac %18.12f %18.12f %18.12f %18.12f\n",vjac[0][0],vjac[0][1],vjac[1][0],vjac[1][1]);
  81.         /*
  82.         printf("%%dvec %18.12f %18.12f\n",dvec[0],dvec[1]);
  83.         */
  84.  
  85.         /* see if the jacobian is identity */
  86.         identity = 1;
  87.         for(i=0;i<dim2;i++){
  88.             for(j=0;j<dim2;j++){
  89.                 if(i=j){
  90.                     if(vjac[i][j]!=1.){
  91.                         identity=0;
  92.                         break;
  93.                     }
  94.                 }
  95.                 else if (vjac[i][j]!=0.){
  96.                     identity = 0;
  97.                     break;
  98.                 }
  99.             }
  100.             if(identity==0)
  101.                 break;
  102.         }
  103.  
  104.         if(identity){
  105.             system_mess_proc(1,"symp_sqze: Jacobian=Identity. Use explicit symp integration!");
  106.             stop = 1;
  107.             return;
  108.         }
  109.         
  110.         stop = (int) linsol(dvecn,vjac,dvec,dim2);
  111.  
  112.         if(stop) {
  113.             system_mess_proc(1,"symp_sqze: sigular matrix");
  114.             goto done;
  115.         }
  116.         /*
  117.         printf("dvecn %18.10e %18.10e dvec %18.10e %18.10e\n",dvecn[0],dvecn[1],dvec[0],dvec[1]);
  118.         */
  119.     
  120.         /* testing for adequacy of roots */
  121.         iswit =1;
  122.         for(i=0;i<dim2;i++){
  123.             if(fabs(dvecn[i]) >eps)
  124.                 iswit = 0;
  125.             vnew[i][1]=vnew[i][0];
  126.             vnew[i][0]=vnew[i][0]+dvecn[i];
  127.         }
  128.         if(iswit)
  129.             goto done;
  130.         iswit =1;
  131.         for(i=0;i<dim2;i++){
  132.             if( fabs(dvec[i]) > epsf)
  133.                 iswit = 0;
  134.         }
  135.         if(iswit)
  136.             goto done;
  137.         /*
  138.         printf("%%v%d %18.10e %18.10e %18.10e %18.10e\n",it,vnew[0][0],vnew[0][1],vnew[1][0],vnew[1][1]);
  139.         */
  140.         printf("\n");
  141.     }
  142.  
  143. done:
  144.     for(i=0;i<dim2;i++){
  145.         i2 = i*2;
  146.         *(pvx1+i2) = vnew[i][0];
  147.         *(pvx2+i2) = vnew[i][1];
  148.     }
  149.     err = 1.e-16;
  150.     for(i=0;i<dim2;i++){
  151.         if(err > dvec[i])
  152.                   err = dvec[i];
  153.     }
  154.     (void) free_dmatrix(vnew,0,dim2-1,0,1);
  155.     (void) free_dmatrix(vjac,0,dim2-1,0,dim2-1);
  156.     (void) free_dmatrix(fiter,0,dim2-1,0,dim2);
  157.     (void) free_dvector(dvec,0,dim2-1);
  158.     (void) free_dvector(dvecn,0,dim2-1);
  159.     (void) free_dvector(vfull,0,dim-1);
  160.     (void) free_dvector(vfull2,0,dim-1);
  161. }
  162.