home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / general / fractal / kaos.lha / fixptlib / usrfun3.c < prev    next >
Encoding:
C/C++ Source or Header  |  1990-01-27  |  1.4 KB  |  57 lines
  1. /*
  2. Compute Jacobian and value of nonlinear eq (f^r(x) - x - xoffset)
  3. by finite difference method
  4. -----------------------------
  5. NOTE:    nonlinear equations are not handled in periodic coords
  6. Input:    r=period of an orbit, n=phase space dimension,
  7.     x[n],x2[n] = coords of two current guesses
  8. Output:    alpha[n][n]=DF(x) (Jacobian), beta[n]=F(x)
  9. */
  10.  
  11. void usrfun3(x,x2,alpha,beta,r,n)
  12. double *x,*x2,**alpha,*beta;
  13. int r,n;
  14. {
  15.     int i,k;
  16.     extern int enable_period;
  17.     extern double time,*t_va,*t_vf,*xoffset;
  18.  
  19.     /* time is set to 0: Fix this later for other problems */
  20.     time = 0;
  21.     /* xoffset need to be installed later */
  22.     for(i=0;i<n;i++) xoffset[i] = 0;
  23.  
  24.     /* F(x) */
  25.     for(i=0;i<n;i++) t_vf[i] = x[i];
  26.     fode_user(t_vf,time,r,n);
  27.     /* solve F(x)=0 if r=0 otherwise x(r * T + t) - x(t) */
  28.     if(r==0)
  29.         for(i=0;i<n;i++) t_va[i] = t_vf[i] - xoffset[i];
  30.     else
  31.         for(i=0;i<n;i++) t_va[i] = t_vf[i] - x[i] - xoffset[i];
  32.  
  33.     /* DF(x) by finite difference */
  34.     for(k=0;k<n;k++){
  35.         for(i=0;i<n;i++) t_vf[i] = x[i];
  36.         t_vf[k] = x2[k];
  37.         /* r=0: t_vf= f(x), r>=1: t_vf = t_vf(time+T*r) */
  38.         fode_user(t_vf,time,r,n);
  39.         /* For equilibrium, dist = f(x) = 0 */
  40.         if(r==0) {
  41.             /* no offset */
  42.         }
  43.         /* For periodic orbits, dist = f^r(x) - x - xoffset */
  44.         else {
  45.             for(i=0;i<n;i++){
  46.                 if(i==k)
  47.                     t_vf[i] -= (x2[i] + xoffset[i]);
  48.                 else
  49.                     t_vf[i] -= (x[i] + xoffset[i]);
  50.             }
  51.         }
  52.         for(i=0;i<n;i++) alpha[i][k] = (t_vf[i]-t_va[i])/(x2[k]-x[k]);
  53.     }
  54.  
  55.     for(i=0;i<n;i++)beta[i]=t_va[i];
  56. }
  57.