- #1

- 25

- 0

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!

You are using an out of date browser. It may not display this or other websites correctly.

You should upgrade or use an alternative browser.

You should upgrade or use an alternative browser.

- Thread starter Tom Piper
- Start date

- #1

- 25

- 0

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!

- #2

saltydog

Science Advisor

Homework Helper

- 1,582

- 3

Tom Piper said:

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.

- #3

- 25

- 0

- #4

saltydog

Science Advisor

Homework Helper

- 1,582

- 3

Let me first make sure I have the equations right:

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

[tex]\frac{dy}{dt}=-x[/tex]

[tex]\frac{dz}{dt}=d(x^2+y^2-z)[/tex]

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?

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

[tex]\frac{dy}{dt}=-x[/tex]

[tex]\frac{dz}{dt}=d(x^2+y^2-z)[/tex]

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?

Last edited:

- #5

- 25

- 0

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 :)

- #6

- 25

- 0

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

a = -0.2

b = -4

a = -0.2

b = -4

- #7

saltydog

Science Advisor

Homework Helper

- 1,582

- 3

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.

Code:

```
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]];
```

Last edited:

- #8

- 25

- 0

>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.

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.

- #9

saltydog

Science Advisor

Homework Helper

- 1,582

- 3

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:

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:

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

Link to feigenbaum plot

Code:

```
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:

Code:

```
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

Last edited:

- #10

- 25

- 0

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.saltydog said: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:

Code:`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:

Code:`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

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.

- #11

saltydog

Science Advisor

Homework Helper

- 1,582

- 3

- #12

- 25

- 0

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.saltydog said:

Share: