QS with an infinite elliptic-saddle and p+q+r

The normal form is given by

change=(((changeSys10/.{ch->0,cd->1,cm->1,cf->df-1})/.df->cf)/.{cf->((ce+df)/dc),cc->dc})/.{df->cf,dc->cc};comp

eq1/.change

eq2/.change

We perform the study of the invariants in terms of the parameter cm and, at the end, we take it equal to one.

Some needed functions

Var={ce,cf,cm,cc} (* We put the last variable for the slices and the penultimate variable for the parameter we shall consider equal to one *)

f[x_]:={ce==x[[1]],cf==x[[2]],cm==x[[3]],cc==x[[4]]}

Thread[Var=={1,2,3,4}]

Thread[Var->{1,2,3,4}]

Fsol[x_,y_]:={ce,cf,cm,cc}/.Flatten[Solve[Union[x,y],{ce,cf,cm,cc}]]

To0[x_]:=x==0

Thread[Var==0]

Map[To0,Var]

ListaBif={};NListaBif={{}}

Lista={{},{},{}}

EstudioSing[varpar_,var_]:={listaeq={};Do[listaeq=Union[listaeq,{ Numerator[Factor[varpar[[j]]-var[[j]]]]==0    }],{j,1,Length[var]}];If[listaeq[[1]]==True,listaeq=Drop[listaeq,1]];Print[listaeq];Do[Do[mmm[j1,j2]=D[listaeq[[j1]][[1]],var[[j2]]   ],{j2,1,Length[var]}],{j1,1,Length[listaeq]}];Clear[M];M=Array[mmm,{Length[listaeq],Length[var]}];Print[MatrixForm[M]];Solve[Union[Map[To0,Flatten[Minors[M,Length[listaeq]]]],listaeq]]}[[1]]

CreaSistema[l1_,l2_,var_]:={listaeq={};Do[listaeq=Union[listaeq,{Numerator[Factor[var[[i]]-l1[[i]]]]}],{i,1,Length[var]}];Do[listaeq=Union[listaeq,{Numerator[Factor[var[[i]]-l2[[i]]]]}],{i,1,Length[var]}];Solve[Map[To0,listaeq]]}

CreaSistema0[l1_,l2_,var_]:={listaeq={};Do[listaeq=Union[listaeq,{Numerator[Factor[var[[i]]-l1[[i]]]]}],{i,1,Length[var]}];Do[listaeq=Union[listaeq,{Numerator[Factor[var[[i]]-l2[[i]]]]}],{i,1,Length[var]}];listaeq=listaeq/.Var[[Length[Var]-1]]->0;Solve[Map[To0,listaeq]]/.Var[[Length[Var]-1]]->0}

CreaSistema2[l1_,l2_,var_]:={listaeq={};Do[listaeq=Union[listaeq,{Numerator[Factor[var[[i]]-l1[[i]]]]}],{i,1,Length[var]}];Do[listaeq=Union[listaeq,{Numerator[Factor[var[[i]]-l2[[i]]]]}],{i,1,Length[var]}];SolveAlways[Map[To0,listaeq],Var[[Length[Var]-1]]]}

Disgregate[x_]:={Do[sol3=Union[sol3,{x[[i2]]}],{i2,1,Length[x]}];sol3}

LocateSqrt[x_]:=Flatten[{sol3={};sol4=Flatten[Disgregate[x]];Do[Do[Flatten[Disgregate[sol4[[i]]]],{i,1,Length[sol4]}];If[sol4==sol3,Break,Print[sol3]];sol4=sol3,{i2,1,10000}];sol5=Cases[sol4,_Power]}]

Bifurcation algebraic sets

change=(((changeSys10/.{ch->0,cd->1,cf->df-1})/.df->cf)/.{cf->((ce+df)/dc),cc->dc})/.{df->cf,dc->cc};comp

eq1/.change

eq2/.change

COMPUTEWS1

When a singularity goes to infinity:

μ0

μ1//FullSimplify

S1=cc cm;

ListaBif=Union[ListaBif,{%}];NListaBif=Union[NListaBif,{"μ1"}]

S1 has two components

solS11=Flatten[Solve[S1[[1]]==0]]

O2x1=Var/.%

Lista[[3]]=Union[Lista[[3]],{%}]

solS12=Flatten[Solve[S1[[2]]==0]]

O2x2=Var/.%

Lista[[3]]=Union[Lista[[3]],{%}]

When we have a coalescence of finite singularities:

𝔻//FullSimplify

We observe that we have to consider only the first and the third components, since the second one appeared before.

S2=(cc-cf) cm (cf+cc (-1+2 cm));

