ME 406Example of a Phase Portrait with Multiple EquilibriasysidMathematica 4.1.2, DynPac 10.66, 3ê6ê2002plotreset; intreset;In this notebook, we construct a phase portrait for the sytem given below. (1)x° = -y - 1ÅÅÅÅ2x3+ 1ÅÅÅÅ2x4+ 3ÅÅÅÅ2y4- 3ÅÅÅÅ4y5 , y° = 4x - 6y3+ x4+3y4-x5 .We begin by defining the system for DynPac.setstate@8x, y<D; setparm@8<D;slopevec = 9-y -1ÅÅÅÅ2 x3+1ÅÅÅÅ2 x4+3ÅÅÅÅ2 y4-3ÅÅÅÅ4 y5,4x- 6y3+ x4+ 3y4- x5=;sysname = "Portrait";We use nfindpolyeq to look for equilibrium states.eqstates = nfindpolyeq881.7484, 2.<, 80.234465 + 1.3507 Â,2.<,80.234465 - 1.3507 Â,2.<, 8-1.21733, 2.<, 81., 1.75211<,80.383225 - 1.76687 Â, 0.743413 + 1.15521 Â<,80.383225 + 1.76687 Â, 0.743413 - 1.15521 Â<,81.76687 - 0.383225 Â, 1.15521 + 0.743413 Â<,81.76687 + 0.383225 Â, 1.15521 - 0.743413 Â<,8-0.383225 - 1.76687 Â, -0.743413 + 1.15521 Â<,8-0.383225 + 1.76687 Â, -0.743413 - 1.15521 Â<,8-1.76687 - 0.383225 Â, -1.15521 + 0.743413 Â<,8-1.76687 + 0.383225 Â, -1.15521 - 0.743413 Â<, 81., 1.17221<,8-0.234482 + 1.08109 Â, -0.706835 + 0.454868 Â<,8-0.234482 - 1.08109 Â, -0.706835 - 0.454868 Â<,81.08109 - 0.234482 Â, -0.454868 + 0.706835 Â<,81.08109 + 0.234482 Â, -0.454868 - 0.706835 Â<,80.234482 + 1.08109 Â, 0.706835 + 0.454868 Â<,80.234482 - 1.08109 Â, 0.706835 - 0.454868 Â<,8-1.08109 - 0.234482 Â, 0.454868 + 0.706835 Â<,8-1.08109 + 0.234482 Â, 0.454868 - 0.706835 Â<,81., -0.462164 + 0.659993 Â<, 81., -0.462164 - 0.659993 Â<, 80., 0.<<We have found five real equilibria. We name them.eq1 = eqstates@@1DD81.7484, 2.<eq2 = eqstates@@4DD8-1.21733, 2.<eq3 = eqstates@@5DD81., 1.75211<eq4 = eqstates@@14DD81., 1.17221<eq5 = eqstates@@25DD80., 0.<Now we classify all of these equilibria.classify2D@eq1Dunstable - saddleclassify2D@eq2Dunstable - saddleclassify2D@eq3Dunstable - nodeclassify2D@eq4Dunstable - saddleclassify2D@eq5Dstable HLL, indeterminate HNLL - centerWe see that four of the equilibria are unstable, and the fifth -- the origin -- has a stability which is not determinedby linearization. As we showed earlier, the stability of the equilibrium at the origin is established by a Liapunovfunction, namelyVlap = 4x2+ y2;The orbital derivative of Vlap is2 portrait.nbSimplify@orbdt@VlapDD-2 H-1 + xLH-2 + yLHx4+ 3y4Land this is clearly negative definite in a neighborhood of the origin. We begin the phase portrait by constructing the stable and unstable manifolds about the saddle points. Weset the plotting window:plrange = 88-4.1, 4.1<, 8-4.1, 4.1<<;We set range checking and integration in only one time direction.rangeflag = True; ranger = plrange; bothdirflag = False;eig1 = [email protected], -3.85689<,880.422889, -0.906182<, 8-0.793775, -0.608212<<<eigvec11 = First@[email protected], -0.906182<eigvec12 = Last@[email protected], -0.608212<Notice how different the eigenvalues are in magnitude. Motion out along the unstable manifold will be muchmore rapid than motion in along the stable manifold. We construct the unstable manifold first.∂∂∂∂= 0.005;initset = 8eq1 + ∂∂∂∂* eigvec11, eq1 - ∂∂∂∂* eigvec11<;t0 = 0.0; h = 0.002; nsteps = 5000;arrowflag = True; arrowvec = 81 ê 3<;labshift = 20;portrait.nb 3graph11 = portrait@initset, t0, h, nsteps, 1, 2D;-4-2 24x-4-224yPortraitNow the stable manifold.initset = 8eq1 + ∂∂∂∂* eigvec12, eq1 - ∂∂∂∂* eigvec12<;t0 = 0.0; h =-0.02; nsteps = 5000;arrowvec = 84 ê 5<;4 portrait.nbgraph12 = portrait@initset, t0, h, nsteps, 1, 2D;-4-2 24x-4-224yPortraitNow we look at the second equilibrium, eq2.eig2 = [email protected], -11.0901<,880.347401, -0.937717<, 8-0.927012, -0.375032<<<eigvec21 = First@[email protected], -0.937717<eigvec22 = Last@[email protected], -0.375032<First the unstable manifold.initset = 8eq2 + ∂∂∂∂* eigvec21, eq2 - ∂∂∂∂* eigvec21<;t0 = 0.0; h = 0.002; nsteps = 5000;portrait.nb 5arrowflag = True; arrowvec = 81 ê 2<;graph21 = portrait@initset, t0, h, nsteps, 1, 2D;-4-2 24x-4-224yPortraitNow the stable manifold.initset = 8eq2 + ∂∂∂∂* eigvec22, eq2 - ∂∂∂∂* eigvec22<;t0 = 0.0; h =-0.002; nsteps = 5000;6 portrait.nbgraph22 = portrait@initset, t0, h, nsteps, 1, 2D;-4-2 24x-4-224yPortraitNext we look at the last saddle, eq4.eig4 = [email protected], 1.21753<, 88-0.232615, 0.972569<, 80.910893, 0.412642<<<eigvec41 = First@[email protected], 0.972569<eigvec42 = Last@[email protected], 0.412642<First the unstable manifold.initset = 8eq4 + ∂∂∂∂* eigvec42, eq4 - ∂∂∂∂* eigvec42<;t0 = 0.0; h = 0.002; nsteps = 5000;arrowflag = True; arrowvec = 81 ê 2<;portrait.nb 7graph42 = portrait@initset, t0, h, nsteps, 1, 2D;-4-2 24x-4-224yPortraitNow the stable manifold.bothdirflag = False;initset = 8eq4 + ∂∂∂∂* eigvec41, eq4 - ∂∂∂∂* eigvec41<;t0 = 0.0; h =-0.002; nsteps = 85000, 5000<;arrowvec = 89 ê 10<;8 portrait.nbgraph41 = portrait@initset, t0, h, nsteps, 1, 2D;-4-2 24x-4-224yPortraitBefore looking at the unstable node eq3 and the stable spiral at the origin, we first combine the graphs wehave so far.portrait.nb 9firstbigraph =show@graph11, graph12, graph21, graph22 , graph41, graph42D;-4-2 24x-4-224yPortraitNow we look at the unstable node eq3. Some of the integral curves leaving this point have already beenconstructed. A little experimenting shows that one more curve is useful.eig3 = [email protected], 2.22907<, 880.499355, -0.866398<, 8-0.920325, 0.391154<<<eigvec31 = First@[email protected], -0.866398<eigvec32 = Last@[email protected], 0.391154<initset = 8eq3 + ∂∂∂∂* eigvec32<;t0 = 0.0; h = 0.002; nsteps = 5000;arrowvec = 81 ê 2<;10 portrait.nbgraph32 = portrait@initset, t0, h, nsteps, 1, 2D;-4-2 24x-4-224yPortraitportrait.nb 11bigraph = show@firstbigraph, graph32D;-4-2 24x-4-224yPortrait We add a few more orbits to fill in the picture. initset = 883, 1.75<, 82, 4<, 82, -3<, 8-3, 0<, 8-4, 2<, 8-2, -2<<;t0 = 0.0; h = 0.00025; nsteps = 5000;bothdirflag = True;arrowvec = 83 ê 5<;12 portrait.nblastpiece = portrait@initset, t0, h, nsteps, 1, 2D;-4-2 24x-4-224yPortraitNow we put the pictures together.portrait.nb 13newbigraph = show@bigraph, lastpieceD;-4-2 24x-4-224yPortraitOne more orbit.bothdirflag = True; init = 80, 1.25<;h = 0.00025; nsteps = 4000; t0 = 0.0;solast = integrate@init, t0, h, nstepsD;14 portrait.nblastpiece = phaser@solastD;-4-2 24x-4-224yPortraitportrait.nb 15newbigraph2 = show@newbigraph, lastpieceD;-4-2