1

I am trying to implement in C the algorithm by Parker and Chua as described in the book "Practical Numerical Algorithms for Chaotic Systems" (chapter 6).

In that chapter, the authors describe a method to plot manifolds for discrete time systems (aka maps). I am trying to translate the pseudocode given into a C code for standard map (Chirikov-Taylor). I already did it for Fortran and now C gives the same issue as Fortran 90.

The idea is to display points on (unstable, for while) manifold such that the distance between them are not far from a given tolerance. If so, then points are linearly interpolated. As authors says, like a window moving and revealing the manifold.

The problem is that the output gives some points on manifold, but after few points the output hangs at one value.

I am warned about C and Fortran, specially about arrays indexation. But I can't see what is wrong with the method or my translation from pseudocode...

I know there are plenty of methods for manifolds visualization and papers about it, but if I am not able to construct one from a given pseudocode, I doubt if my other implementations are in the right way... so I prefer to make it work.

Well, let me show a simplified version of my code, that should work (in C):

  x_bar = alpha * vx;
  y_bar = alpha * vy;

  //choose linear approximation as initial condition
  x[0] = x_bar;
  p[0] = y_bar;

  //iterate once
  x[1] = map_x(x[0], p[0]);
  p[1] = map_p(x[0], p[0]);

  //iterate again
  x[2] = map_x(x[1], p[1]);
  p[2] = map_p(x[1], p[1]);

  x_new = x[1];
  y_new = p[1];

  do
    {
      dist = sqrt( pow((x[2] - x[1]), 2.) + pow((p[2] - p[1]), 2.) ); 
      
      if(dist < 0.001) //absolute criteria for while...
    {
      //accept the point and move the window
      printf("%f\t %f\n", x[0], p[0]);
      
      x[0] = x[1];
      p[0] = p[1];
      
      x[1] = x[2];
      p[1] = p[2];
      
      i++;
    }
      else
    {
      //linear interpolation required
      x_new = (x[0] + x_new) / 2.;
      y_new = (p[0] + y_new) / 2.;
      
      x_aux = map_x(x_new, y_new);
      y_aux = map_p(x_new, y_new);
      
      x[2] = x_aux;
      p[2] = y_aux;
    }
    }
  while(i < 10);

In the code above, vx and vy are the eigenvectors of the Jacobian matrix, evaluated at the fixed point (x, p) = (0, 0). My functions 'map_x' and 'map_p' are working well and returns the one step iteration of the Chirikov-Taylor map.

I attached a figure produced by first steps of my code. Consider r = (x, p).

Fig: first points produced by the algorithm

Thank you in advance.

baseaxis
  • 11
  • 1
  • While probably not the source of your issues, the expression you use to calculate `dist` would be better served by the [**hypot**](https://en.cppreference.com/w/c/numeric/math/hypot) function for numerical stability. Also, it would be better to use a less 'arbitrary' epsilon value, e.g., `if (dist < FLT_EPSILON) ...` or `DBL_EPSILON`, as defined in: `` – Brett Hale Aug 30 '20 at 08:20
  • Thank you very much, @BrettHale. I really appreciated your suggestions. However, to be perfectly honest, it seems that the map require that $x^2 + y^2$ is used instead of distance, because of the variable update. For that, I chose hypot * hypot now to see how things run. Concerning the epsilon issue, it is a nice suggestion also, but the pseudocode I am studying says that relative error and absolute error are user choices. Furthermore, the authors use a combination of relative/absolute criteria, like: $d < L * E_r + E_a$, being $E_r$ the relative error, $E_a$ the absolute error and $L$ a norm. – baseaxis Sep 02 '20 at 15:13

1 Answers1

1

It seems that if the dist < 0.001 you aren't updating the x[2] or p[2] but you are setting x[1]=x[2] and p[1]=p[2]. This means the next time it loops the distance will be 0 because you're doing x[2]-x[1] and p[2]-p[1] but they're equal to each other so you'll get dist = 0. And 0 < 0.001 so you'll just keep repeating this and be stuck at the same point.

You're probably missing something like:

if (dist < 0.001) {
// other stuff
x[2] = map_x(x[1], p[1]);
p[2] = map_p(x[1], p[1]);
} else { // ...
g23
  • 666
  • 3
  • 9
  • Thank you very much, @g23. Despite your great answer, the problem is with else and not with if. When I print the value of distance at each iteration, this distance hangs at the same value. My hypothesis is that, since the variable is being updated, they act like a discrete map with a convergence. This is a problem, because it was a translation from the pseudocode. Still trying... But thank you very much anyway. Your answer made me look all the loops and found some issues. – baseaxis Sep 02 '20 at 16:11