ME 406Example of Stable and Unstable ManifoldssysidMathematica 6.0.3, DynPac 11.01, 1ê13ê2009‡IntroductionIn this notebook, we study an example given in Differential Equations and Dynamical Systems,Lawrence Perko, second edition, Springer-Verlag, 1996, p. 111. The equations are (1)x° = -x - y2, y° = x2 + y .This system has an equilibrium point at (0,0). We will construct both the linear and nonlinear stable and unsta-ble manifolds in the vicinity of the origin. Then we will construct the global stable and unstable manifolds. ‡System DefinitionWe define the system for DynPac.setstate[{x,y}];setparm[{}];slopevec = 9-x - y2, y + x2=;parmval = {};sysname = "nlsystem";intreset;plotreset;imsize=250;‡Equilibrium PointWe start by analyzing the equilibrium at the origin. We first use classify2D, and then we use eigers toobtain the eigenvectors and eigenvalues at the origin.classify2D[{0,0}]Abbreviations used in classify2D.L = linear, NL = nonlinear, R2 = repeated root.Z1 = one zero root, Z2 = two zero roots.This message printed once.unstable - saddleeigers = eigsys@80, 0<D88-1, 1<, 881, 0<, 80, 1<<<‡Stable and Unstable Manifolds for Nonlinear SystemWe will use blue for stable manifolds, red for unstable manifolds.The stable manifold corresponds to the negative eigenvalue. To construct it, we integrate backwards intime, starting very near the equilibrium, displaced from the equilibrium along the eigenvector associated with theeigenvalue -1. The solution going backwards in time could grow very rapidly, so we avoid the possibility ofoverflow and crash by using range checking. We limit the range to a box from -1 to 1 in both variables. We addan arrow to the mid-point of each orbit. We remove the axes (and put a frame around the picture) so that thetangencies can be better seen when we compare linear and nonlinear systems.ranger = {{-1.4,1.4},{-1.4,1.4}};rangeflag = True;t0 = 0.0;h = -0.02;nsteps = 1000;initset = {0.005*{1,0},-0.005*{1,0}};arrowflag = True;arrowvec = 81 ê 2<;asprat = 1.0; plrange = 88-1.1, 1.1<, 8-1.1, 1.1<<;axon = False; frameon = True;setcolor@8Blue<D;2 manifold.nbstableport = portrait[initset,t0,h,nsteps,1,2] -1.0-0.50.00.51.0-1.0-0.50.00.51.0xynlsystemNow we construct the unstable manifold.setcolor@8Red<D;initset = {0.005*{0,1},-0.005*{0,1}};h = 0.02;manifold.nb 3unstableport = portrait[initset,t0,h,nsteps,1,2] -1.0-0.50.00.51.0-1.0-0.50.00.51.0xynlsystemNow we combine these.nlman = show@stableport, unstableportD-1.0-0.50.00.51.0-1.0-0.50.00.51.0xynlsystem‡Stable and Unstable Manifolds for Linearized System4 manifold.nb‡Stable and Unstable Manifolds for Linearized SystemNow we construct the linearized system and then go through the same manifold constructions for it. Wefirst save the nonlinear system under the name nlsystem.savesys[nlsystem];We continue to use {x,y} for the state vector, and both parmval and parmvec remain empty lists. Thusfor the linear system, we must redefine only the slope function. It is the derivative matrix at equilibrium timesthe state vector.slopevec = Dot[dermatval[{0,0}],{x,y}]8-x, y<sysname = "linsystem";Now we construct the linear stable manifold.ranger = {{-1.4,1.4},{-1.4,1.4}};rangeflag = True;t0 = 0.0;h = -0.02;nsteps = 1000;initset = {0.005*{1,0},-0.005*{1,0}};arrowflag = True;arrowvec = 81 ê 2<;asprat = 1.0; plrange = 88-1.1, 1.1<, 8-1.1, 1.1<<;axon = False;setcolor@8Blue<D;manifold.nb 5linstableport = portrait[initset,t0,h,nsteps,1,2] -1.0-0.50.00.51.0-1.0-0.50.00.51.0xylinsystemNow we construct the unstable manifold.setcolor@8Red<D;initset = {0.005*{0,1},-0.005*{0,1}};h = 0.02;6 manifold.nblinunstableport = portrait[initset,t0,h,nsteps,1,2] -1.0-0.50.00.51.0-1.0-0.50.00.51.0xylinsystemWe combine the last two pictures.lman = show@linstableport, linunstableportD-1.0-0.50.00.51.0-1.0-0.50.00.51.0xylinsystemmanifold.nb 7Now for the finale, showing both the linear and nonlinear manifolds and illustrating the tangencies at theequilibrium point.labon = "linsystem & nlsystem";show@lman, nlmanD-1.0-0.50.00.51.0-1.0-0.50.00.51.0xylinsystem & nlsystem‡Global Manifolds for Nonlinear SystemNow we restore the nonlinear system, and we try to find the global stable and unstable manifolds. Wewill get some surprises in that effort.restoresys@nlsystemD;We enlarge our plotting window, and we begin with the construction of the stable manifold. We start on theeigenvectors very near the origin and integrate backwards in time.eps = 0.01;initset = 8eps * 81, 0<, eps * 8-1, 0<<;t0 = 0.0; h = -0.02; nsteps = 1000;plrange = 88-2.1, 2.1<, 8-2.1, 2.1<<; asprat = 1.0;rangeflag = True; ranger = plrange; bothdirflag = False;setcolor@8Black<D;8 manifold.nblabon = "Stable Manifold";arrowflag = True; arrowvec = 81 ê 2<;globestab = portrait[initset,t0,h,nsteps,1,2]-2-1012-2-1012xyStable ManifoldNow we carry this out for the unstable manifold. For this, we start near the equilibrium on the eigenvectors andintegrate forward in time.initset = 8eps * 80, 1<, eps * 80, -1<<;h = 0.02;labon = "Unstable Manifold";manifold.nb 9globeunstab = portrait[initset,t0,h,nsteps,1,2] -2-1012-2-1012xyUnstable ManifoldThe biggest surprise here is the homoclinic loop rooted in the origin: it is part of both the global unstable andglobal stable manifold! A little thought shows that this doesn't violate the rules, whether or not we imaginedsuch a thing could happen. If we start on the loop and go forward in time, we arrive at the origin, hence it is partof the stable manifold. If we start on the loop and go backward in time, we also arrive at the origin, so it is partof the unstable manifold. We can show the two manifolds together:labon = "Unstable and Stable Manifolds";10 manifold.nbshow@globestab, globeunstabD-2-1012-2-1012xyUnstable and Stable ManifoldsAlthough the homoclinic loop is not periodic, it is pretty clear that there is a periodic solution or solutions in itsinterior, and also an equilibrium point. The equilibrium point is easy to find -- it is {-1,-1}, also confirmed byDynPac:findpolyeq:8-1, -1<, 80, 0<, :H-1L1ê3, -H-1L2ê3>, :-H-1L2ê3, H-1L1ê3>>We classify the new equilibrium.classify2D@8-1, -1<Dstable HLL, indeterminate HNLL - centerLet's try a few initial conditions inside the loop. initset = 88-1, -1.4<, 8-1, -0.7<<;t0 = 0.0; h = 0.02; nsteps = 500;arrowvec = 81 ê 2<;labon = "Periodic Solutions";manifold.nb 11periodic = portrait[initset,t0,h,nsteps,1,2]