Calculations for Limit Cycles in a Two-Species Reaction with Wolfram Mathematica

Created by Brigita Ferčec, Ilona Nagy, Valery Romanovski, Gábor Szederkényi and János Tóth.
Last modified: 27 June, 2018.

1. Preparation for the analysis of the model

Decreasing the number of parameters

Equations (15) and (16)

Quit

du=-K2 u-2K4 (u^2)+K3 v+K5 (u^2)v;
dv=K1+K4 (u^2)-K3 v-K5 (u^2)v;
dx1=(du/(a τ))/.{u->a x, v->b y}//Expand
dy1=(dv/(b τ))/.{u->a x, v->b y}//Expand

Equation (17)

Solve[{((a b K5)/τ)==1,(((a^2) K5)/τ)==1,((a K4)/τ)==1,(((a^2)K4)/(b τ))==1},{τ,a,b}]

dx2=dx1//.{τ->(a^2)K5, b->a,a->(K4/K5)}//Expand
dy2=dy1//.{τ->(a^2)K5, b->a,a->(K4/K5)}//Expand

Equation (18)

dx3=dx2/.{((K1 (K5^2))/(K4^3))->k1,((K2 K5 )/(K4^2))->k2,((K3 K5 )/(K4^2))->k3}//Expand
dy3=dy2/.{((K1 (K5^2))/(K4^3))->k1,((K2 K5 )/(K4^2))->k2,((K3 K5 )/(K4^2))->k3}//Expand

2. Limit cycles of system (18)

Transformations

We investigate the behaviour of the trajectories at the end of the O x axis In sections 2.a) - 2.e) and at the end of the O y axis in sections 3.a) - 3.b).

1. Singular points of (18) are in the first and second quadrant

Quit

xd=-k2 x-2(x^2)+k3 y+(x^2)y;
yd=k1+(x^2)-k3 y-(x^2)y;
sol=Solve[{xd==0,yd==0},{x,y}]//FullSimplify

((4 (k1^2)+k1 k2 (k2+(Sqrt[4 k1+(k2^2)]))+4 k1 k3+k2 (k2-(Sqrt[4 k1+(k2^2)])) k3)/(2 ((k1^2)+2 k1 k3+k3 ((k2^2)+k3))))==((2k1(Sqrt[(k2^2)+4k1]))/(k2(k3-k1)+(Sqrt[(k2^2)+4k1])(k3+k1)))//Simplify

2. a) Substitution u=(y/x), z=(1/x) and time rescaling dt->(1/(z^2))dt

System (19)

ClearAll[xd,yd,x,y,k1,k2,k3,u,z,uu,zz,ud,zd];
xd=-k2 x-2(x^2)+k3 y+(x^2)y;
yd=k1+(x^2)-k3 y-(x^2)y;
uu=(y/x);zz=(1/x);
ud=(z^2)(D[uu,x]xd+D[uu,y]yd//.{x->(1/z),y->x u})//Expand
zd=(z^2)(D[zz,x]xd+D[zz,y]yd//.{x->(1/z),y->x u})//Expand

Singular points of (19)

Solve[{ud==0,zd==0},{u,z}]

2. b) Blow up the singular point (0,0) in (19) by X=u, Y=(z/u)

System (21)

ClearAll[u,z,ud,zd,k1,k2,k3,X,Y,Xd,Yd,XX,YY];
ud=-u-(u^2)+z+2 u z+k2 u (z^2)-k3 u (z^2)-k3 (u^2) (z^2)+k1 (z^3);
zd=-u z+2 (z^2)+k2 (z^3)-k3 u (z^3);
XX=u;YY=(z/u);
Xd=D[XX,u]ud+D[XX,z]zd//.{u->X,z->u Y}//Factor
Yd=D[YY,u]ud+D[YY,z]zd//.{u->X,z->u Y}//Factor

Singular points of (21)

Solve[{Xd==0,Yd==0},{X,Y}]//FullSimplify

2. c) Moving the singular point (0,0) to (0,1) using the substitution w=X, v=Y-1 and time rescaling dt->-dt

System (22)

