Maple Text Commands for Chapter 9 -
LOCAL AND GLOBAL BIFURCATIONS
Section 9.1: Computation of focal values for a Lienard system.
| > | restart:kstart:=2:kend:=11:
pp:=array(1..20):qq:=array(1..20):vv:=array(1..20):vx:=array(0..20): vy:=array(0..20):xx:=array(0..20,0..20): yy:=array(0..20,0..20):uu:=array(0..20,0..20): z:=array(0..20):ETA:=array(1..20): pp[1]:=y:qq[1]:=-x:vv[2]:=(x^2+y^2)/2:vx[2]:=x:vy[2]:=y: for j1 from 0 to 20 do for j2 from 0 to 20 do xx[j1,j2]:=0:yy[j1,j2]:=0:end do:end do: # Insert the coefficients for a specific Lienard system # xx[0,1]:=1:xx[2,0]:=-a2:xx[3,0]:=-a3:xx[4,0]:=-a4: yy[1,0]:=-1:yy[2,0]:=-b2:yy[3,0]:=-b3: for kloop from kstart to kend do kk:=kloop: dd1:=sum(pp[i]*vx[kk+2-i]+qq[i]*vy[kk+2-i],i=2..kk-1): pp[kk]:=sum(xx[kk-i,i]*x^(kk-i)*y^i,i=0..kk): qq[kk]:=sum(yy[kk-i,i]*x^(kk-i)*y^i,i=0..kk): vv[kk+1]:=sum(uu[kk+1,i]*x^(kk+1-i)*y^i,i=0..kk+1): d1:=y*diff(vv[kk+1],x)-x*diff(vv[kk+1],y)+pp[kk]*vx[2]+qq[kk]*vy[2]+dd1: dd:=expand(d1): if irem(kk,2)=1 then dd:=dd-ETA[kk+1]*(x^2+y^2)^((kk+1)/2): fi: dd:=numer(dd):x:=1: for i from 0 to kk+1 do z[i]:=coeff(dd,y,i); od: if kk=2 then seqn:=solve({z[0],z[1],z[2],z[3]},{uu[3,0],uu[3,1],uu[3,2],uu[3,3]}): elif kk=3 then seqn:=solve({z[0],z[1],z[2],uu[4,2],z[3],z[4]},{uu[4,0],uu[4,1],uu[4,2 ],uu[4,3],uu[4,4],ETA[4]}): elif kk=4 then seqn:=solve({z[0],z[1],z[2],z[3],z[4],z[5]},{uu[5,0],uu[5,1],uu[5,2],uu[5,3],uu[5,4],uu[5,5]}): elif kk=5 then seqn:=solve({z[0],z[1],z[2],uu[6,2]+uu[6,4],z[3],z[4],z[5],z[6]},{uu[6 ,0],uu[6,1],uu[6,2],uu[6,3],uu[6,4],uu[6,5],uu[6,6],ETA[6]}): elif kk=6 then seqn:=solve({z[0],z[1],z[2],z[3],z[4],z[5],z[6],z[7]},{uu[7,0],uu[7,1],uu[7,2],uu[7,3],uu[7,4],uu[7,5],uu[7,6],uu[7,7]}): elif kk=7 then seqn:=solve({z[0],z[1],z[2],z[3],z[4],uu[8,4],z[5],z[6],z[7],z[8]},{uu[8 ,0],uu[8,1],uu[8,2],uu[8,3],uu[8,4],uu[8,5],uu[8,6],uu[8,7],uu[8,8],ETA[8]}): elif kk=8 then seqn:=solve({z[0],z[1],z[2],z[3],z[4],z[5],z[6],z[7],z[8],z[9]},{uu[9,0],uu[9,1],uu[9,2],uu[9,3],uu[9,4],uu[9,5],uu[9,6],uu[9,7],uu[9,8],uu[9,9]}): elif kk=9 then seqn:=solve({z[0],z[1],z[2],z[3],z[4],uu[10,4]+uu[10,6],z[5],z[6],z[7],z[8],z[9],z[10]},{uu[10,0],uu[10,1],uu[10,2],uu[10,3],uu[10,4],uu[10,5],uu[10,6],uu[10,7],uu[10,8],uu[10,9],uu[10,10],ETA[10]}): elif kk=10 then seqn:=solve({z[0],z[1],z[2],z[3],z[4],z[5],z[6],z[7],z[8],z[9],z[10],z[11]},{uu[11,0],uu[11,1],uu[11,2],uu[11,3],uu[11,4],uu[11,5],uu[11,6],uu[11,7],uu[11,8],uu[11,9],uu[11,10],uu[11,11]}): elif kk=11 then seqn:=solve({z[0],z[1],z[2],z[3],z[4],z[5],z[6],uu[12,6],z[7],z[8],z[9],z[10],z[11],z[12]},{uu[12,0],uu[12,1],uu[12,2],uu[12,3],uu[12,4],uu[12,5],uu[12,6],uu[12,7],uu[12,8],uu[12,9],uu[12,10],uu[12,11],uu[12,12],ETA[12]}):fi: |
| > | assign(seqn):x:='x':i:='i':
vv[kk+1]:=sum(uu[kk+1,i]*x^(kk+1-i)*y^i,i=0..kk+1): vx[kk+1]:=diff(vv[kk+1],x):vy[kk+1]:=diff(vv[kk+1],y): ETA[kk+1]:=ETA[kk+1]:od: print(L1=numer(ETA[4])):a3:=2*a2*b2/3:print(L2=numer(ETA[6])); print(L3=numer(ETA[8]));print(L4=numer(ETA[10]));print(L5=numer(ETA[12])); |
|
|
|
|
|
|
|
|
|
(1) |
| > | with(Groebner): |
Example 4: Division algorithm for multivariate polynomials. Fix a pure lexicographical order 
, divide the polynomial 
by the ordered list of polynomials
and
| > | Reduce(x^4+y^4+z^4,[x^2+y,z^2*y-1,y-z^2],plex(x,y,z)); |
 |
