(*Date: February 28, 2014.*)
(*Author: Nazmi Burak Budanur*)
(*Study of the two-mode system.*) 
(*Covered ChaosBook exercises (as of ver 14.5.2): 10.7, 10A.6-10 and 10A.12-14*)

(*2-mode velocity function:*)

v[{x1_, y1_ ,x2_, y2_}] := {((Subscript[μ, 1])+(Subscript[a, 1])(x1^2 + y1^2) + (Subscript[b, 1])(x2^2 + y2^2))x1 + (Subscript[c, 1])(x1 x2 + y1 y2),((Subscript[μ, 1])+(Subscript[a, 1])(x1^2 + y1^2) + (Subscript[b, 1])(x2^2 + y2^2))y1 + (Subscript[c, 1])(x1 y2 - x2 y1),((Subscript[μ, 2])+(Subscript[a, 2])(x1^2 + y1^2) + (Subscript[b, 2])(x2^2 + y2^2))x2 + (Subscript[c, 2])(x1^2 - y1^2) + (Subscript[e, 2])y2, ((Subscript[μ, 2])+(Subscript[a, 2])(x1^2 + y1^2) + (Subscript[b, 2])(x2^2 + y2^2))y2 + 2(Subscript[c, 2])x1 y1 - (Subscript[e, 2])x2};

(*Paremeters that are 0 or 1*)par01 = {(Subscript[μ, 2])-> 1, (Subscript[e, 2])-> 1, (Subscript[a, 1])-> -1, (Subscript[b, 1])-> 0, (Subscript[b, 2])-> 0, (Subscript[c, 2])-> 1} ;

pars = {(Subscript[μ, 1])-> -2.8, (Subscript[a, 2])-> -2.66, (Subscript[c, 1])-> -7.75};

vpars[{x1_,y1_,x2_,y2_}] := v[{x1, y1, x2, y2}] /. par01 /. pars

x0={0.439966,0,-0.386267,0.070204};

(*10.7 - Integrating the two-mode system:*)

sol = NDSolve[{{x1sol'[t], y1sol'[t],x2sol'[t], y2sol'[t]  } == vpars[{x1sol[t], y1sol[t], x2sol[t], y2sol[t]}],{x1sol[0], y1sol[0],x2sol[0], y2sol[0] }== x0}, {x1sol,y1sol,x2sol,y2sol}, {t,0,400}, MaxSteps-> ∞, Method-> {"ExplicitRungeKutta", "DifferenceOrder"-> 8}, InterpolationOrder-> All];

(*Plot:*)

ParametricPlot3D[Evaluate[{x1sol[t], y1sol[t], x2sol[t]}/.sol ], {t,0,400}, PlotPoints-> 4000, Boxed-> False, AxesLabel->{(Subscript[x, 1]),(Subscript[y, 1]),(Subscript[x, 2])}]

(*10A.6 - Invariant subspace of the two-mode system: *)
FullSimplify[v[{0,0, x2, y2}]]

(*SO(2) Lie algebra generator:*)

T = {{0,1,0,0},{-1,0,0,0},{0,0,0,2},{0,0,-2,0}};

(*Phase velocity:*)
vϕ [xRed_] := (v[xRed].(T.xp))/((T.xRed).(T.xp))

(*Velocity within the slice:*)
vred[xRed_] := v[xRed] - vϕ[xRed](T.xRed)

(*10A.7a - Phase velocity for the 1st mode slice:*)

slicep = {1,0,0,0};

vϕ1stmode = vϕ[{x1, y1, x2, y2}] /. xp -> slicep //FullSimplify

(*Note that numerator = v[[2]]*)

(*10A.7b, c, d: Slice hyperplane border: x1=0, d - 2 = 2 dimensional, it is the invariant subspace (x1=0, y1=0), symmetry reduced trajectory cannot cross the slice hpyerplane border*)

(*Reduced velocity function:*)

vred1stmode = vred[{x1, y1, x2, y2}] /. {xp-> slicep, y1 -> 0} // FullSimplify

(*10A.10 - Relative equilibria:*)

vredsimp = vred1stmode /. par01

condx2=Solve[vredsimp[[1]]==0, x2]

condy2 = Solve[vredsimp[[4]]==0,y2]

fx1 = vredsimp[[3]]/.condy2[[1]]/.condx2[[1]]//FullSimplify

fx1A = fx1/. x1^2 -> A

Solve[fx1A == 0, A]/.pars