ClearAll[X,Y,Xd,Yd,k1,k2,k3,v,w,vd,wd,vv,ww,P,Q];
Xd=-X (1+X-Y-2 X Y-k2 (X^2) (Y^2)+k3 (X^2) (Y^2)+k3 (X^3) (Y^2)-k1 (X^2) (Y^3));
Yd=-Y (-1+Y-k3 (X^2) (Y^2)+k1 (X^2) (Y^3));
ww=X;vv=Y-1;
wd=-(D[ww,X]Xd+D[ww,Y]Yd//.{X->w,Y->v+1})//Expand
vd=-(D[vv,X]Xd+D[vv,Y]Yd//.{X->w,Y->v+1})//Expand

Solution of v+Q(w,v)=0

P[w_,v_]:=wd;
Q[w_,v_]:=vd;
P[w,v]/.v->0//Factor
Q[w,v]/.v->0//Factor

ClearAll[sol];
sol=Solve[v+Q[w,v]==0,v];

P[w,v/.sol[[1]]]==P[w,v/.sol[[2]]]==P[w,v/.sol[[3]]]==P[w,v/.sol[[4]]]

P[w,v/.sol[[1]]]

2. d) Figures for system (22) when k1>k3, k1=k3, k1<k3

Quit

P[w_,v_]:=-v w-(w^2)-2 v (w^2)-k1 (w^3)-k2 (w^3)+k3 (w^3)-3 k1 v (w^3)-2 k2 v (w^3)+2 k3 v (w^3)-3 k1 (v^2) (w^3)-k2 (v^2) (w^3)+k3 (v^2) (w^3)-k1 (v^3) (w^3)+k3 (w^4)+2 k3 v (w^4)+k3 (v^2) (w^4);
Q[w_,v_]:=v+(v^2)+k1 (w^2)-k3 (w^2)+4 k1 v (w^2)-3 k3 v (w^2)+6 k1 (v^2) (w^2)-3 k3 (v^2) (w^2)+4 k1 (v^3) (w^2)-k3 (v^3) (w^2)+k1 (v^4) (w^2);

k1>k3

StreamPlot[{P[w,v],Q[w,v]}/.{k1->4,k2->1,k3->1},{w,-2,2},{v,-2,2},StreamScale->0.1,Axes->True]

k1=k3

StreamPlot[{P[w,v],Q[w,v]}/.{k1->1,k2->1,k3->1},{w,-2,2},{v,-2,2},StreamScale->0.1,Axes->True]

k1<k3

StreamPlot[{P[w,v],Q[w,v]}/.{k1->1,k2->1,k3->2},{w,-2,2},{v,-2,2},StreamScale->0.1,Axes->True]

2. e) Figure for system (19) after the blow-down

StreamPlot[{-u-(u^2)+z+2 u z+k2 u (z^2)-k3 u (z^2)-k3 (u^2) (z^2)+k1 (z^3),-u z+2 (z^2)+k2 (z^3)-k3 u (z^3)}/.{k1->1,k2->1,k3->1},{u,-1,3},{z,-1,3},StreamScale->0.1,Axes->True]

3. a) Substitution u=(x/y), z=(1/y) and time rescaling dt->(1/(z^2))dt

System (23)