ListaBif=Union[ListaBif,{%}];NListaBif=Union[NListaBif,{"𝔻"}]

solS21=Flatten[Solve[S2[[1]]==0]]

O2x3=Var/.%

Lista[[3]]=Union[Lista[[3]],{%}]

solS23=Flatten[Solve[S2[[3]]==0]]

O2x4=Var/.%

Lista[[3]]=Union[Lista[[3]],{%}]

We verify the presence of finite weak singularities

{𝔗4,𝔗3,𝔗2,𝔗1,σ,𝔅1}//FullSimplify

𝔅1//FullSimplify

S3=((((-1+cc) cc+ce+cf) (cc+(cc^2)-ce-cf-2 cc cm) (cc cf+(cc^2) (-1+cm)+ce cm))/cc);

ListaBif=Union[ListaBif,{%}];NListaBif=Union[NListaBif,{"𝔅1"}]

solS3=Solve[S3==0]//FullSimplify

We observe that S3 has three components

O2x5=Var/.solS3[[1]]

Lista[[3]]=Union[Lista[[3]],{%}];

O2x6=Var/.solS3[[2]]

Lista[[3]]=Union[Lista[[3]],{%}];

O2x7=Var/.solS3[[3]]

Lista[[3]]=Union[Lista[[3]],{%}];

We search for possible existence of invariant straight lines

WB1//FullSimplify

S4=ce cm ((cc^2)-cc (cf-2 cm)-2 (ce+cf) cm) (cc (-cc+cf)+2 ((-1+cc) cc+ce+cf) cm+4 cc (cm^2));

ListaBif=Union[ListaBif,{%}];NListaBif=Union[NListaBif,{"WB1"}]

Note that S4 has 4 components and the second one appeared before

solS41=Flatten[Solve[S4[[1]]==0]]

O2x8=Var/.%

Lista[[3]]=Union[Lista[[3]],{%}];

solS43=Solve[S4[[3]]==0]

O2x9=Var/.solS43[[1]]

Lista[[3]]=Union[Lista[[3]],{%}];

O1x1=Var/.solS43[[2]]

Lista[[2]]=Union[Lista[[2]],{%}];

O1x2=Var/.solS43[[3]]

Lista[[2]]=Union[Lista[[2]],{%}];

solS44=Solve[S4[[4]]==0]//FullSimplify

The second and third solutions appeared before

O2x10=Var/.solS44[[1]]

Lista[[3]]=Union[Lista[[3]],{%}];

Now we are interested in verify when there occurs a coalescence of infinite singularities. We already know that η=0. The next comitant that we must analyse is WM. Calculations yield:

WM//FullSimplify

S5=cc+2cm;

ListaBif=Union[ListaBif,{%}];NListaBif=Union[NListaBif,{"WM"}]

O2x11=Flatten[Var/.Solve[S5==0]]

Lista[[3]]=Union[Lista[[3]],{%}];

We also want to detect the node-focus bifurcation

WW4//FullSimplify

WW7//FullSimplify

Solve[WW7==0]

S6= (2 (cc^3)+(cc^4)+(cc^2) (1+2 ce-2 cf)-2 cc (ce+cf)+(((ce+cf))^2)) cm ((cc^2) (((cc-cf))^2)-2 cc (cc-cf) ((-1+cc) cc+ce+cf) cm+(-4 (cc^3)+(cc^4)+(ce^2)+2 (cc^2) (ce+2 cf)) (cm^2)) ((cc^4)+(((ce+cf))^2)+(cc^2) (-2 ce+2 cf+(((1-2 cm))^2))+2 cc (ce+cf) (-1+2 cm)+(cc^3) (-2+4 cm))

ListaBif=Union[ListaBif,{%}];NListaBif=Union[NListaBif,{"WW7"}]

Solve[S6==0]

This algebraic set has four components, the second one appeared before

solS61=Solve[S6[[1]]==0]//FullSimplify

O2x12=Var/.solS61[[1]]

Lista[[2]]=Union[Lista[[2]],{%}];

O2x13=Var/.solS61[[2]]

Lista[[2]]=Union[Lista[[2]],{%}];

solS63=Solve[S6[[3]]==0,ce]//FullSimplify

solS63disc=Solve[Discrim[S6[[3]],ce]==0,cf]//FullSimplify

sol63indisc=solS63[[1]]/.solS63disc//FullSimplify

O1x4=Var/.sol63indisc[[1]]/.solS63disc[[1]]

Lista[[2]]=Union[Lista[[2]],{%}]

O1x5=Var/.sol63indisc[[2]]/.solS63disc[[2]]

Lista[[2]]=Union[Lista[[2]],{%}]

