Mecflu - caso de oscila??o n?o-linear

Mecflu - caso de oscila??o n?o-linear

Nonlinear Oscillations and Chaos


d d

(a) (b) (c) m ms x

The unextended length of each spring is , as shown in (a). In order to attach the mass m, each spring must be stretched a distance d, as indicated in (b). When the mass is moved a distance x, as in (c), the force acting on the mass (neglecting gravity) is

Then, kx kx Fx x x d

d x kx kx

128 CHAPTER 4 Expanding the radical in powers of 22xA and retaining only the first two terms, we have dx Fx kx d kx

kdkd x

The potential is given by

4-2. Using the general procedure explained in Section 4.3, the phase diagram is constructed as follows:

x E


diagram is given in (b):


4-4. Differentiation of Rayleigh’s equation above yields

implies that ya x ya x

ya x

When these expressions are substituted in (1), we find y yaa ba a a

by b b y y by

Multiplying by 0 3b y a and rearranging, we arrive at van der Pol’s equation:


tan x

Ox π⁄ 2 x + x +1

The procedure is to use this approximate solution as a starting point and to substitute 138xπ=

value, then use the result as a new starting point and repeat the calculation. This procedure leads to the following values:

Thus, the solution is x = 1.3319. Parts b) and c) are solved in exactly the same way with the results:

b) x = 1.9151 c) x = 0.9271

4-6. For the plane pendulum, the potential energy is

If the total energy is larger than , all values of θ are allowed, and the pendulum revolves continuously in a circular path. The potential energy as a function of θ is shown in (a) below. 2mgA

(a) Since T = E – U(θ), we can write


and, therefore, the phase paths are constructed by plotting

versus θ. The phase diagram is shown in (b) below. E = 2mg E = 3mg


4-7. Let us start with the equation of motion for the simple pendulum:

where 2g≡Aω. Put this in terms of the horizontal component by setting sinyx≡=Aθ.

Solving for θ and taking time derivatives, we obtain 2

Since we are keeping terms to third order, we need to get a better handle on the term. Help comes from the conservation of energy:

where 0θ is the maximum angle the pendulum makes, and serves as a convenient parameter that describes the total energy. When written in terms of , the above equation becomes (with the obvious definition for ) y 0y

Substituting (4) into (2), and the result into (1) gives

Using the binomial expansion of the square roots and keeping terms up to third order, we can obtain for the x equation of motion

4-8. For x > 0, the equation of motion is

Thus, the phase path is a parabola with a vertex on the x-axis at x = A and symmetrical about both axes as shown below.

t = 0

A x

Because of the symmetry, the period τ is equal to 4 times the time required to move from x = A to x = 0 (see diagram). Therefore, from (2) we have

4-9. The proposed force derives from a potential of the form kx x a Ux

which is plotted in (a) below.


For small deviations from the equilibrium position (x = 0), the motion is just that of a harmonic oscillator.

For energies , the particle cannot reach regions with x < –a, but it can reach regions of x > a if . For the possibility exists that the particle can be trapped near x = a. 6E E<

A phase diagram for the system is shown in (b) below.


4-10. The system of equations that we need to solve are

The values of ω that give chaotic orbits are 0.6 and 0.7. Although we may appear to have chaos for other values, construction of a Poincaré plot that samples at the forcing frequency show that they all settle on a one period per drive cycle orbit. This occurs faster for some values of ω than others. In particular, when ω = 0.8 the plot looks chaotic until it locks on to the point . The phase plot for ω = 0.3 shown in the figure was produced by numerical integration of the system of equations (1) with 100 points per drive cycle. The box encloses the point on the trajectory of the system at the start of a drive cycle. In addition, we also show

Poincaré plot for the case ω = 0.6 in figure, integrated over 8000 drive cycles with 100 points per cycle.

4-1. The three-cycle does indeed occur where indicated in the problem, and does turn chaotic near the 80th iteration. This value is approximate, however, and depends on the precision at which the calculations are performed. The behavior returns to a three-cycle near the 200th iteration, and stays that way until approximately the 270th iteration, although some may see it continue past the 300th.



x = 0.4 x = 0.75

These plots are created in the manner described in the text. They are created with the logistic equation

The first plot has the seed value as asked for in the text. Only one additional seed has been done here () as it is assumed that the reader could easily produce more of these plots after this small amount of practice.

136 CHAPTER 4 4-13.

x = 0.7 x = 0.700000001 iteration x = 0.7 x = 0.7000000001 iteration