ClearAll[xd,yd,x,y,k1,k2,k3,u,z,uu,zz,ud,zd];
xd=-k2 x-2(x^2)+k3 y+(x^2)y;
yd=k1+(x^2)-k3 y-(x^2)y;
uu=(x/y);zz=(1/y);
ud=(z^2)(D[uu,x]xd+D[uu,y]yd//.{y->(1/z),x->y u})//Expand
zd=(z^2)(D[zz,x]xd+D[zz,y]yd//.{y->(1/z),x->y u})//Expand

3. b) Blow-up in the u direction. Substitution u=X, z=X Y and time rescaling dt->X dt

System (25)

ClearAll[u,z,ud,zd,k1,k2,k3,X,Y,Xd,Yd,XX,YY];
ud=(u^2)+(u^3)-2 (u^2) z-(u^3) z+k3 (z^2)-k2 u (z^2)+k3 u (z^2)-k1 u (z^3);
zd=(u^2) z-(u^2) (z^2)+k3 (z^3)-k1 (z^4);
XX=u;YY=(z/u);
Xd=(1/X)(D[XX,u]ud+D[XX,z]zd//.{u->X,z->u Y})//Factor
Yd=(1/X)(D[YY,u]ud+D[YY,z]zd//.{u->X,z->u Y})//Factor

Theorem 3

4. a) Singular points of (18) if x0=1.

Quit

ClearAll[xd,yd,x,y,k1,k2,k3,sol,x1,y1,x1d,y1d];
xd=-k2 x-2(x^2)+k3 y+(x^2)y;
yd=k1+(x^2)-k3 y-(x^2)y;
Solve[{xd==0,yd==0}/.x->1,{k1,y}]//FullSimplify

4. b) System (26): the singular point (if x0=1) is shifted into (0,0)

ClearAll[x0,y0,x,y,xd,yd,x1,y1,x1d,y1d,k1,k2,k3,k4,k5];
xd=-k2 x-2(x^2)+k3 y+(x^2)y;
yd=k1+(x^2)-k3 y-(x^2)y;
x0=1;y0=((2+k2)/(1+k3));
k1=k2+1;
xx1=x-x0;yy1=y-y0;
x1d=D[xx1,x]xd+D[xx1,y]yd/.{x->x1+x0,y->y1+y0}//Factor
y1d=D[yy1,x]xd+D[yy1,y]yd/.{x->x1+x0,y->y1+y0}//Factor

Solve[{x1d==0,y1d==0},{x1,y1}]

4. c) The Jacobian at the origin

Jac=D[{x1d,y1d},{{x1,y1}}];
JacOrigin=Jac/.{x1->0,y1->0}//FullSimplify

trace=Tr[JacOrigin]//Factor
Solve[trace==0,k2]//Factor

k2=-((1+6 k3+(k3^2))/(-1+k3));
Eigenvalues[JacOrigin]//FullSimplify//PowerExpand(*if k3<1 then these are pure imaginary*)

4. d) System (29)

The matrix S that transforms JacOrigin into Jordan canonical form (change from x1, y1 to u,v)

a11=JacOrigin[[1,1]]//Factor;
a21=JacOrigin[[2,1]]//Factor;
alfa=0;
beta=-(((1+k3) (Sqrt[3+k3]))/(Sqrt[1-k3]));(*if k3<1 then this is real*)
S={{a11-alfa,-beta},{a21,0}}

ClearAll[uu,vv,sol];
{uu,vv}=Inverse[S].{x1,y1}//FullSimplify
sol=Solve[{u,v}==Inverse[S].{x1,y1},{x1,y1}]//FullSimplify

ClearAll[ud,vd];
ud=D[uu,x1]x1d+D[uu,y1]y1d/.sol[[1]]//FullSimplify//Expand;
vd=D[vv,x1]x1d+D[vv,y1]y1d/.sol[[1]]//FullSimplify//Expand;

JacOriginNew=D[{ud,vd},{{u,v}}]/.{u->0,v->0}//FullSimplify

Eigenvalues[JacOriginNew]

The transformed system after multiplying by S

ud
vd

Time rescaling dt->-(((1+k3) (Sqrt[3+k3]))/(Sqrt[1-k3]))dt

ClearAll[BB];
BB=JacOriginNew[[2,1]]

ClearAll[Ud,Vd];
Ud=ud/BB//FullSimplify
Vd=vd/BB//FullSimplify

Coefficient[Ud,v]/.u->0//FullSimplify
Coefficient[Vd,u]/.v->0//FullSimplify

System (29)

Ud//Factor

Vd//Factor

5. a) Lyapunov’s theorem, RHS=g1((((x^2)+(y^2)))^2)+g2((((x^2)+(y^2)))^3)

System (29)

Quit