(2) |
| > | Reduce(x^4+y^4+z^4,[y-z^2,z^2*y-1,x^2+y],plex(x,y,z)); |
 |
(3) |
Example 5: S-polynomials. Fix a pure lexicographical order 
, suppose that 
determine 
| > | SPolynomial(x-13*y^2-12*z^3,x^2-x*y+92*z,plex(x,y,z)); |
 |
(4) |
Example 6: Compute the reduced Groebner basis of 
with respect to
the pure lexicographical order 
.
| > | GB=Basis([x+y^2-x^3,4*x^3-12*x*y^2+x^4+2*x^2*y^2+y^4],plex(x,y,z)); |
  |
(5) |
Example 7: Compute the Groebner basis of the first five Lyapunov quantities for a Lienard system.
| > | restart:with(Groebner): |
| > | GB:=Basis([-a1,2*a2*b2-3*a3,5*b2*(2*a4-b3*a2),-5*b2*(92*b2^2*a4-99*b3^2*a2+1520*a2^2*a4-760*a2^3*b3-46*b2^2*b3*a2+198*b3*a4),-b2*(14546*b2^4*a4+105639*a2^3*b3^2+96664*a2^3*b2^2*b3-193328*a2^2*b2^2*a4-891034*a2^4*a4+445517*a2^5*b3+211632*a2*a4^2-317094*a2^2*b3*a4-44190*b2^2*b3*a4+22095*b2^2*b3^2*a2-7273*b2^4*b3*a2+5319*b3^3*a2-10638*b3^2*a4)],plex(a1,a2,a3,a4,b2,b3)); |
![[`+`(`*`(3, `*`(a3, `*`(b3))), `-`(`*`(4, `*`(b2, `*`(a4))))), `+`(`*`(2, `*`(a2, `*`(b2))), `-`(`*`(3, `*`(a3)))), a1]](images/Chap9_27.gif){kind=small} |
(6) |
Figure 9.1: Global bifurcation of a limit cycle for 
| > | with(plots):with(DEtools):
deq1:=diff(x(t),t)=y(t)+10*x(t)*(0.1-(y(t))^2): deq2:=diff(y(t),t)=-x(t)+C: bifdeq1:=(parameter)->subs(C=parameter,deq1): bifdeq2:=(parameter)->subs(C=parameter,deq2): Homoclinic:=seq(DEplot({bifdeq1('i/100-0.3'),bifdeq2('i/100-0.3')},[x(t),y(t)],0..100,[[x(0)=0.1,y(0)=0]],y=-1.5..1.5,x=-1.5..1.5,arrows=NONE,stepsize=0.1,linecolour=blue),i=0..100): Homoclinic:=subs(THICKNESS(3)=THICKNESS(0),[Homoclinic]): display(Homoclinic,insequence=true); |
 |
Example 10: A homoclinic bifurcation for 
| > | restart:with(plots):with(DEtools):
deq1:=diff(x(t),t)=y(t): deq2:=diff(y(t),t)=x(t)+(x(t))^2-x(t)*y(t)+lambda*y(t): bifdeq1:=(parameter)->subs(lambda=parameter,deq1): bifdeq2:=(parameter)->subs(lambda=parameter,deq2): Homoclinic:=seq(DEplot({bifdeq1('i/80-1.5'),bifdeq2('i/80-1.5')},[x(t),y(t)],0..80,[[x(0)=-0.4,y(0)=0]],y=-1..1,x=-2..0,arrows=NONE,stepsize=0.1,linecolour=blue),i=0..80): Homoclinic:=subs(THICKNESS(3)=THICKNESS(0),[Homoclinic]): display(Homoclinic,insequence=true); |
 |
End of Chapter 9 Commands