The plots are created by iteration on the initial values of (i) 0.7, (i) 0.700000001, and (ii) 0.7000000001, using the equation

A subset of the iterates from (i) and (i) are plotted together, and clearly diverge by n = 39. The plot of (i) and (i) clearly diverge by n = 43.


x = 0.9 x = 0.9000001 fractional difference iteration


The fractional difference is defined as xyx−, and clearly exceeds 30% when n = 30.

4-15. A good way to start finding the bifurcations of the function f(α,x) = α sin πx is to plot its bifurcation diagram.

One can expand regions of the diagram to give a rough estimate of the location of a bifurcation. Its accuracy is limited by the fact that the map does not converge very rapidly near the bifurcation point, or more precisely, the Lyapunov exponent approaches zero. One may continue undaunted, however, with the help of a graphical fractal generating software application, to estimate quite a few of the period doublings nα.Using Fractint for Windows, and Equation (4.47) to compute the Feigenbaum constant, we can obtain the following:

One can see that although we should obtain a better value of δ as n increases, numerical precision and human error quickly degrade the quality of the calculation. This is a perfectly acceptable answer to this question.

One may compute the nα to higher accuracy by other means, all of which are a great deal more complicated. See, for example, Exploring Mathematics with Mathematica, which exploits the vanishing Lyapunov exponent. Using their algorithm, one obtains the following:

Note that these are shown here only as reference, and the student may not necessarily be expected to perform to this degree of sophistication. The above values are only good to about

, but this time only limited by machine precision. Another alternative in computing the

Feigenbaum constant, which is not requested in the problem, is to use the so-called ”supercycles,” or super-stable points , which are defined by

The values obey the same scaling as the bifurcation points, and are much easier to compute since these points converge faster than for other α (the Lyapunov exponent goes to –∞). See, for example, Deterministic Chaos: An Introduction by Heinz Georg Schuster or Chaos and Fractals: New

Frontiers of Science by Peitgen, Jürgens and Saupe. As a result, the estimates for δ obtained in this way are more accurate than those obtained by calculating the bifurcation points.

4-16. The function y = f(x) intersects the line y = x at 0xx=, i.e. is defined as the point where 0x

Now define 0nnxxε≡−. If we have very close to , then 1x0x1ε should be very small, and we

If the approximation (1) remains valid from the initial value, we have 1nεβε+ .

b) Clearly, when 1β< we have stability since

Similarly we have a divergent sequence when 1β>, although it will not really be exponentially divergent since the approximation (1) becomes invalid after some number of iterations, and normally the range of allowable is restricted to some subset of the real numbers. nx


α = 0.4 α = 0.7iteration

The first plot (with α = 0.4) converges rather rapidly to zero, but the second (with α = 0.7) does appear to be chaotic.


The tent map always converges to zero for α < 0.5. Near α = 0.5 it takes longer to converge, and that is the artifact seen in the figure. There exists a “hole” in the region 0.5 < α < 0.7 (0.7 is approximate), where the iterations are chaotic but oscillate between an upper and lower range of values. For α > 0.7, there is only a single range of chaos, which becomes larger until it fills the range (0,1) at α = 1.

4-19. From the definition in Equation (4.52) the Lyapunov exponent is given by

1 lim ln

i xdf nd x

The tent map is defined as

2for 0 21for 121

This gives 2dfdxα=, so we have

As indicated in the discussion below Equation (4.52), chaos occurs when λ is positive: 12α> for the tent map.


xy 4-21.

The shape of this plot (the attractor) is nearly identical to that obtained in the previous problem. In Problem 4-20, however, we can clearly see the first few iterations (0,0), (1,0), (–0.4,0.3),

NONLINEAR OSCILLATIONS AND CHAOS 141 whereas the next iteration (1.076,–0.12) is almost on the attractor. In this problem the initial value is taken to be on the attractor already, so we do not see any transient points.

4-2. The following. system of differential equations were integrated numerically