Ud=(1/(2 (Sqrt[3-2 k3-(k3^2)])))(3 (u^2)-4 k3 (u^2)+(k3^2) (u^2)+2 (u^3)-2 (k3^2) (u^3)-2 (Sqrt[3-2 k3-(k3^2)]) v+2 (Sqrt[3-2 k3-(k3^2)]) u v-2 k3 (Sqrt[3-2 k3-(k3^2)]) u v+4 (Sqrt[3-2 k3-(k3^2)]) (u^2) v+4 k3 (Sqrt[3-2 k3-(k3^2)]) (u^2) v-3 (v^2)-4 k3 (v^2)-(k3^2) (v^2)+6 u (v^2)+8 k3 u (v^2)+2 (k3^2) u (v^2));
Vd=(1/(2 (-1+k3) (3+k3)))(-6 u+4 k3 u+2 (k3^2) u-11 (u^2)+13 k3 (u^2)-(k3^2) (u^2)-(k3^3) (u^2)-6 (u^3)-2 k3 (u^3)+6 (k3^2) (u^3)+2 (k3^3) (u^3)-10 (Sqrt[3-2 k3-(k3^2)]) u v+8 k3 (Sqrt[3-2 k3-(k3^2)]) u v+2 (k3^2) (Sqrt[3-2 k3-(k3^2)]) u v-12 (Sqrt[3-2 k3-(k3^2)]) (u^2) v-16 k3 (Sqrt[3-2 k3-(k3^2)]) (u^2) v-4 (k3^2) (Sqrt[3-2 k3-(k3^2)]) (u^2) v+3 (v^2)+19 k3 (v^2)+9 (k3^2) (v^2)+(k3^3) (v^2)-18 u (v^2)-30 k3 u (v^2)-14 (k3^2) u (v^2)-2 (k3^3) u (v^2));

