ME 406Trapping a Limit CyclesysidMathematica 6.0.3, DynPac 11.01, 1ê13ê2009intreset;plotreset; imsize = 250;üIntroductionIn this notebook, we localize a limit cycle by orbit trapping. The example we consider, given below, is taken from Nonlinear Differential Equations and Dynamical Systems, 2nd edition, Ferdinand Verhulst, Springer, 1996, example 4.8, page 48, with a correction of a sign error. x°= y - x Ix2+ y2- 2 x - 3M , y°= -x - y Ix2+ y2- 2 x - 3M .We now define the system for Mathematica. We include a parameter a so that variations of the system can be studied easily. For our system, a = 1.setstate@8x, y<D;slopevec = 8y - a * x * Hx ^ 2 + y^2 - 2 * x - 3L, -x - a * y * Hx^2 + y ^ 2 - 2 * x - 3L<;setparm@8a<D;parmval = 81<;sysname = "LimCyc 4.8";plrange = 88-5, 5<, 8-5, 5<<;üEquilibrium PointWe look for equilibrium states using findpolyeq.findpolyeq:80, 0<, :-Â - 3 a2 a,1 - 3 Â a2 a>, :Â - 3 a2 a,1 + 3 Â a2 a>>Mathematica has found three equilibria of which only one is real. As shown in class, there are no other real equilibria. (Important caution: we never know for sure without independent confirmation whether a root-finder such as findpolyeq has found all of the equilibrium states.) Let's look at the stability of the equilibrium at the origin:orbtrap.nb 1Mathematica has found three equilibria of which only one is real. As shown in class, there are no other real equilibria. (Important caution: we never know for sure without independent confirmation whether a root-finder such as findpolyeq has found all of the equilibrium states.) Let's look at the stability of the equilibrium at the origin:classify2D@80, 0<DAbbreviations used in classify2D.L = linear, NL = nonlinear, R2 = repeated root.Z1 = one zero root, Z2 = two zero roots.This message printed once.unstable - spiralThis is typical of what we would find inside a limit cycle. In a search for periodic solutions, we move next to the Bendixson test.üApplication of Bendixson TestFrom the Bendixson test, we know that periodic solutions can occur only in regions where the divergence of the slope vector changes sign or is identically zero. We use the function signcontour to map the algebraic sign of the divergence here. The function dival returns the divergence with current values of the parameters substituted.bendix = dival@slopevecD6 + 6 x - 4 x2- 4 y2We turn off the graphics axes and turn on the frame.axon = False; frameon = True;orbtrap.nb 2graph1 = signcontour@bendixDThe divergence is positive in the white region (which can be shown to be a circle with center at x = 3/4, y = 0 and radius 33/4). By the Bendixson criterion, any periodic solution will have to enclose both white and grey regions. This is only a necessary condition, and we now move to orbit trapping techniques to establish that there is a periodic solution.üOrbit TrappingBecause they are easy, we will start with circles. We define our trapping curve as a circle of radius b with center at the origin. We first define an arc and then a curve from the arc.trapperarc@b_D := 88b * Cos@uD, b * Sin@uD<, 8u, 0, 2 * Pi<<;trappercurve@b_D := 8trapperarc@bD<Now we use the routine orbcross to see whether a given circle traps orbits. We start with a radius b = 0.5.orbcross@[email protected]: No.Crossings in negative N direction: No.Crossings in positive N direction: Yes.Thus the circle with radius 0.5 traps orbits outside. If we can find a larger circle which traps orbits inside, we will have established the existence of a periodic solution. We try b = 5.0.orbtrap.nb 3orbcross@[email protected]: No.Crossings in negative N direction: Yes.Crossings in positive N direction: No.We have been lucky. We now know that there is a periodic solution somewhere between b = 0.5 and b = 5.0. Let's refine our guesses.orbcross@[email protected]: Yes.Crossings in negative N direction: No.Crossings in positive N direction: Yes.This is optimum because it still traps, but there are tangencies. For a slightly larger b, we would get crossings in both directions. Let's verify that:orbcross@[email protected]: No.Crossings in negative N direction: Yes.Crossings in positive N direction: Yes.Now we locate more precisely an outer trapping curve, working back from b = 5.0. orbcross@[email protected]: No.Crossings in negative N direction: Yes.Crossings in positive N direction: No.orbcross@[email protected]: Yes.Crossings in negative N direction: Yes.Crossings in positive N direction: No.This also appears optimum, and we verify it by a slightly smaller b:orbcross@[email protected]: No.Crossings in negative N direction: Yes.Crossings in positive N direction: Yes.Now we have a periodic orbit trapped between b = 1.0 and b = 3.0. Let's plot these two curves in red:orbtrap.nb 4setcolor@8Red<D;graph2 = plotcurve@[email protected];graph3 = plotcurve@[email protected];graph4 = show@graph1, graph3, graph2DNow let's add the limit cycle in blue.nsteps = 500;h = 0.025;t0 = 0;initvec = 82, 0<;sol1 = limcyc@initvec, t0, h, nstepsD;arrowflag = True;arrowvec = 81 ê 2<;setcolor@8Blue<D;orbtrap.nb 5graph5 = phaser@sol1D-4-2024-4-2024xyLimCyc 4.8 8a< = 81<Now we put it all together:graph6 = show@graph4, graph5DAs a final calculation, we add a few orbits approaching the limit cycle, also in blue.orbtrap.nb 6initset = 88-2, 0<, 80, 2<, 83, 3<, 8-3, -3<, 83, -3<, 8-3, 3<, 81, 1<, 81, -1<<;nsteps = 100;arrowvec = 81 ê 8<;graph7 = portrait@initset, t0, h, nsteps, 1, 2D-4-2024-4-2024xyLimCyc 4.8 8a< = 81<orbtrap.nb 7graph8 = show@graph6, graph7DThe approach orbits show clearly the stability of the limit cycle.orbtrap.nb