We point out that if we do Solve[Discrim[S6[[4]],ce]==0,cf] and repeat the last procedure (as we did for S6[[3]]) we will obtain a solution that already appeared before.

Now we analyse the algebraic set due to a presence of an invariant parabola

WPAR=(cc+2 (cc^2)-2 ce-cf) (ce+cf+cc cf)(-cc-2 (cc^2)-2 ce-cf);

S8=WPAR

ListaBif=Union[ListaBif,{%}];NListaBif=Union[NListaBif,{"WPAR"}]

We observe that S8 has three components

O2x14=Var/.Flatten[Solve[S8[[1]]==0]]

Lista[[3]]=Union[Lista[[3]],{%}];

O2x15=Var/.Flatten[Solve[S8[[2]]==0]]

Lista[[3]]=Union[Lista[[3]],{%}];

O2x16=Var/.Flatten[Solve[S8[[3]]==0]]

Lista[[3]]=Union[Lista[[3]],{%}];

Finally we search for the infinite singularity E-H

WN//FullSimplify

O2x17=Var/.{cm->0}

Lista[[3]]=Union[Lista[[3]],{%}];

O2x18=Var/.{cc->-cm}

Lista[[3]]=Union[Lista[[3]],{%}];

Then, we have the following list of regions (given by the invariants).

Lista;

Singularities of the elements from Lista[[3]]

Lista[[3]]

Length[Lista[[3]]]

Now we will determine the singularities of each element from Lista[[3]]

Do[sol={EstudioSing[Lista[[3]][[i]],Var]};If[Length[sol[[1]]]>0,Do[If[Length[sol[[1,i1]]]<Length[Var],Lista[[Length[Var]-Length[sol[[1,i1]]]]]=Union[Lista[[Length[Var]-Length[sol[[1,i1]]]]],{Var/.sol[[1,i1]]}]],{i1,1,Length[sol[[1]]]}]],{i,1,Length[Lista[[3]]]}]

Lista;

Intersection between each two distinct elements from Lista[[3]]

Lista[[3]]

Length[Lista[[3]]]

Now we will determine all the intersection between the elements from Lista[[3]]. In order to do this, in what follows, we take each pair of distinct elements from Lista[[3]],  transform them to an equation (each one equal to zero) and then we solve the respective system of equations:

Do[Do[sol=CreaSistema[Lista[[3]][[i1]],Lista[[3]][[i2]],Var];If[Length[sol[[1]]]>0,Do[If[Length[sol[[1,i]]]<Length[Var],Lista[[Length[Var]-Length[sol[[1,i]]]]]=Union[Lista[[Length[Var]-Length[sol[[1,i]]]]],{Var/.sol[[1,i]]}]],{i,1,Length[sol[[1]]]}]],{i2,i1+1,Length[Lista[[3]]]}],{i1,1,Length[Lista[[3]]]-1}]

Then, we have finished describing the intersection between each pair of distinct elements from Lista[[3]].

Lista;

Tangencies between distinct elements from Lista[[3]] and planes given by c= cte and m=1

Lista[[3]]

Length[Lista[[3]]]

Detecting the existence of square roots and extracting the discriminants

Do[Print[Lista[[3,i]]];sol=Solve[Thread[Var==Lista[[3,i]]]];Do[sol2=sol[[i1,1,1]] Denominator[sol[[i1,1,2]]]-Numerator[sol[[i1,1,2]]];sol3={};Do[sol3=Union[sol3,{sol2[[i2]]}],{i2,1,Length[sol2]}];sol4=LocateSqrt[sol2];Do[If[sol4[[i3,2]]==1/2,sol5=Solve[sol4[[i3,1]]==0];Do[Print[{Var/.Simplify[sol[[i1]]/.sol5[[i4]]]/.sol5[[i4]]}];cont=Length[Var]-Length[Union[Thread[Var==(Var/.Union[Simplify[sol[[i1]]/.sol5[[i4]]],sol5[[i4]]])]]-1];Lista[[cont]]=Union[Lista[[cont]],{Var/.Simplify[sol[[i1]]/.sol5[[i4]]]/.sol5[[i4]]}],{i4,1,Length[sol5]}],Print[0]],{i3,1,Length[sol4]}],{i1,1,Length[sol]}],{i,1,Length[Lista[[3]]]}]

Lista;

Singularities of the elements from Lista[[2]]

Length[Lista[[2]]]

Now we will determine the singularities of each element from Lista[[2]]

In order to obtain the singularities of the curves, we take each element from Lista[[2]], transform it to an equation and then verify the existence (or not) of singularities. In what follows, the index “i” represents the position of each element from Lista[[2]]