pp=Ud/.{u->x,v->y};
qq=Vd/.{u->x,v->y};
Ser[s_]:=Plus@@Table[(x^i)(y^(s-i))p[i,s-i],{i,0,s}];
Hom[s_]:=Table[p[s-i,i],{i,0,s}];
hh=(x^2)+  (y^2)+Sum[Ser[i],{i,3,10}];
Lie=D[hh,x] pp+D[hh,y] qq//Expand;
RHS=g1((((x^2)+(y^2)))^2)+g2((((x^2)+(y^2)))^3)//Expand;
coefflist[m_,polynomial_]:=Module[{poly=polynomial,list,d=m,listfirst,full,complistfirst,complist},
list=Select[CoefficientRules[poly,{x,y}],MemberQ[Table[{d-i,i},{i,0,d}],#[[1]]]&];
listfirst=Table[list[[i,1]],{i,Length[list]}];
full=Table[{d-i,i},{i,0,d}];
complistfirst=DeleteCases[full,x_/;MemberQ[listfirst,x]];
complist=Table[complistfirst[[i]]->0,{i,Length[complistfirst]}];
Sort[Union[list, complist],#1[[1,1]]>#2[[1,1]]&]
]

degree 1

ClearAll[m,LS1,aa,bb];m=1;
aa=coefflist[m,Lie]
bb=coefflist[m,RHS]
LS1=Table[aa[[i,2]]-bb[[i,2]],{i,Length[aa]}]

degree 2

ClearAll[m,LS2,aa,bb];m=2;
aa=coefflist[m,Lie]//FullSimplify
bb=coefflist[m,RHS]//FullSimplify
LS2=Table[aa[[i,2]]-bb[[i,2]],{i,Length[aa]}]

degree 3

ClearAll[m,LS3,aa,bb];m=3;
aa=coefflist[m,Lie]//Simplify;
bb=coefflist[m,RHS]//Simplify;
LS3=Table[aa[[i,2]]-bb[[i,2]],{i,Length[aa]}]

sol3=Solve[LS3==0,Hom[3]]//FullSimplify

{p[3,0],p[2,1],p[1,2],p[0,3]}={p[3,0],p[2,1],p[1,2],p[0,3]}/.sol3[[1]]//Factor

LS3//FullSimplify

degree 4: p[4,0] is arbitrary

ClearAll[m,LS4,aa,bb];m=4;
aa=coefflist[m,Lie]//FullSimplify;
bb=coefflist[m,RHS]//FullSimplify;
LS4=Table[aa[[i,2]]-bb[[i,2]],{i,Length[aa]}]//FullSimplify//Factor

sol4=Solve[LS4==0,{g1,p[3,1],p[2,2],p[1,3],p[0,4]}]//Factor

{g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}={g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}/.sol4[[1]];

LS4//Factor

5. b) Lyapunov’s theorem, RHS=g1 (x^4)+g2 (x^6)+...

Here the result for (Subscript[g, 1]) is a constant multiple of the previous result.

System (29)

Quit

Ud=(1/(2 (Sqrt[3-2 k3-(k3^2)])))(3 (u^2)-4 k3 (u^2)+(k3^2) (u^2)+2 (u^3)-2 (k3^2) (u^3)-2 (Sqrt[3-2 k3-(k3^2)]) v+2 (Sqrt[3-2 k3-(k3^2)]) u v-2 k3 (Sqrt[3-2 k3-(k3^2)]) u v+4 (Sqrt[3-2 k3-(k3^2)]) (u^2) v+4 k3 (Sqrt[3-2 k3-(k3^2)]) (u^2) v-3 (v^2)-4 k3 (v^2)-(k3^2) (v^2)+6 u (v^2)+8 k3 u (v^2)+2 (k3^2) u (v^2));
Vd=(1/(2 (-1+k3) (3+k3)))(-6 u+4 k3 u+2 (k3^2) u-11 (u^2)+13 k3 (u^2)-(k3^2) (u^2)-(k3^3) (u^2)-6 (u^3)-2 k3 (u^3)+6 (k3^2) (u^3)+2 (k3^3) (u^3)-10 (Sqrt[3-2 k3-(k3^2)]) u v+8 k3 (Sqrt[3-2 k3-(k3^2)]) u v+2 (k3^2) (Sqrt[3-2 k3-(k3^2)]) u v-12 (Sqrt[3-2 k3-(k3^2)]) (u^2) v-16 k3 (Sqrt[3-2 k3-(k3^2)]) (u^2) v-4 (k3^2) (Sqrt[3-2 k3-(k3^2)]) (u^2) v+3 (v^2)+19 k3 (v^2)+9 (k3^2) (v^2)+(k3^3) (v^2)-18 u (v^2)-30 k3 u (v^2)-14 (k3^2) u (v^2)-2 (k3^3) u (v^2));

pp=Ud/.{u->x,v->y};
qq=Vd/.{u->x,v->y};
Ser[s_]:=Plus@@Table[(x^i)(y^(s-i))p[i,s-i],{i,0,s}];
Hom[s_]:=Table[p[s-i,i],{i,0,s}];
hh=(x^2)+  (y^2)+Sum[Ser[i],{i,3,10}];
Lie=D[hh,x] pp+D[hh,y] qq//Expand;
RHS=g1 (x^4)+g2 (x^6)//Expand;
coefflist[m_,polynomial_]:=Module[{poly=polynomial,list,d=m,listfirst,full,complistfirst,complist},
list=Select[CoefficientRules[poly,{x,y}],MemberQ[Table[{d-i,i},{i,0,d}],#[[1]]]&];
listfirst=Table[list[[i,1]],{i,Length[list]}];
full=Table[{d-i,i},{i,0,d}];
complistfirst=DeleteCases[full,x_/;MemberQ[listfirst,x]];
complist=Table[complistfirst[[i]]->0,{i,Length[complistfirst]}];
Sort[Union[list, complist],#1[[1,1]]>#2[[1,1]]&]
]

degree 1

ClearAll[m,LS1,aa,bb];m=1;
aa=coefflist[m,Lie]
bb=coefflist[m,RHS]
LS1=Table[aa[[i,2]]-bb[[i,2]],{i,Length[aa]}]

degree 2

ClearAll[m,LS2,aa,bb];m=2;
aa=coefflist[m,Lie]//FullSimplify
bb=coefflist[m,RHS]//FullSimplify
LS2=Table[aa[[i,2]]-bb[[i,2]],{i,Length[aa]}]

degree 3

ClearAll[m,LS3,aa,bb];m=3;
aa=coefflist[m,Lie]//Simplify;
bb=coefflist[m,RHS]//Simplify;
LS3=Table[aa[[i,2]]-bb[[i,2]],{i,Length[aa]}]

sol3=Solve[LS3==0,Hom[3]]//FullSimplify

{p[3,0],p[2,1],p[1,2],p[0,3]}={p[3,0],p[2,1],p[1,2],p[0,3]}/.sol3[[1]]//Factor

LS3//FullSimplify

degree 4: p[4,0] is arbitrary

ClearAll[m,LS4,aa,bb];m=4;
aa=coefflist[m,Lie]//FullSimplify;
bb=coefflist[m,RHS]//FullSimplify;
LS4=Table[aa[[i,2]]-bb[[i,2]],{i,Length[aa]}]//FullSimplify//Factor

sol4=Solve[LS4==0,{g1,p[3,1],p[2,2],p[1,3],p[0,4]}]//Factor

{g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}={g1,p[4,0],p[3,1],p[2,2],p[1,3],p[0,4]}/.sol4[[1]];

LS4//Factor

Plotting the limit cycles

Preparations

Quit

SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ParametricPlot,ListPlot,ListLinePlot};

SetDirectory[NotebookDirectory[]];
SetOptions[#,AxesStyle->Arrowheads[Automatic]]&/@{Plot,ListPlot,ListLinePlot,ListLogLogPlot,ParametricPlot,DateListPlot,DiscretePlot};

LaunchKernels[];

The function creating the plots

ClearAll[k,p,q,x,y,g];
(Subscript[k, 3])=(3/10);(* 0 < k3 < 1 *)
(Subscript[k, 4])=1;
(Subscript[k, 5])=1;
(Subscript[k, 2])=-((((Subscript[k, 3])^2)+4 (Subscript[k, 3]) (Subscript[k, 4])+2 (Subscript[k, 3]) (Subscript[k, 5])+((Subscript[k, 5])^2))/((Subscript[k, 3])-(Subscript[k, 5])))-(1/10000);
(Subscript[k, 1])=(Subscript[k, 2])+(Subscript[k, 4]);
p[x_,y_]:=-(Subscript[k, 2]) x-2(Subscript[k, 4]) (x^2)+(Subscript[k, 3]) y+(Subscript[k, 5]) (x^2)y;
q[x_,y_]:=(Subscript[k, 1])+(Subscript[k, 4]) (x^2)-(Subscript[k, 3]) y-(Subscript[k, 5]) (x^2)y;

ClearAll[nsol,ev,plotter];nsol=First@NSolve[Join@@Thread/@{{p[x,y],q[x,y]}==0,{x,y}>0},{x,y},20];
ev=Eigenvalues[D[{p[x,y],q[x,y]},{{x,y}}]/.nsol];
plotter[τ_,shift_,ag_:Automatic,pg_:Automatic,pp_:1000,ar_:Automatic,opts___]:=Module[{startingpoint,sys,solution},
startingpoint=({x,y}/.nsol)+shift;
sys:=NDSolveValue[Join[{u'[t]==p[u[t],v[t]],v'[t]==q[u[t],v[t]]},Thread[{u[0],v[0]}==startingpoint]],
{u,v},{t,τ},AccuracyGoal->ag,PrecisionGoal->pg,opts];
trafo[point_]:=(1/(shift[[2]]))(point-startingpoint);
solution[t_]:=trafo[Through[sys[t]]];
 {ParametricPlot[Evaluate[solution[t]],{t,0,τ},Epilog->{Red,PointSize[0.05],Point[{0,0}],Point[trafo[{x,y}/.nsol]]},PlotRange->All,PlotPoints->pp,AspectRatio->ar,AxesLabel->{x,y},LabelStyle->Directive[14],ImageSize->250],
Plot[Evaluate[solution[t][[1]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,x},LabelStyle->Directive[12],ImageSize->250],
Plot[Evaluate[solution[t][[2]]],{t,0,τ},PlotRange->All,PlotPoints->pp,AxesLabel->{t,y},LabelStyle->Directive[12],ImageSize->250]}]

Figures with components of the solution

Limit cycle, trajectory going inward

Distance from the singular point: 1

Figure1=plotter[200,{0,1}]

Limit cycle, trajectory going outward

Distance from the singular point: (10^(-5))

Figure2=plotter[30,{0,0.00001}]