The van der Pol oscillator An example for a nonlinear oscillator possessing a unique limit cycle solution 1) The tangent field plot restart: with(plots): with(DEtools): alpha:=2: ps1:=diff(x(t),t)=y(t); ps2:=diff(y(t),t)=-alpha*(x(t)*x(t)-1)*y(t)-x(t); dfieldplot([ps1,ps2],[x,y],t=0..1,x=-4..+4,y=-5..+5,dirgrid=[20,20],color=blue,axes=boxed,title=`van der Pol oscillator: Tangent field plot`); 2) Integration of the differential equation The van der Pol equation is an example of a stiff system of differential equations (roughly: the solutions show several widely differing time scales) which makes numerical solution more difficult. The option 'stiff=true' tells MAPLE to use an integration method suitable for such problems (the Rosenbrock method). restart: with(plots): ps1:=diff(x(t),t)=y(t): ps2:=diff(y(t),t)=-alpha*(x(t)*x(t)-1)*y(t)-x(t): alpha:=20.: # system parameter x0:=0.5: y0:=0.: # initial conditions sol:=dsolve({ps1,ps2,x(0)=x0,y(0)=y0},{x(t),y(t)},type=numeric,stiff=true,output=listprocedure): plot1:=odeplot(sol,[t,x(t)],0..300,numpoints=5000,color=blue): display(plot1,view=[0..200,-5..+5],axes=boxed,labels=[t,x], title= `van der Pol oscillator`); 3) Phase space plot plot2:=odeplot(sol,[x(t),y(t)],0..100,numpoints=30000,color=blue): display(plot2,view=[-4..+4,-30..+30],axes=boxed,labels=[x,y], title= `van der Pol oscillator`); 4) Modified formulation of van der Pol's equation To analyze the limit \316\261>>1 an alternative formulation of the differential equations is used. In this formulation the limit cycle consists of four segments: two horizontal straight lines and two segments of the cubic function y=g(x). restart: with(plots): g:= x-> -x+x^3/3: pmod1:=diff(x(t),t)=alpha*(z(t)-g(x(t))): pmod2:=diff(z(t),t)=-x(t)/alpha: alpha:=20: # system parameter x0:=-2.: z0:=2.: # initial conditions sol:=dsolve({pmod1,pmod2,x(0)=x0,z(0)=z0},{x(t),z(t)},type=numeric,stiff=true,output=listprocedure): plot2:=odeplot(sol,[x(t),z(t)],0..200,numpoints=1000,color=blue,thickness=2): gfunct:=plot(g(x),x=-3..+3,color=red,thickness=3): display([plot2,gfunct],view=[-3..+3,-3..+3],axes=boxed,labels=[x,z], title= `Modified phase plane`); Floquet stability analysis The stability of periodic trajectories according to Floquet can be analyzed as follows: a) Calculate the reference trajectory LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYlLUYvNiVRInJGJ0YyRjUvJSxwbGFjZWhvbGRlckdGNC9GNlEnbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9JK21zZW1hbnRpY3NHRiRRJ2F0b21pY0YnLUkobWZlbmNlZEdGJDYkLUYjNiQtRi82JVEidEYnRjJGNUY/Rj8tSSNtb0dGJDYtUSIuRidGPy8lJmZlbmNlR1EmZmFsc2VGJy8lKnNlcGFyYXRvckdGVS8lKXN0cmV0Y2h5R0ZVLyUqc3ltbWV0cmljR0ZVLyUobGFyZ2VvcEdGVS8lLm1vdmFibGVsaW1pdHNHRlUvJSdhY2NlbnRHRlUvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0Zeb0Y/ b) Set up the linearized equation of mation LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYuLUkmbWZyYWNHRiQ2KC1JI21vR0YkNi1RMCZEaWZmZXJlbnRpYWxEO0YnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy8lJmZlbmNlR1EmdW5zZXRGJy8lKnNlcGFyYXRvckdGNy8lKXN0cmV0Y2h5R0Y3LyUqc3ltbWV0cmljR0Y3LyUobGFyZ2VvcEdGNy8lLm1vdmFibGVsaW1pdHNHRjcvJSdhY2NlbnRHRjcvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZGLUYjNiZGLi1GLzYtUSJ+RidGMi9GNlEmZmFsc2VGJy9GOUZPL0Y7Rk8vRj1GTy9GP0ZPL0ZBRk8vRkNGT0ZERkctSSNtaUdGJDYlUSJ0RicvJSdpdGFsaWNHUSV0cnVlRicvRjNRJ2l0YWxpY0YnRjIvJS5saW5ldGhpY2tuZXNzR1EiMUYnLyUrZGVub21hbGlnbkdRJ2NlbnRlckYnLyUpbnVtYWxpZ25HRl5vLyUpYmV2ZWxsZWRHRk9GSy1GVzYlUScmIzk4MTtGJy9GZW5GT0YyRkstRi82LVEiPUYnRjJGTkZQRlFGUkZTRlRGVS9GRVEsMC4yNzc3Nzc4ZW1GJy9GSEZbcEZLLUZXNiVRIk1GJ0ZaRmduRktGY29GS0ZLRjI=with the Jacobian matrix M. c) Integrate this differential equation over one oscillation period T to generate a fundamental system of N solutions LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYmLUklbXN1YkdGJDYmLUkjbWlHRiQ2JlEnJiM5ODE7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYlLUYvNiVRImlGJy9GM0Y3L0Y5USdpdGFsaWNGJy8lLHBsYWNlaG9sZGVyR0Y3RjgvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9JK21zZW1hbnRpY3NHRiRRJ2F0b21pY0YnLUkobWZlbmNlZEdGJDYkLUYjNiQtRi82JVEiVEYnRkBGQUY4RjgtSSNtb0dGJDYtUSIuRidGOC8lJmZlbmNlR0Y0LyUqc2VwYXJhdG9yR0Y0LyUpc3RyZXRjaHlHRjQvJSpzeW1tZXRyaWNHRjQvJShsYXJnZW9wR0Y0LyUubW92YWJsZWxpbWl0c0dGNC8lJ2FjY2VudEdGNC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRmFvRjg= d) Evaluate the eigenvalues LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYmLUkjbWlHRiQ2JlEnJiM5NTU7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYlLUYvNiVRImlGJy9GM0Y3L0Y5USdpdGFsaWNGJy8lLHBsYWNlaG9sZGVyR0Y3RjgvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9JK21zZW1hbnRpY3NHRiRRJ2F0b21pY0YnLUkjbW9HRiQ2LVEifkYnRjgvJSZmZW5jZUdGNC8lKnNlcGFyYXRvckdGNC8lKXN0cmV0Y2h5R0Y0LyUqc3ltbWV0cmljR0Y0LyUobGFyZ2VvcEdGNC8lLm1vdmFibGVsaW1pdHNHRjQvJSdhY2NlbnRHRjQvJSdsc3BhY2VHUSYwLjBlbUYnLyUncnNwYWNlR0ZpbkY4of the matrix of fundamental solutions C = (LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYqLUklbXN1YkdGJDYmLUkjbWlHRiQ2JlEpJnZhcnBoaTtGJy8lJ2l0YWxpY0dRJmZhbHNlRicvJTZzZWxlY3Rpb24tcGxhY2Vob2xkZXJHUSV0cnVlRicvJSxtYXRodmFyaWFudEdRJ25vcm1hbEYnLUYjNiUtSSNtbkdGJDYkUSIxRidGOC8lMGZvbnRfc3R5bGVfbmFtZUdRJVRleHRGJ0Y4LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvSSttc2VtYW50aWNzR0YkUSdhdG9taWNGJy1JKG1mZW5jZWRHRiQ2JC1GIzYlLUYvNiVRIlRGJy9GM0Y3L0Y5USdpdGFsaWNGJ0ZBRjhGOC1JI21vR0YkNi1RIixGJ0Y4LyUmZmVuY2VHRjQvJSpzZXBhcmF0b3JHRjcvJSlzdHJldGNoeUdGNC8lKnN5bW1ldHJpY0dGNC8lKGxhcmdlb3BHRjQvJS5tb3ZhYmxlbGltaXRzR0Y0LyUnYWNjZW50R0Y0LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdRLDAuMzMzMzMzM2VtRictRiw2JkYuLUYjNiUtRj42JFEiMkYnRjhGQUY4RkRGR0ZKLUYvNiNRIUYnRkFGOA==). If all eigenvalue satisfies |LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYoLUklbXN1YkdGJDYmLUkjbWlHRiQ2KFEnJiM5NTU7RicvJSVib2xkR1EldHJ1ZUYnLyUnaXRhbGljR0Y0LyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR0Y0LyUsbWF0aHZhcmlhbnRHUSxib2xkLWl0YWxpY0YnLyUrZm9udHdlaWdodEdRJWJvbGRGJy1GIzYoLUYvNidRImlGJ0YyRjVGOUY8RjJGNS8lLHBsYWNlaG9sZGVyR0Y0RjlGPC8lL3N1YnNjcmlwdHNoaWZ0R1EiMEYnL0krbXNlbWFudGljc0dGJFEnYXRvbWljRictRi82I1EhRidGMkY1RjlGPA==| \342\211\244 1 the reference trajectory is stable. The nonlinear differential equations restart: with(plots): with(DEtools): with(LinearAlgebra): ps1:=diff(xr(t),t)=yr(t); ps2:=diff(yr(t),t)=-alpha*(xr(t)*xr(t)-1)*yr(t)-xr(t); alpha := 2.; To make sure that we are on the limit cycle we start with arbitrary initial values and integrate over some time interval. xr0:=1.: yr0:=1.: sol:=dsolve({ps1,ps2,xr(0)=xr0,yr(0)=yr0},{xr(t),yr(t)},type=numeric,maxfun=200000,stiff=true,relerr=1e-11,output=listprocedure): start := sol(50*alpha); Now we start with these initial conditions xr0:=rhs(start[2]): yr0:=rhs(start[3]): sol:=dsolve({ps1,ps2,xr(0)=xr0,yr(0)=yr0},{xr(t),yr(t)},type=numeric,stiff=true,relerr=1e-11,output=listprocedure): x1:=subs(sol,xr(t)): y1:=subs(sol,yr(t)): plot1:=odeplot(sol,[t,xr(t)],0..4*alpha,numpoints=1000,color=blue): display(plot1,view=[0..4*alpha,-5..+5],axes=boxed,labels=[t,x], title= `van der Pol oscillator`); Since we need the oscillation period we look for the solution of the equation LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYyLUklbXN1YkdGJDYmLUkjbWlHRiQ2JVEieEYnLyUnaXRhbGljR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdpdGFsaWNGJy1GIzYmLUYvNiVRInJGJ0YyRjUvJSxwbGFjZWhvbGRlckdGNC8lMGZvbnRfc3R5bGVfbmFtZUdRJVRleHRGJy9GNlEnbm9ybWFsRicvJS9zdWJzY3JpcHRzaGlmdEdRIjBGJy9JK21zZW1hbnRpY3NHRiRRJ2F0b21pY0YnLUkobWZlbmNlZEdGJDYkLUYjNiUtRi82JVEiVEYnRjJGNUY/RkJGQi1JI21vR0YkNi1RIn5GJ0ZCLyUmZmVuY2VHUSZmYWxzZUYnLyUqc2VwYXJhdG9yR0ZYLyUpc3RyZXRjaHlHRlgvJSpzeW1tZXRyaWNHRlgvJShsYXJnZW9wR0ZYLyUubW92YWJsZWxpbWl0c0dGWC8lJ2FjY2VudEdGWC8lJ2xzcGFjZUdRJjAuMGVtRicvJSdyc3BhY2VHRmFvLUZTNi1RIj1GJ0ZCRlZGWUZlbkZnbkZpbkZbb0Zdby9GYG9RLDAuMjc3Nzc3OGVtRicvRmNvRmhvRlJGUi1GLDYmRi4tRiM2JUY6Rj9GQkZERkctRks2JC1GIzYlLUkjbW5HRiQ2JFEiMEYnRkJGP0ZCRkJGUi1GUzYtUSIuRidGQkZWRllGZW5GZ25GaW5GW29GXW9GX29GYm9GUkZSLUknbXNwYWNlR0YkNiYvJSdoZWlnaHRHUSYwLjBleEYnLyUmd2lkdGhHUSYwLjBlbUYnLyUmZGVwdGhHRl5xLyUqbGluZWJyZWFrR1EobmV3bGluZUYnLUYvNiNRIUYnRj9GQg==The fsolve command of MAPLE is not too clever, so we have to specifiy a good search interval. T:=fsolve(x1(t)=xr0,t=7..8); T/alpha; The solution really is periodic: x1(T);y1(T); The linearized differential equations. Note that the reference trajectory (x1(t),y1(t)) enters through the time-dependent coefficients. ld1:=diff(phi1(t),t)=phi2(t); ld2:=diff(phi2(t),t)=(-1-2*alpha*x1(t)*y1(t))*phi1(t) - alpha*(x1(t)*x1(t)-1)*phi2(t); The first fundamental solution of the linear system: sollin:=dsolve({ld1,ld2,phi1(0)=1.,phi2(0)=0.},{phi1(t),phi2(t)},type=numeric,relerr=1e-11,output=listprocedure): phi1T := sollin(T); Let us plot this solution. plotlin1:=odeplot(sollin,[t,phi1(t)],0..T,numpoints=5000,color=blue): plotlin2:=odeplot(sollin,[t,phi2(t)],0..T,numpoints=5000,color=red): display(plotlin1,plotlin2,view=[0..T,-20..+20],axes=boxed,labels=[t,phi], title= `Fundamental solution 1`); The second fundamental solution: sollin:=dsolve({ld1,ld2,phi1(0)=0.,phi2(0)=1.},{phi1(t),phi2(t)},type=numeric,relerr=1e-11,output=listprocedure): phi2T := sollin(T); plotlin3:=odeplot(sollin,[t,phi1(t)],0..T,numpoints=5000,color=blue): plotlin4:=odeplot(sollin,[t,phi2(t)],0..T,numpoints=5000,color=red): display(plotlin3,plotlin4,view=[0..T,-20..+20],axes=boxed,labels=[t,phi], title= `Fundamental solution 2`); The two solution vectors are put into the 'monodromic matrix' C: C := Matrix(2,2,[[rhs(phi1T[2]),rhs(phi2T[2])],[rhs(phi1T[3]),rhs(phi2T[3])]]); Eigenvalues(C); One of the eigenvalues is LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYpLUkjbWlHRiQ2JlEoZXhhY3RseUYnLyUnaXRhbGljR1EmZmFsc2VGJy8lNnNlbGVjdGlvbi1wbGFjZWhvbGRlckdRJXRydWVGJy8lLG1hdGh2YXJpYW50R1Enbm9ybWFsRictSSNtb0dGJDYuUSJ+RidGMkY1LyUmZmVuY2VHRjEvJSpzZXBhcmF0b3JHRjEvJSlzdHJldGNoeUdGMS8lKnN5bW1ldHJpY0dGMS8lKGxhcmdlb3BHRjEvJS5tb3ZhYmxlbGltaXRzR0YxLyUnYWNjZW50R0YxLyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGTC1JJW1zdWJHRiQ2Ji1GLDYmUScmIzk1NTtGJ0YvRjJGNS1GIzYlLUkjbW5HRiQ2JFEiMUYnRjUvJSxwbGFjZWhvbGRlckdGNEY1LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvSSttc2VtYW50aWNzR0YkUSdhdG9taWNGJy1GOTYtUSI9RidGNUY8Rj5GQEZCRkRGRkZIL0ZLUSwwLjI3Nzc3NzhlbUYnL0ZORmFvRlctRjk2LUY7RjVGPEY+RkBGQkZERkZGSEZKRk1GNQ==(within numerical accuracy) as it should be for an autonomous system. The second eigenvalue LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYlLUklbXN1YkdGJDYmLUkjbWlHRiQ2JlEnJiM5NTU7RicvJSdpdGFsaWNHUSZmYWxzZUYnLyU2c2VsZWN0aW9uLXBsYWNlaG9sZGVyR1EldHJ1ZUYnLyUsbWF0aHZhcmlhbnRHUSdub3JtYWxGJy1GIzYlLUkjbW5HRiQ2JFEiMkYnRjgvJSxwbGFjZWhvbGRlckdGN0Y4LyUvc3Vic2NyaXB0c2hpZnRHUSIwRicvSSttc2VtYW50aWNzR0YkUSdhdG9taWNGJy1JI21vR0YkNi1RIn5GJ0Y4LyUmZmVuY2VHRjQvJSpzZXBhcmF0b3JHRjQvJSlzdHJldGNoeUdGNC8lKnN5bW1ldHJpY0dGNC8lKGxhcmdlb3BHRjQvJS5tb3ZhYmxlbGltaXRzR0Y0LyUnYWNjZW50R0Y0LyUnbHNwYWNlR1EmMC4wZW1GJy8lJ3JzcGFjZUdGZ25GOA==is much smaller than 1, proving that the periodic solution of the van der Pol oscillator has very high stability. LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=