Strange Limit Cycle: Unveiling a 3D Dynamics

[Tom Piper]

The 3D dynamical system introduced in the following site;

http://www.geocities.com/tontokohirorin/mathematics/limitcycle/limitcycle2.htm

shows not only chaotic behaviour, but also converges to a strange-shaped limit cycle under a certain parameter setting!

 

[saltydog]

The 3D dynamical system introduced in the following site;

http://www.geocities.com/tontokohirorin/mathematics/limitcycle/limitcycle2.htm

shows not only chaotic behaviour, but also converges to a strange-shaped limit cycle under a certain parameter setting!

Hello Tom. I think that's quite consistent with deterministic chaos in general: a chaotic system has dense periodic points (infinite number of limit cycles) and depending on the parameters chosen for the trajectory, the dynamics will exhibit a chaotic (strange) attractor or collapse (attract to) one of the many limit cycles. This is exhibited by the corresponding Feigenbaum plot of the system which shows where the limit cycles are and where the chaotic regime lies. For example, the Lorenz and Rossler attractors also have limit cycles depending on the parameters chosen.

 

[Tom Piper]

I agree with you. Generally in dynamical system the chaotic behavior will be changed to the recurrent motion converging to a simple limit cycle by tuning parameters. However if my calculation is correct, the Lyapunov exponent on the "limit cycle" appeared on the dynamical system shown in my webpage takes a certain positive value, i.e. that "limit cycle" is unstable! Actually any integral curve starting from the vicinity to that limit cycle diverges locally. Nevertheless globally that integral curve converges to that limit cycle. It seems very strange (at least) for me...

 

[saltydog]

Let me first make sure I have the equations right:

$$\frac{dx}{dt}=x\left[a+b(z-1)\right]+y+cx(x^2+y^2)$$

$$\frac{dy}{dt}=-x$$

$$\frac{dz}{dt}=d(x^2+y^2-z)$$

with a,b,c,d parameters.

My first thought is have you drawn a Feigenbaum plot? What are the values of the parameters which yield a chaotic trajectory and what are some values which yield a limit cycle? Also what is the basin of attracton: For others, that's the local region around the attractor if the initial values of x,y,and z are started in the basin, the trajectory is drawn to the attractor. Initial values of x,y z outside of the basin result in a trajectory which flies off to infinity.

Is this your equation or a study of one which has documentation elsewhere?

 

[Tom Piper]

>My first thought is have you drawn a Feigenbaum plot?

No, I haven't done Feigenbaum plot for that dynamical system.

>What are the values of the parameters which yield a chaotic trajectory and what >are some values which yield a limit cycle?

As shown in my website, when setting the parameters in your notations as,

a = 0.2
b = 4
c = 0.1
d = 0.1

and setting the initial values as x = 2, y = 2, and z = 2, then the trajectory shows chaotic behavior. If c = 0, then that trajectory converges to a simple limit cycle. However if c = 0.15, then the trajectory converges to a strange "limit cycle" as indicated in my earlier post.

>Also what is the basin of attracton?

That's the problem. The region: {(x,y,z); x^2+y^2-z < 0} is a part of the basin, but obviously there are other (unknown) regions where the trajectory falls down to the attractor. As indicated in the last part of my webpage, it is necessary to do further search for them.

>Is this your equation or a study of one which has documentation elsewhere?

No, it's my original equation :)

 

[Tom Piper]

Sorry, for the values of a and b shown above are incorrect. The right values are;

a = -0.2
b = -4

 

[saltydog]

Hello Tom. Got Mathematica? The code below generated the attached plot. Starting to look like a Feigenbaum plot. The parameter range is for d ranging from 0.093 to 0.1. The second plot gives a better idea what the first plot is doing. The second plot is a cross-section of the attractor at the z-axis for d=0.1. The first plot then just plots the y-intercept of this plot as d ranges through the values indicated.

