Dear @AmusingYeti 

Thank you for the link, which seems more focus on the bifurcation diagram of the one-dimensional map, whereas my model is continuous odes... 

I have uploaded a pure text version of the same code 23kB. The pointploted figure also makes my system extremely slow.

Thank you.

Very kind wishes,

Wang Zhe

Hi @AmusingYeti 

Just attached, thank you for reminding me.

Very kind wishes,

Wang Zhe

Dear @tomleslie 

Thank you so very much for your help all this while! 

Very kind wishes,

Wang Zhe

Dear @tomleslie,

Via studying your code, I plotted the u1(t) against t and found that the resulting figures look very different. Could you please comment on the following? plot_vs_odeplot.mw

Thank you.

Very kind wishes,

Wang Zhe

Dear @vv,

Thank you for the reply! 

Representing say that the probability of Phi to be in the region [0,0.1] is 0.15, in the region [0.1.0.2] is 0.05, and so on. Since I found the value of Phi tends to be near -1, 0, and 1, I am thinking that perhaps I can study the dynamics of the system by looking at the probability distribution of Phi along [-1,1], starting with a list of initial conditions.

Very kind wishes,

Wang Zhe 

Dear @tomleslie,

Thank you for the reply and the code!!! It takes an almost triple amount of time to run on my machine...(:

The function Phi satisfies -1<Phi<1, and it seems that the value of Phi tends to accumulate at the values -1, 0, and 1. I am thinking to start with a list of initial conditions and organize the resulting values of Phi into a single array. Then I would like to plot the probability distribution of Phi along the range [-1,1] and calculate the flatness or skewness of the distribution.

Since the {sys} has reflectional symmetry: f(z)=f(-z), the expected PDF for Phi is symmetric to the zero, see attached.

Thank you.

Very kind wishes,

Wang Zhe

Dear @MapleSystems 

Thank you so very much! This is exactly what I want.

Have a nice weekend!

Very kind wishes,

Wang Zhe

Dear @Markiyan Hirnyk 

Thank you.

Very kind wishes,

Wang Zhe

Dear @Markiyan Hirnyk 

Thank you for your reply!

It looks much better now. But what I do not understand is why curves in 3d look much better than those in a 2d plot? 

Sometimes, even the curves in 3d look less satisfactory, I would like to know if there is a way to improve it further?

restart;
with(DEtools): with(plots): with(plottools):
a := -1: b := -3: c := 3: d := 1: omega := 1: v1 := 1: f := -4: epsilon := 0.1e-1:
sys := diff(u1(t), t) = v1*u1(t)-(omega+u2(t)^2)*u2(t)+u1(t)*(a*(u1(t)^2+u2(t)^2)+b*z(t)^2), 
       diff(u2(t), t) = (omega+u1(t)^2)*u1(t)+v1*u2(t)+u2(t)*(a*(u1(t)^2+u2(t)^2)+b*z(t)^2),
       diff(z(t), t) = z(t)*(-v1+c*(u1(t)^2+u2(t)^2)+z(t)^2)+epsilon*z(t)*(v2+f*z(t)^4): 

solC := dsolve({eval(sys, v2 = 2.0500014987), 
                u1(0)=.6, u2(0)=.6, z(0)=.1}, 
                type = numeric, method = rkf45, maxfun = 0, 
                range = 450 .. 1050): 

p1 := odeplot(solC, 
             [t, sqrt(u1(t)^2+u2(t)^2), z(t)], 
              refine = 5, 
              color = burgundy, 
              thickness = 1,
              orientation=[0,90,0]);
p2 := odeplot(solC, 
             [u1(t),u2(t), z(t)], 
              refine = 5, 
              color = burgundy, 
              thickness = 1,
              orientation=[41,75,0]); 
 

Thank you!

Very kind wishes,

Wang Zhe

Dear @tomleslie

Yes, is is absolutely amazing to look at, thank you for sharing!!! You have found, "accidentally", the very essence of the chaos theory: stretching and folding! 

What we see is how the Duffing attractor is formed in the phase space. For each period, trajectories are firstly stretched and elongated, and then fold onto itself, resulting in the formation of an invariant set in the phase space, the Duffing attractor.  

A similar phenomenon can be found in the formation of many other attractors, representing say, the Lorenz attractor. Although I failed to visualize it myself, I found a fantastic animation on YouTube, started at 7:40.

https://www.youtube.com/watch?v=aAJkLh76QnM&index=7&list=PLelIK3uylPMHTEZ0hEx3PshdSx6awKmxa&t=7s

If you have any idea how this can be done for the Lorenz system, could you please let me know?

Thank you. 

Very kind wishes,

Wang Zhe

Dear @tomleslie

Fantastic! By following your approach I have visualized the period doubling phenomenon in the Lorenz system.

Yes, I admit that it does not make much sense by investing the intersection points on the line u(t)=v(t) in the Duffing system, except for sigma=0. As you have already indicated that the addition of the term -sigma*v(t) directs the vector field towards one of the fixed points on the x-axis so that it is unlikely to have any intersection. Nevertheless, the result (no intersection) suggests that there is no periodic motion in the Duffing system for Gamma=0 and sigma=0.25.

For a two-dimensional dynamical system, the line u(t)=v(t) actually defines a Poincare section such that the intersection point represents the periodic orbit. Since the Poincare map is orientation preservative, we look at only the positive half of the line. One intersection point indicates a period-one orbit, two intersections indicate a period-two orbit, and so on. By decreasing r from 230 in the Lorenz system, we can observe a period doubling cascade route to chaos, see attached.

So the next thing we can do for the Lorenz system is to plot x(t) against the system parameter r, and the resulting would be the bifurcation diagram!

I cannot thank you enough for your help!!!

Very kind wishes,

Wang Zhe

Duffing_Poincare_return_section_u(t)=v(t).mw

Lorenz_Poincare_return_section_x(t)=y(t).mw 

Dear @tomleslie

You are absolutely awesome! It is a real pleasure to read your code and comments! I did not that it could be done in this way. Thank you!

I have not ignored your code for the autonomous system. However, in your code, you discretized the solution trajectory with respect to time, and the resulting solution points are contained in the two-dimensional (u(t),v(t)) plane; whereas I want to plot the intersection points of the solution trajectory and a specific line: u(t)=v(t), so that the expected results are discrete points on the one-dimensional line u(t)=v(t), see attached, Figure 6.5(c).

Let us back to the autonomous Duffing oscillator by setting Gamma=0, we have

restart:
with(DEtools): with(plots):
sigma:=0.25: omega:=1:
Duffing:=diff(u(t),t)=v(t),
           diff(v(t),t)=u(t)-u(t)^3-sigma*v(t):
solA:=dsolve( {Duffing,u(0)=0.5,v(0)=0.5},
              {u(t),v(t)},
              numeric,
              method=rkf45,
              maxfun=10^7,
              output=Array([seq(i,i=500..1500)])):

With dsolve() we obtained information for [t, u(t), v(t)] and let us focused on [u(t), v(t)]. So, what I want to plot is actually the points that are both contained in the matrix [u(t),v(t)] and on the line, u(t)=v(t).

Moreover, I would like to know if there is a way to connect the discrete points in pointplot() by a smooth solid curve. I have tried the option: connect=true, and increase the sampling intervals to i=0.01, but the curve still looks furry.

Thank you!

Very kind wishes,

Wang Zhe

Period_doubling_in_the_Lorenz_system.pdf

Dear @tomleslie 

Thank you so very much for the attached file, it is much more efficient and elegant than mine. However, when Gamma=0, t=2*Pi*i/omega is no longer the Poincare section. In this case, the Duffing oscillator is two-dimensional and the corresponding Poincare section should be a line, representing say, u=v.

I would like to study the intersections as trajectories crossing the line u=v transversely, but I have totally no idea how to define the seq().

I have specified the problem in the attached file, could you please help me.

Thanks a lot.

Very kind wishes,

Wang Zhe 

Duffing_Poincare_return_map.mw

Dear @tomleslie 

Thank you for the suggestion.

Meanwhile, I would like to investigate the Poincare return section, and I was wondered if you may point me out.

The following is a code for drawing a Poincare return section for the Duffing oscillator, and I would like to know if I can modify the code for the autonomous system.

restart; with(DEtools); with(plots);
Duffing := diff(u(t), t) = v(t), diff(v(t), t) = u(t)-u(t)^3-sigma*v(t)+Gamma*cos(omega*t);
sigma := .25; omega := 1; Gamma := .3;
soln := dsolve({Duffing, u(0) = .5, v(0) = .5}, {u(t), v(t)}, type = numeric, method = rkf45, maxfun = 10^7, output = procedurelist);
pt := array(0 .. 5000); u1 := array(0 .. 5000); v1 := array(0 .. 5000);
imax := 1000;
for i from 0 to imax do u1[i] := eval(u(t), soln(2*Pi*i/omega)); v1[i] := eval(v(t), soln(2*Pi*i/omega)) end do;
pts := [[u1[n], v1[n]]$n = 150 .. imax];
pointplot(pts);

Thanks in advance.

Very kind wishes,

Wang Zhe