ME 406Example of a Phase Portrait with Multiple EquilibriasysidMathematica 6.0.3, DynPac 11.01, 1ê13ê2009plotreset; intreset; imsize = 250;In this notebook, we construct a phase portrait for the sytem given below. (1)x° = -y - 12x3+ 12x4+ 32y4- 34y5 , y° = 4x - 6y3+ x4+3y4-x5 .We begin by defining the system for DynPac.setstate@8x, y<D; setparm@8<D;slopevec = :-y -12 x3+12 x4+32 y4-34 y5, 4 x - 6 y3+ x4+ 3 y4- x5>;sysname = "Portrait";We use nfindpolyeq to look for equilibrium states.eqstates = nfindpolyeq88-1.21733, 2.<, 8-0.383225 - 1.76687 Â, -0.743413 + 1.15521 Â<,8-0.383225 + 1.76687 Â, -0.743413 - 1.15521 Â<,80.383225 - 1.76687 Â, 0.743413 + 1.15521 Â<,80.383225 + 1.76687 Â, 0.743413 - 1.15521 Â<,80.234465 - 1.3507 Â, 2.<, 80.234465 + 1.3507 Â, 2.<,8-1.08109 - 0.234482 Â, 0.454868 + 0.706835 Â<,8-1.08109 + 0.234482 Â, 0.454868 - 0.706835 Â<,81.08109 - 0.234482 Â, -0.454868 + 0.706835 Â<,81.08109 + 0.234482 Â, -0.454868 - 0.706835 Â<,81., -0.462164 + 0.659993 Â<, 81., -0.462164 - 0.659993 Â<,81.76687 - 0.383225 Â, 1.15521 + 0.743413 Â<,81.76687 + 0.383225 Â, 1.15521 - 0.743413 Â<,8-1.76687 - 0.383225 Â, -1.15521 + 0.743413 Â<,8-1.76687 + 0.383225 Â, -1.15521 - 0.743413 Â<,80.234482 - 1.08109 Â, 0.706835 - 0.454868 Â<,80.234482 + 1.08109 Â, 0.706835 + 0.454868 Â<,8-0.234482 - 1.08109 Â, -0.706835 - 0.454868 Â<,8-0.234482 + 1.08109 Â, -0.706835 + 0.454868 Â<,81.7484, 2.<, 81., 1.75211<, 81., 1.17221<, 80., 0.<<Length@eqstatesD25We have found five real equilibria. We name them.eq1 = eqstates@@22DD81.7484, 2.<eq2 = eqstates@@1DD8-1.21733, 2.<eq3 = eqstates@@23DD81., 1.75211<eq4 = eqstates@@24DD81., 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, namely2 portrait.nbVlap = 4 x2+ y2;The orbital derivative of Vlap isSimplify@orbdt@VlapDD-2 H-1 + xL H-2 + yL Ix4+ 3 y4Mand 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.We set 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<;portrait.nb 3graph11 = portrait@initset, t0, h, nsteps, 1, 2D-4-224x-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-224x-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;arrowflag = True; arrowvec = 81 ê2<;portrait.nb 5graph21 = portrait@initset, t0, h, nsteps, 1, 2D-4-224x-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-224x-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-224x-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-224x-4-224yPortraitBefore looking at the unstable node eq3 and the stable spiral at the origin, we first combine the graphswe have so far.portrait.nb 9firstbigraph =show@graph11, graph12, graph21, graph22, graph41, graph42D-4-224x-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-224x-4-224yPortraitbigraph = show@firstbigraph, graph32D-4-224x-4-224yPortrait We add a few more orbits to fill in the picture. portrait.nb 11initset = 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<;lastpiece = portrait@initset, t0, h, nsteps, 1, 2D-4-224x-4-224yPortraitNow we put the pictures together.newbigraph = show@bigraph, lastpieceD;One more orbit.bothdirflag = True; init = 80, 1.25<; h = 0.00025; nsteps = 4000; t0 = 0.0;solast = integrate@init, t0, h, nstepsD;12 portrait.nblastpiece = phaser@solastD-4-224x-4-224yPortraitnewbigraph2 = show@newbigraph, lastpieceD;As a final touch, we place a blue dot on the stable equilibrium and red dots on the saddles, and a green dot on theunstable node.setcolor@8Red, Blue, Green<D;ptsize = 0.03;dotgraph1 = dots@8eq1, eq2, eq4<D;dotgraph2 = dots@8eq5<D;dotgraph3 = dots@8eq3<D;imsize
View Full Document