index = 0;
a = -0.2;
b = -4;
c = 0.1;
psave = Table[{0, 0}, {200000}];
For[d = 0.093, d <= 0.1, d += 0.00001,
    sol1 = 
      NDSolve[{x'[t] == 
            x[t](a + b(z[t] - 1)) + y[t] + c x[t](x[t]^2 + y[t]^2), 
          y'[t] == -x[t], z'[t] == d(x[t]^2 + y[t]^2 - z[t]), x[0] == 2, 
          y[0] == 2, z[0] == 2}, {x, y, z}, {t, 0, 1000}, AccuracyGoal -> 10, 
        MaxSteps -> 50000];
    fnx[t_] := Evaluate[x[t] /. Flatten[sol1]];
    fny[t_] := Evaluate[y[t] /. Flatten[sol1]];
    fnz[t_] := Evaluate[z[t] /. Flatten[sol1]];
    xsign = Sign[fn1[200]];
    For[t2 = 200, t2 <= 1000, t2 += .1,
      If[Sign[fnx[t2]] != xsign,
          index++;
          psave[[index]] = {d, fny[t2]};
          xsign = Sign[fnx[t2]];
          ];
      ];
    ];
index
psave = Take[psave, index];
ListPlot[psave, Prolog -> AbsolutePointSize[1]];

 

[Tom Piper]

>Hello Tom. Got Mathematica?

Unfortunately I don't have it. It seems too expensive to buy for me... Anyway thank you for your analysis using Mathematica. The result of Feigenbaum plot is very interesting. Did you notice the blank around d=0.0985 in that plot? I imagined that blank might imply possible another strange limit cycle and tried to do the simulation around d=0.0985. Just right! I found a strange limit cycle when d=0.0982 shown in the following picture and updated my webpage referring to your contribution.

 

[saltydog]

Hello Tom. I've been working on the Feigenbaum plot of this system for d in the range of 0.093 to 0.1. I'm unable to obtain a better-resolved plot of the bifurcations, the limit cycles, and the chaotic regime. Here's my main problem: The following Mathematica code runs the numeric ODE solver for your system:

d = 0.09831;
sol1 = NDSolve[{x'[t] == 
          x[t](-0.2 - 4(z[t] - 1)) + y[t] + 0.1 x[t](x[t]^2 + y[t]^2), 
        y'[t] == -x[t], z'[t] == d(x[t]^2 + y[t]^2 - z[t]), x[0] == 1, 
        y[0] == 1, z[0] == 1}, {x, y, z}, {t, 0, 2000}, MaxSteps -> 275000];
fnx[t_] := Evaluate[x[t] /. Flatten[sol1]];
fny[t_] := Evaluate[y[t] /. Flatten[sol1]];
pp1 = ParametricPlot[Evaluate[{x[t], y[t]} /. sol1], {t, 1500, 2000}, 
    PlotPoints -> 2000, PlotRange -> {{-5, 5}, {-5, 5}}]

I know it's cryptic. It's just plotting x(t) and y(t) parametrically. This is the first plot below. Note where the curve crosses the y-axis. It is these points, plotted as a function of d which makes up the Feigenbaum plot above. If you were to look at the Feigenbaum plot at d=0.09831, you'd see three points above the center line, some right at the center, and then 3 more below the center line reflecting the points where the first plot crosses the y-axis.

Now, when I change d to 0.09832, I get the second plot. This is NOT correct and reflects some abberation of round-off error or some other problem with Mathematica as the trajectory veers off the limit cycle. Note in the above NDSolve code, I let t go from 0 to 2000. If I change it to 3000, I get the third plot which IS the correct limit cycle. However, leaving it at 3000 soon causes problems with other points. I believe there is a beautiful Feigenbaum plot underneath all this spurious behavior. I'm not sure how to proceed further. It requires more study. Just for the record, I'll post the code to generate the Feigenbaum Plot which is a bit messy:

plotlist = Table[{0, 0}, {100000}];
index = 0;
For[d=0.093,d<=0.1,d+=0.000001,
 sol1 = 
      NDSolve[{x'[t] == 
            x[t](-0.2 - 4(z[t] - 1)) + y[t] + 0.1 x[t](x[t]^2 + y[t]^2), 
          y'[t] == -x[t], z'[t] == d(x[t]^2 + y[t]^2 - z[t]), x[0] == 1, 
          y[0] == 1, z[0] == 1}, {x, y, z}, {t, 0, 1000}, MaxSteps -> 5000000,
         MaxStepSize -> 0.001];
fnx[t_] := Evaluate[x[t] /. Flatten[sol1]];
fny[t_] := Evaluate[y[t] /. Flatten[sol1]];
    Clear[ftable];
    ftable = Table[{t, fnx[t]}, {t, 500, 1000, 0.01}];
fsize = Length[ftable];
cursign = Sign[ftable[[1, 2]]];
oldpt = ftable[[1, 1]];
For[i = 2, i <= fsize, i++,
      newpt = ftable[[i, 1]];
      If[Sign[ftable[[i, 2]]] != cursign,
        curpt = oldpt + (newpt - oldpt)/2;
        While[Abs[fnx[curpt]] > 0.001,
          If[Sign[fnx[curpt]] != cursign,
            newpt = curpt;
            ,
            oldpt = curpt;
            ];
          curpt = oldpt + (newpt - oldpt)/2;
          ];
        cursign = -cursign;
        index++;
        plotlist[[index]] = {d, fny[curpt]};
        ];
      ];
    ];

plotlist = Take[plotlist, index];
lp3 = ListPlot[plotlist, PlotRange -> {{0.093, 0.1}, {-3, 3}}, 
      Prolog -> AbsolutePointSize[1]];

Here's a better picture, still not good enough, of the Feigenbaum plot:

Link to feigenbaum plot

 

[Tom Piper]

Hello Tom. I've been working on the Feigenbaum plot of this system for d in the range of 0.093 to 0.1. I'm unable to obtain a better-resolved plot of the bifurcations, the limit cycles, and the chaotic regime. Here's my main problem: The following Mathematica code runs the numeric ODE solver for your system:

d = 0.09831;
sol1 = NDSolve[{x'[t] == 
          x[t](-0.2 - 4(z[t] - 1)) + y[t] + 0.1 x[t](x[t]^2 + y[t]^2), 
        y'[t] == -x[t], z'[t] == d(x[t]^2 + y[t]^2 - z[t]), x[0] == 1, 
        y[0] == 1, z[0] == 1}, {x, y, z}, {t, 0, 2000}, MaxSteps -> 275000];
fnx[t_] := Evaluate[x[t] /. Flatten[sol1]];
fny[t_] := Evaluate[y[t] /. Flatten[sol1]];
pp1 = ParametricPlot[Evaluate[{x[t], y[t]} /. sol1], {t, 1500, 2000}, 
    PlotPoints -> 2000, PlotRange -> {{-5, 5}, {-5, 5}}]

I know it's cryptic. It's just plotting x(t) and y(t) parametrically. This is the first plot below. Note where the curve crosses the y-axis. It is these points, plotted as a function of d which makes up the Feigenbaum plot above. If you were to look at the Feigenbaum plot at d=0.09831, you'd see three points above the center line, some right at the center, and then 3 more below the center line reflecting the points where the first plot crosses the y-axis.

Now, when I change d to 0.09832, I get the second plot. This is NOT correct and reflects some abberation of round-off error or some other problem with Mathematica as the trajectory veers off the limit cycle. Note in the above NDSolve code, I let t go from 0 to 2000. If I change it to 3000, I get the third plot which IS the correct limit cycle. However, leaving it at 3000 soon causes problems with other points. I believe there is a beautiful Feigenbaum plot underneath all this spurious behavior. I'm not sure how to proceed further. It requires more study. Just for the record, I'll post the code to generate the Feigenbaum Plot which is a bit messy:

plotlist = Table[{0, 0}, {100000}];
index = 0;
For[d=0.093,d<=0.1,d+=0.000001,
 sol1 = 
      NDSolve[{x'[t] == 
            x[t](-0.2 - 4(z[t] - 1)) + y[t] + 0.1 x[t](x[t]^2 + y[t]^2), 
          y'[t] == -x[t], z'[t] == d(x[t]^2 + y[t]^2 - z[t]), x[0] == 1, 
          y[0] == 1, z[0] == 1}, {x, y, z}, {t, 0, 1000}, MaxSteps -> 5000000,
         MaxStepSize -> 0.001];
fnx[t_] := Evaluate[x[t] /. Flatten[sol1]];
fny[t_] := Evaluate[y[t] /. Flatten[sol1]];
    Clear[ftable];
    ftable = Table[{t, fnx[t]}, {t, 500, 1000, 0.01}];
fsize = Length[ftable];
cursign = Sign[ftable[[1, 2]]];
oldpt = ftable[[1, 1]];
For[i = 2, i <= fsize, i++,
      newpt = ftable[[i, 1]];
      If[Sign[ftable[[i, 2]]] != cursign,
        curpt = oldpt + (newpt - oldpt)/2;
        While[Abs[fnx[curpt]] > 0.001,
          If[Sign[fnx[curpt]] != cursign,
            newpt = curpt;
            ,
            oldpt = curpt;
            ];
          curpt = oldpt + (newpt - oldpt)/2;
          ];
        cursign = -cursign;
        index++;
        plotlist[[index]] = {d, fny[curpt]};
        ];
      ];
    ];

plotlist = Take[plotlist, index];
lp3 = ListPlot[plotlist, PlotRange -> {{0.093, 0.1}, {-3, 3}}, 
      Prolog -> AbsolutePointSize[1]];

Here's a better picture, still not good enough, of the Feigenbaum plot:

Link to feigenbaum plot

Sorry for my response has been delaying. Oh, that's a nice drawing for that limit cycle! By the way I've been continuing to find any other limit cycle in this dynamical system. And I found it at last! Strange enough, the parameters for that limit cycle are exactly same as the previous ones, i.e.

a = 0.2
b = 4
c = 0.1
d = 0.1

Using those parameters, I got the limit cycle shown in the right side of the attached file after 6 million iteration. As the reference, I put in the left side the trajectory after 4.2 million iteration. Namely that limit cycle is a very weak attractor. (Exactly saying, that limit cycle should be called a tube because it has a finite sectional area, although it is very small.)

I updated my webpage to explain those things as,
http://www.geocities.com/tontokohirorin/mathematics/limitcycle/limitcycle2.htm
Also I showed in that webpage the Feigenbaum Plot you've obtained. If you feel it is inappropriate, please let me know. I'll withdraw it from the webpage immediately.

 

[saltydog]

Tom, that's nice of you to refer to me. I have a suggestion for you. Remember how I told you that some type of error was entering into the calculations for the feigenbaum plots preventing me from plotting what I think is a very beautiful filigree feigenbaum plot? The error as I recall seemed to look like the same error I obtained when attempting to use a non-symplectic algorithm for a Hamiltonian system. If this is of interest to you, refer to the DE forum under "Symplectic Algorithm", ignore the little mistake I made at the start of the thread and look at the two plots for the non-linear pendulum. That's very similar to the errors I was getting when I generated the feigenbaum plot for your system. Perhaps if a symplectic algorithm was used, the feigenbaum diagram would be more resolved. Check it out.

 

[Tom Piper]

Tom, that's nice of you to refer to me. I have a suggestion for you. Remember how I told you that some type of error was entering into the calculations for the feigenbaum plots preventing me from plotting what I think is a very beautiful filigree feigenbaum plot? The error as I recall seemed to look like the same error I obtained when attempting to use a non-symplectic algorithm for a Hamiltonian system. If this is of interest to you, refer to the DE forum under "Symplectic Algorithm", ignore the little mistake I made at the start of the thread and look at the two plots for the non-linear pendulum. That's very similar to the errors I was getting when I generated the feigenbaum plot for your system. Perhaps if a symplectic algorithm was used, the feigenbaum diagram would be more resolved. Check it out.

Thank you for suggestion. I'll visit the forum you've indicated later. I think to estimate the influence of the error in the numerical analysis is difficult particularly for chaotic system since the effect of any small error could become tremendously great after repeating many iterations. By the way I appreciate for your comment to another thread I posted regarding logistic equation. The method employed in the webpage linked in that thread is to represent the logistic equation in a quotient space. Similar method has been applied to analyze the chaotic behavior shown on the dynamical system described above.