Do[sol={EstudioSing[Lista[[2]][[i]],Var]};Print[i,sol];If[Length[sol[[1]]]>0,Do[If[Length[sol[[1,i1]]]<Length[Var],Lista[[Length[Var]-Length[sol[[1,i1]]]]]=Union[Lista[[Length[Var]-Length[sol[[1,i1]]]]],{Var/.sol[[1,i1]]}]],{i1,1,Length[sol[[1]]]}]],{i,1,Length[Lista[[2]]]}]

Lista;

Intersection between each two distinct elements from Lista[[2]]

Length[Lista[[2]]]

Now we will determine all the intersection between the elements from Lista[[2]]. In order to do this, in what follows, we take each pair of distinct elements from Lista[[2]],  transform them to an equation (each one equal to zero) and then we solve the respective system of equations:

Do[Do[sol=CreaSistema[Lista[[2]][[i1]],Lista[[2]][[i2]],Var];If[Length[sol[[1]]]>0,Do[If[Length[sol[[1,i]]]<Length[Var],Lista[[Length[Var]-Length[sol[[1,i]]]]]=Union[Lista[[Length[Var]-Length[sol[[1,i]]]]],{Var/.sol[[1,i]]}]],{i,1,Length[sol[[1]]]}]],{i2,i1+1,Length[Lista[[2]]]}],{i1,1,Length[Lista[[2]]]-1}]

Lista;

Tangencies between distinct elements from Lista[[2]] and planes given by c= cte and m=1

Length[Lista[[2]]]

2) Detecting the existence of square roots and extracting the discriminants

Do[sol=Solve[Thread[Var==Lista[[2,i]]],{Var[[1]]}];Do[sol2=sol[[i1,1,1]] Denominator[sol[[i1,1,2]]]-Numerator[sol[[i1,1,2]]];sol3={};Do[sol3=Union[sol3,{sol2[[i2]]}],{i2,1,Length[sol2]}];sol4=LocateSqrt[sol2];Do[If[sol4[[i3,2]]==1/2,sol5=Solve[sol4[[i3,1]]==0];Do[Print[{Var/.Simplify[sol[[i1]]/.sol5[[i4]]]/.sol5[[i4]]}];Lista[[1]]=Union[Lista[[1]],{Var/.Simplify[sol[[i1]]/.sol5[[i4]]]/.sol5[[i4]]}],{i4,1,Length[sol5]}],Print[0]],{i3,1,Length[sol4]}],{i1,1,Length[sol]}],{i,1,Length[Lista[[2]]]}]

Lista;

List of singular slices

Length[Lista[[1]]]

ListaSlices={}

Do[If[NumberQ[Lista[[1,i,3]]],Print[0],sol=Lista[[1,i,4]]/.Solve[Lista[[1,i,3]]==1];Print[i,sol];If[Element[sol,Reals],ListaSlices=Union[ListaSlices,sol]]],{i,1,Length[Lista[[1]]]}]

ListaSlices

Length[ListaSlices]

Sort[ListaSlices//N]

NListaBif=Drop[NListaBif,-1]

Do[Print[i," ",Lista[[1,i]]," ",Thread[NListaBif->ListaBif/.Thread[Var->Lista[[1,i]]]]],{i,1,Length[Lista[[1]]]}]//FullSimplify

Sort[ListaSlices,Greater]

Length[ListaSlices]

We point out that, since we have shown that the bifurcation diagram has a symmetry, it is important to verify when we have intersections along the curve cf=0, so we must add some new objects (if they exist) in the previous list of slices. In what follows, we consider the list of bifurcation curves and we put cm=1 and cf=0.

bif=ListaBif/.cm->1/.cf->0

Now we calculate the intersection between each pair of distinct elements of the previous list.

Do[sol=Solve[{bif[[i]]==0,bif[[j]]==0}];Print[sol],{i,1,Length[bif]-1},{j,i+1,Length[bif]}]

Observe that the previous computations give us the following list of values of cc.

{0,-2,1/4,-1,1,-9/16,-1/4,1/2,4,2-(Sqrt[3]),2+(Sqrt[3]),-4,-16/9,2};

Sort[{-2,1/4,-1,1,-9/16,-1/4,1/2,4,2-(Sqrt[3]),2+(Sqrt[3]),-4,-16/9,2},Greater]

We define the new list of slices we must consider in our study:

NewListaSlices=Sort[Union[ListaSlices,{0,-2,1/4,-1,1,-9/16,-1/4,1/2,4,2-(Sqrt[3]),2+(Sqrt[3]),-4,-16/9,2}],Greater]

%//Length

It is important to remember that the parameter cc cannot be zero.