x1reqv = Sqrt[0.1935697113530641`]

0.43996557973671546`

x2reqv = x2/.condx2[[1]]/.pars/.x1->x1reqv

-0.3862670595294276`

y2reqv = y2/.condy2[[1]]/.pars/.{x1->x1reqv, x2->x2reqv}

(*Search for relative equilibria in the invariant subspace:*)

vRedInvssp = vred1stmode /. x1 -> 0

Solve[vRedInvssp=={0,0,0,0}, {x2,y2}]//FullSimplify

(*We are only interested in the real roots, so the last two are unimportant.*) (*First two non-zero roots correspond to the relative equilibrium that we previously found using Hilbert basis. Slice does the job:*)

x2^2 + y2^2 /. {x2->-((Sqrt[-(Subscript[b, 2]) (Subscript[e, 2]^2)-4 (Subscript[c, 1]^2) (Subscript[μ, 2])])/(2 (Sqrt[(Subscript[b, 2])]) (Subscript[c, 1]))),y2->-((Subscript[e, 2])/(2 (Subscript[c, 1])))}//Simplify

(*Search for relative equilibria within the slice:*)

vred1stmode[[1]]

(*We divide out the root x1 = 0 as we know that it corresponds to the invariant subspace:*)

f1 = vred1stmode[[1]]

solx1 = Solve[f1 == 0, x1]//FullSimplify

(*We now got a condition on x1^2 for relative equilibria within the slice, we continue by substituting this into the 3rd and 4th velocity components:.*)

f3 = vred1stmode[[3]] /. solx1[[3]]//FullSimplify

soly2 = Solve[f3==0, y2]//FullSimplify

f4 = vred1stmode[[4]] /. solx1[[3]]//FullSimplify

(*Try to f3 and f4 for r2sq = x2^2 + y2^2*)

f3r2 = f3 /. (x2^2 + y2^2) -> r2sq //FullSimplify;

f4r2 = f4 /. (x2^2 + y2^2) -> r2sq//FullSimplify;

Solve[f3r2 == 0, r2sq]//FullSimplify;

(*Cannot see any other helpful simplification*)

(*Simplify parameters:*)

f3simple= f3 /. par01 //FullSimplify

f4simple = f4 /. par01 // FullSimplify

solx2simple=Solve[f4simple == 0, x2]

fy2 = f3simple /. solx2simple[[1]]//FullSimplify

(*fy2 == 0 corresponds to the solution of a 4th order polynomial, hence the solution is exact.*)

(*However, it's long, so at this point, we plug the numbers in and see what we have:*)

soly2Num = Solve[fy2 == 0/. pars, y2]

(*Meaningful solutions are real ones.*)

solx2Num = solx2simple /. pars /. soly2Num

(*Finally, we should check x1^2 for these numbers and find relative equilibria:*)

solx1simple = solx1/.par01

solx1SimpleNum = solx1simple /. pars

(*Note, in this case, x1 depends only on x2*)

solx1Num =solx1SimpleNum /. solx2Num

(*We have two physical relative equilibrium (Satisfying x1>0 and x1 ∈ Reals):*)

xreqv = Re[({x1, 0, x2, y2} /. solx1Num[[4]][[1]][[3]] /. solx2Num[[4]] /. soly2Num[[4]])[[1]]]

(*Let us verify that this is a relative equilibrium by plugging this into the reduced velocity:*)

vred1stmode /. par01 /. pars

vRedNum[{x1_,y1_,x2_,y2_}] := {x1 (-2.8`-(x1^2)-7.75` x2),0,(x1^2)+x2-2.66` (x1^2) x2+(1-15.5` y2) y2,-x2 (1-15.5` y2)+y2-2.66` (x1^2) y2}

vRedNum[xreqv]

(*Seems correct!*)

(*Solve the symmetry reduced equations:*)

solred = NDSolve[{{x1redsol'[t], y1redsol'[t],x2redsol'[t], y2redsol'[t]} == vRedNum[{x1redsol[t], y1redsol[t], x2redsol[t], y2redsol[t]}],{x1redsol[0], y1redsol[0],x2redsol[0], y2redsol[0] }== x0}, {x1redsol,y1redsol,x2redsol,y2redsol}, {t,0,1000}, MaxSteps-> ∞, Method-> {"ExplicitRungeKutta", "DifferenceOrder"-> 8, "StiffnessTest"->False}, InterpolationOrder-> All ];

(*Plot reduced trajectory*)

redPlot = ParametricPlot3D[Evaluate[{x1redsol[t], x2redsol[t], y2redsol[t]}/.solred ], {t,0,1000}, PlotPoints-> 1000, Boxed-> False, AxesLabel->{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 2])}, BoxRatios->{1,1,1} , PlotRange-> {{0,2.5},{-2.5,1},{-0.5,0.6}}]

(*10A.9 - Visualizing the relative equilibrium:*)

ParametricPlot3D[{Evaluate[{x1redsol[t], x2redsol[t], y2redsol[t]}/.solred ], Evaluate[{x1sol[t], x2sol[t], y2sol[t]}/.sol ]} , {t,0,250}, PlotPoints-> 500, Boxed-> False, AxesLabel->{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 2])}, BoxRatios->{1,1,1} , PlotRange-> Automatic, FaceGrids->{{-1,0,0},{0,0,-1}, {0,1,0}}]

(*Red -> Symmetry equivariant state space, Blue -> Reduced state space trajectory*)

(*10A.10 - Stability of Relative Equilibrium:*)

(*Let us start by computing the stability matrix of the 2-mode system in the reduced statespace:*)

Ared = IdentityMatrix[4] (*Dummy assignment*);

xdummy = {x1,y1,x2,y2};

For[i=1, i≤ 4, i++, For[j = 1, j≤ 4, j++, Ared[[i,j]] = D[ vred[xdummy][[i]], xdummy[[j]]]]];

Ared = Ared /. xp -> slicep /. y1->0 /. par01//FullSimplify

Ared/.{x1-> (Subscript[x, 1]), y1-> (Subscript[y, 1]), x2-> (Subscript[x, 2]), y2-> (Subscript[y, 2])}

(*Note:*)
Ared[[2]]

fAredNum[{xx1_,yy1_,xx2_,yy2_}]:=Ared/.pars /.{x1-> xx1, y1-> yy1, x2-> xx2, y2-> yy2}

xreqv

λ = Eigenvalues[fAredNum[xreqv]]

e = Eigenvectors[fAredNum[xreqv]]

(*Note, Re[(Subscript[λ, 3,4])] > 0, hence the unstable direction:*)

unstable = λ[[2]]*e[[2]] + λ[[3]]*e[[3]]

(*This is a real vector:*)

unstable = Re[unstable]/Norm[unstable]

xreqv3D = {xreqv[[1]], xreqv[[3]],xreqv[[4]]}

unstable3D = {unstable[[1]], unstable[[3]], unstable[[4]]};

sum3D = xreqv3D + unstable3D/1.5;

redPlot2 = ParametricPlot3D[Evaluate[{x1redsol[t], x2redsol[t], y2redsol[t]}/.solred ], {t,0,1000}, PlotPoints-> 400, Boxed-> False, AxesLabel->{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 2])}, BoxRatios->{1,1,1} , PlotRange-> {{0,0.80},{-2,0.5},{-0.5,0.4}} ];

Show[redPlot2, Graphics3D[{Red,Arrow[{xreqv3D,sum3D}]}]]

(*Reduced flow, and the unstable direction of the relative equilibrium:*)

(*Now let us try to construct a Poincare section that includes this unstable direction vector:*)

(*We still have the freedom in the orientation of the Poincare section.*)

(*We can choose a normal vector by needing the Poincare section include also the basis in (Subscript[y, 2]) direction:*)

n3D = Cross[{0,0,1}, unstable3D]

n = {n3D[[1]], 0, n3D[[2]], n3D[[3]]} (*Note: this is a particular choice.*)

(*Now, the condition for being on a Poincare section:*)

Poincare[x_] := (x - xreqv).n

(*Interpolated function corresponding to the solution in the reduced state space:*)

fSolred[t_] := Evaluate[{x1redsol[t],y1redsol[t], x2redsol[t], y2redsol[t]}/.solred ][[1]]

(*Find time points of Poincare section intersections:*)

t00 = 0;

tf = 1000;

i=1;

dt = 0.1;

Pt = {} (*Empty array*);

(*Compute Poincare section condition for each time point of the solution, if there is an intersection, then look for a numerically precise root:*)

For[t0=t00, t0<tf-dt, t0=t0+dt, If[Poincare[fSolred[t0]]*Poincare[fSolred[t0+dt]]<0, Pt =Append[Pt, t/. Quiet[FindRoot[Poincare[fSolred[t]],{t,t0}]]]; ]]

(*Discard what seems to correspond to the same intersection:*)

Pt = Union[Pt, SameTest->(Abs[#1 - #2]< 10^-3 &)];

len = Length[Pt];

(*Discard the intersections where the flow velocity and the Poincare section normal points at the opposite directions:*)

For[i=1,i≤len, i++, If[ Sign[n.vRedNum[fSolred[Pt[[i]]]]]==-1, {Pt=Delete[Pt,i]; len=len-1; i--} ]]

Psect = ConstantArray[0, {Length[Pt], 4}] (*Dummy matrix to hold the poincare section intersections*);

Psect3D = ConstantArray[0, {Length[Pt], 3}] (*Dummy matrix to hold the poincare section intersections (3D)*);

PsectRel = ConstantArray[0, {Length[Pt], 3}] (*Dummy matrix to hold the relative poincare section coordinates (3D)*);

For[i = 1, i<= Length[Pt], i++, {Psect[[i]] = fSolred[Pt[[i]]], Psect3D[[i]] ={fSolred[Pt[[i]]][[1]], fSolred[Pt[[i]]][[3]], fSolred[Pt[[i]]][[4]]}, PsectRel[[i]] = Psect3D[[i]] - xreqv3D}];

PsectPlot=ListPointPlot3D[Psect3D, PlotStyle->Directive[PointSize[Medium], Pink]]

a1=ContourPlot3D[Poincare[{x,0,y,z}]==0,{x,xreqv[[1]],0.8},{y,xreqv[[3]],1},{z,-0.2,0.3},AxesLabel->{x,y,z},Mesh->None,ContourStyle->{Directive[Green],Opacity[1]}];

Show[redPlot2, a1, PsectPlot, Graphics3D[{Red,Arrow[{xreqv3D,sum3D}]}]]

(*Now let's project to the basis on the Poincare section. It's done by the following transformation matrix:*)

U = {unstable3D,{0,0,1} , n3D}

(*Relative coordinates of the Poincare section:*)

PsectProj = U.Transpose[PsectRel];

Psect2D = Transpose[PsectProj[[Range[2]]]];

plotPoincare2D = ListPlot[Psect2D, PlotStyle->(Pink)]

(*Now let us interpolate to this data:*)

fPoincare = Interpolation[Psect2D];

plotfPoincare = Plot[fPoincare[x], {x, 0, 0.6}];

Show[plotPoincare2D, plotfPoincare]

(*Function to compute arclengths along the curve:*)

ArcLength[xP_]:= NIntegrate[Sqrt[1+fPoincare'[x]^2] ,{x, 0,xP}]

S = ConstantArray[0,Length[Pt]];

For[i=1, i≤Length[Pt], i++ , S[[i]] = ArcLength[Psect2D[[i]][[1]]]]

(*Produce the return map data:*)

RetMap =Transpose[{S[[1;;Length[S]-1]],S[[2;;Length[S]]]}];

plotRetMap = ListPlot[RetMap, PlotStyle->(Pink), AxesLabel->{(Subscript[s, 𝔫]), (Subscript[s, 𝔫+1])}, AspectRatio->1]

fRetMap = Interpolation[RetMap];

plotfRetMap = Plot[fRetMap[x], {x, 0, 0.6}];

fFixed[x_]:=x;

plotfFixed = Plot[fFixed[x], {x,0,0.6}, PlotStyle->{Black}];

Show[plotRetMap, plotfRetMap,plotfFixed]

(*There seems to be an intersection around (Subscript[s, n]) = 0.42, let's find it precisely:*)

sFixed = s /. FindRoot[fRetMap[s]==fFixed[s], {s,0.42}]

(*Now let's find the corresponding Poincare section coordinate:*)

(*We pick initial guess as the arclength since the vertical components on the 2D Poincare section have quiet small values, the guess is good enough*)

xPS2DFixed = x0 /. FindRoot[ArcLength[x0] == sFixed, {x0, sFixed}]

yPS2DFixed = fPoincare[xPS2DFixed]

(*We now, have our relative periodic orbit candidate in the Poincare section, final step is transforming back to the 4D statespace:*)

U

FixedPoint3D  = Inverse[U].{xPS2DFixed, yPS2DFixed, 0}+xreqv3D

FixedPointssp = {FixedPoint3D[[1]],0, FixedPoint3D[[2]], FixedPoint3D[[3]]}

Poincare[FixedPointssp]

(*To guess the period, we look for the arclength data that is closest to the fixed point arclength we found*)
imin= Position[Abs[S-sFixed],Min[Abs[S-sFixed]]][[1]][[1]]

T0 = Pt[[imin+1]]-Pt[[imin]]

(*Now let us integrate this point:*)

solCycle = NDSolve[{{x1cyclesol'[t], y1cyclesol'[t],x2cyclesol'[t], y2cyclesol'[t]} == vRedNum[{x1cyclesol[t], y1cyclesol[t], x2cyclesol[t], y2cyclesol[t]}],{x1cyclesol[0], y1cyclesol[0],x2cyclesol[0], y2cyclesol[0] }== FixedPointssp}, {x1cyclesol,y1cyclesol,x2cyclesol,y2cyclesol}, {t,0,T0}, MaxSteps-> ∞, Method-> {"ExplicitRungeKutta", "DifferenceOrder"-> 8, "StiffnessTest"->False}, InterpolationOrder-> All ];

fSolcycle[t_] := Evaluate[{x1cyclesol[t],y1cyclesol[t], x2cyclesol[t], y2cyclesol[t]}/.solCycle ][[1]]

cyclePlot = ParametricPlot3D[Evaluate[{x1cyclesol[t], x2cyclesol[t], y2cyclesol[t]}/.solCycle ], {t,0,3.7}, PlotPoints-> 100, Boxed-> False, AxesLabel->{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 2])}, BoxRatios->{1,1,1} , PlotRange-> {{0,2.5},{-2.5,1},{-0.5,0.6}}]

(*10A.12 - Precise relative periodic orbits:*)

rpo={{0.4525719,0.0,0.0509257,0.0335428},{0.4517771,0.0,0.0202026,0.0405222},{0.4514665,0.0,0.0108291,0.0424373},{0.4503967,0.0,-0.0170958,0.0476009}};

Trpo = {3.6415120,7.3459412, 14.6795175, 18.3874094};

g = {1,2,3,4};
nperiods =3;

For[i = 1, i≤4, i++,
Tn  =Trpo[[i]]*nperiods;
solRPOred = NDSolve[{{x1RPOsolred'[t], y1RPOsolred'[t],x2RPOsolred'[t], y2RPOsolred'[t]} == vRedNum[{x1RPOsolred[t], y1RPOsolred[t], x2RPOsolred[t], y2RPOsolred[t]}],{x1RPOsolred[0], y1RPOsolred[0],x2RPOsolred[0], y2RPOsolred[0] }== rpo[[i]]}, {x1RPOsolred,y1RPOsolred,x2RPOsolred,y2RPOsolred}, {t,0,Tn}, MaxSteps-> ∞, Method-> {"ExplicitRungeKutta", "DifferenceOrder"-> 8, "StiffnessTest"->False}, InterpolationOrder-> All ];

solRPO = NDSolve[{{x1RPOsol'[t], y1RPOsol'[t],x2RPOsol'[t], y2RPOsol'[t] } == vpars[{x1RPOsol[t], y1RPOsol[t], x2RPOsol[t], y2RPOsol[t]}],{x1RPOsol[0], y1RPOsol[0],x2RPOsol[0], y2RPOsol[0] }== rpo[[i]]}, {x1RPOsol,y1RPOsol,x2RPOsol,y2RPOsol}, {t,0,Tn}, MaxSteps-> ∞, Method-> {"ExplicitRungeKutta", "DifferenceOrder"-> 8}, InterpolationOrder-> All];

PlotRPO= ParametricPlot3D[Evaluate[{x1RPOsol[t], x2RPOsol[t], y2RPOsol[t]}/.solRPO], {t,0,Tn}, PlotPoints-> 100, Boxed-> False, AxesLabel->{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 2])}, BoxRatios->{1,1,1} , PlotRange-> {{-2,2},{-2.5,2},{-2,1.8}}];

PlotRPORed= ParametricPlot3D[Evaluate[{x1RPOsolred[t], x2RPOsolred[t], y2RPOsolred[t]}/.solRPOred], {t,0,Tn}, PlotPoints-> 100, Boxed-> False, AxesLabel->{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[y, 2])}, BoxRatios->{1,1,1} , PlotRange-> {{0,2.5},{-2.5,1},{-0.5,0.6}}, PlotStyle->{Red}];

g[[i]] = Show[ PlotRPO, PlotRPORed];
]

GraphicsGrid[{{g[[1]], g[[2]]}, {g[[3]], g[[4]]}}, ImageSize-> 600,Spacings-> 0]