using different values of B in the range [9.8,13.4], and with a variety of initial conditions. The integration range is over a large number of drive cycles, throwing away the first several before starting to store the data in order to reduce the effects of the transient response. For the case B = 9.8, we have a one period per three drive cycle orbit. The phase space plot (line) and Poincaré section (boxes) for this case are overlaid and shown in figure (a). All integrations are done here with 100 points per drive cycle. One can experiment with B and determine that the system becomes chaotic somewhere between 9.8 and 9.9. The section for B = 10.0, created by integrating over 8000 drive cycles, is shown in figure (b). If one further experiments with different values of B, and one is also lucky enough to have the right initial conditions, (0,0) is one that works, then a transition will be found for B in the range (1.6,1.7). As an example of the different results one can get depending on the initial conditions, we show two plots in figure (c). One is a phase plot, overlaid with its section, for B = 12.0 and the initial condition (0,0). Examination of the time evolution reveals that it has one period per cycle. The second plot is a Poincaré section for the same B but with the initial condition (10,0), clearly showing chaotic motion. Note that the section looks quite similar to the one for B = 10.0. Another transition is in the range (13.3,13.4), where the orbits become regular again, with one period per drive cycle, regardless of initial conditions. The phase plot for B = 13.4 looks similar to the one with B = 12.0 and initial condition (0,0).

To summarize, we may enumerate the above transition points by B, , and .

Circumventing the actual task of computing where these transition points are, we do know that , 1, and 13.3

We should remind ourselves, though, that the above list only applies for B in the range we have examined here. We do not know the behavior when B < 9.8 and B > 13.4, without going beyond the scope of this problem.

y (a)

y (b)

y (c) y (d)


The results one should get from doing this problem should be some subset of the results shown in figures (a), (b), and (c) (for K = 0.8, 3.2, and 6.4, respectively). These were actually generated using some not-so-random initial points so that a reasonably complete picture could be made. What look to be phase paths in the figures are actually just different points that come from iterating on a single initial condition. For example, in figure (a), an ellipse about the origin (just pick one) comes from iterating on any one of the points on it. Above the ellipses is chaotic orbit, then a five ellipse orbit (all five come from a single initial condition), etc. The case for K = 3.2 is similar except that there is an orbit outside of which the system is always undergoing chaotic motion. Finally, for K = 6.4 the entire space is filled with chaotic orbits, with the exception of two small lobes. Inside of these lobes are regular orbits (the ones in the left are separate from the ones in the right).

q⁄π p ⁄π (a)

q⁄π p ⁄π (b)

q⁄π p ⁄π (c)

4-24. a) The Van de Pol equation is

dx dx xa x dt dt ωµ+=−

Now look for solution in the form

dx du bt dt dt ωω=−+ and dx d u bt

dt dt ωω=−+

Putting these into the Van de Pol equation, we obtain

() () ( ) cos ( ) 2 ( ) cos sin du t du t ut b t u t bu t t a b t

From this one can see that u(t) is of orderµ (i.e. ~()uOµ), which is assumed to be small here. Keeping only terms up to order µ, the above equation reads

sin sin 3 4 du t ut b t t a b t b ba t t

0t (where we have used the identity 2 0004sincossinsin3tttωωω=+ω )

This equation has 2 frequencies (0ω and 30ω), and is complicated. However, if then the term 2b = a

0sintω disappears and the above equation becomes


ba ttµµutωωωω=−=−

So, finally putting this form of u(t) into (1), we obtain one of the exact solutions of Van de Pol equation:

a utattµωωω=−

b) See phase diagram below. Since 0.05µ=is very small, then actually the second term in the expression of u(t) is negligible, and the phase diagram is very close to a circle of radius b = 2a = 2.

4-25. We have used Mathematica to numerically solve and plot the phase diagram for the van de Pol equation. Because 0.07µ= is a very small value, the limit cycle is very close to a circle of radius b = 2a = 2.

a) In this case, see figure a), the phase diagram starts at the point (x = 1, x′ = 0) inside the limit cycle, so the phase diagram spirals outward to ultimately approach the stable solution presented by the limit cycle (see problem 4-24 for exact expression of stable solution).

b) In this case, see figure b), the phase diagram starts at the point (x = 3, x′ = 0) outside the limit cycle, so the phase diagram spirals inward to ultimately approaches the stable solution presented by the limit cycle (see problem 4-24 for exact expression of stable solution).

4-26. We have used Mathematica to numerically solve and plot the phase diagram for the van de Pol equation. Because 0.5µ= is not a small value, the limit cycle is NOT close to a circle (see problem 4-24 above).

a) In this case, see figure a), the phase diagram starts at the point (x = 1, x′ = 0) inside the limit cycle, so the phase diagram spirals outward to ultimately approach the stable solution presented by the limit cycle (see problem 4-24 for exact expression of stable solution).

NONLINEAR OSCILLATIONS AND CHAOS 147 b) In this case (see figure below), the phase diagram starts at the point (x = 3, x′ = 0) outside the limit cycle, so the phase diagram spirals inward to ultimately approaches the stable solution presented by the limit cycle (see problem 4-24 for exact expression of stable solution).