ME 406Examples of Periodic SolutionsIn[403]:=sysidMathematica 6.0.3, DynPac 11.02, 2ê6ê2009In[404]:=plotreset; intreset; imsize = 250;üIntroductionIn this notebook we look at several examples of periodic solutions of autonomous systems of two equations. We visualize the solutions in the phase plane, where periodic solutions have a distinctive topological signature: their orbits are simple closed curves. As we shall see from the examples, there are a variety of possible outcomes from the search for periodic solutions. For some systems there are no periodic solutions. For others, every solution is periodic. For still others, there are periodic solutions in some regions of the phase plane but not everywhere. Finally there is the surprising and interesting phenomenon of the limit cycle, which is an isolated periodic solution. Our examples will illustrate all of this behavior. In later notebooks we will cover some of the tools used here in a more systematic search for periodic solutions. üExample 1: Linear System with All Solutions PeriodicWe have already used the linear oscillator a number of times as an example, and we know that the undamped oscillator has periodic solutions. We repeat that example here. In scaled form the equations are x° = y, y°= -x .It is easy to solve the initial value problem for this system analytically. If we let x(0) = a, and x°(0) = b, we get the solutionx(t) = a cos(t) + b sin(t), y(t) = -a sin(t) + b cos(t) .Thus all solutions are periodic with the same period, namely 2p. From the above solution, we find by direct calculation that [x(t)D2 + [y(t)D2 = a2 + b2, so that all of the orbits are circles. Of course if we didn't know any of this, we could find a numerical approximation to the truth by solving the system numerically and plotting the orbits. Let's do that. We define the system for DynPac.In[405]:=setstate@8x, y<D; setparm@8<D; slopevec = 8y, -x<;In[406]:=sysname = "LinOsc";We now construct a single solution, with initial point {1,0}.perexamp.nb 1We now construct a single solution, with initial point {1,0}.In[407]:=t0 = 0.0; h = 0.02; nsteps = 350; initvec = 81, 0<;In[408]:=sol1 = integrate@initvec, t0, h, nstepsD;In[409]:=plrange = 88-2, 2<, 8-2, 2<<; asprat = 1.0;In[410]:=arrowflag = True; arrowvec = 81 ê3, 3 ê4<;In[411]:=graph1 = phaser@sol1DOut[411]=-2-112x-2-112yLinOscWe see that the orbit is a simple closed curve. How do we find the period of this numerical solution? We use the function period, applied to the solution.In[412]:=period@sol1DOut[412]=6.28This function returns an answer accurate to within one time step, and we see that the result is consistent with the exact answer of 2p. We construct several other solutions, find their periods and then plot them together.perexamp.nb 2In[413]:=sol2 = [email protected], 0<, t0, h, nstepsD;In[414]:=period@sol2DOut[414]=6.28In[415]:=sol3 = [email protected], 0<, t0, h, nstepsD;In[416]:=period@sol3DOut[416]=6.28In[417]:=sol4 = [email protected], 0<, t0, h, nstepsD;In[418]:=period@sol4DOut[418]=6.28In[419]:=phaser@8sol1, sol2, sol3, sol4<DOut[419]=-2-112x-2-112yLinOscThe numerical work we have just finished suggests that all solutions are periodic and that they have the same period, 6.28. Although we can never prove such results with numerical work, the numerical techniques are our primary exploratory tools, and at the very least suggest what properties we might want to establish by analysis. perexamp.nb 3The numerical work we have just finished suggests that all solutions are periodic and that they have the same period, 6.28. Although we can never prove such results with numerical work, the numerical techniques are our primary exploratory tools, and at the very least suggest what properties we might want to establish by analysis. It is worth noting that this system, which has all solutions periodic, is immediately converted to a system with no periodic solutions if we have any amount of damping, no matter how small. There is an important concept called structural stability behind that statement. We will discuss it later in the course.üExample 2: Nonlinear System with All Solutions PeriodicFor our next example, we consider a nonlinear oscillator governed by the well-known Duffing's equation. We consider the simplest case of no damping. The equation is x– + x + x3 = 0 .A physical model for this system is an oscillator with a nonlinear spring, for which the restoring spring force is given by x + x3. In vibration theory such a spring is called a hard spring, because the spring force per unit extension, 1 + 3x2, increases as the displacement increases. We convert the equation to a system. x° = y , y° = -x -x3. It is actually possible to solve this equation analytically in terms of elliptic functions. Although we won't carry this out, we will look at the first step. The second order equation is of a form known as a conservative equation. It admits a conservation law which is an energy equation. The general form of a conservation equation isx– +f(x) = 0 .If we multiply through by x°, we can integrate the result to get12x°2 + V(x) = constant ,where the potential V is given by Ÿf HxL „ x . For the Duffing equation, the result is12x°2 + 12x2 + 14x4 = 12x02+ 14x04 ,where we have parametrized the solution by initial conditions x(0) = x0, x°(0) = 0. We may solve for x° in terms of x:x° = Ix02- x2M I1 +12Ax02+ x2EM .The above equation is separable and can be further manipulated to generate the solution in terms of elliptic functions. We will use it in a more limited way shortly to derive an expression for the period of the motion.We define this system for DynPac, generate four solutions, plot them to discover that they are periodic, and then calculate the periods.In[420]:=setstate@8x, y<D; setparm@8<D; slopevec = 9y, -x - x3=;perexamp.nb 4In[421]:=sysname = "Duffing";We generate a solution with x0= 1 and then plot it.In[422]:=sol1 = integrate@81, 0<, 0.0, 0.02, 400D;In[423]:=period@sol1DOut[423]=4.76In[424]:=arrowvec = 81 ê3<;In[425]:=plrange = 88-3, 3<, 8-3, 3<<;In[426]:=phaser@sol1DOut[426]=-3-2-1123x-3-2-1123yDuffingWe see a closed curve for the periodic orbit, but it is not circular. In the limit of small amplitude the cubic nonlinearity in the force should be unimportant, the orbit should be close to circular, and the period should be close to 2p. Let's check that by doing a second integration for x0= 0.25.In[427]:=sol2 = [email protected],