Você está na página 1de 28

# PY 502, Computational Physics (Fall 2012) Numerical Solutions of Classical Equations of Motion

## Anders W. Sandvik, Department of Physics, Boston University

Introduction

Classical equations of motion, i.e., Newtons laws, govern the dynamics of systems ranging from very large, such as solar systems and galaxies, to very small, such as molecules in liquids and gases (in cases where quantum mechanical uctuations can be neglected, which is often the case). In between these extremes, Newtons equations of motion apply, litterally, to everything that moves. Exact analytical solutions of the equations of motion exist only for simple systems, of the types that are discussed in elementary classical mechanics courses, and therefore numerical integration methods are very important in practice. Here we will discuss some commonly used dierential equation solvers and use them to study the dynamics of mechanical systems, including ones that exhibit chaotic dynamics. We will limit the discussion to a single moving body, although the methods can be easily generalized to many-body systems as welldynamics of many-body systems will be discussed later in connection with molecular dynamics simulations. While some of the numerical schemes that we will discuss here are particularly suitable for integrating classical equations of motions, we will also described methods, such as the classic Runge-Kutta algorithm, that are more generally applicable to a large class of ordinary dierential equations. The discussion of systems with chaotic dynamics, although here introduced in the context of classical mechanics, is also of relevance more broadly in studies of nonlinear dynamics.

## Basic algorithms for equations of motion

Consider a single object (here regarded as a point particle) with mass m moving in one dimension. With its time-dependent position denoted x(t), the dierential equation governing its dynamics is x (t) = 1 F [x(t), x (t), t], m (1)

where F is the total force acting on the particle, and x and x are the rst and second time derivatives of x. We have indicated that the force may depend on x, x and t. These dependencies typically come from from a position dependent static potential, a velocity dependent damping (friction), and a time-dependent driving force, but there are other natural posibilities as well, e.g., a position dependent friction. To study the system numerically, it is convenient to rewrite the second-order dierential equation (1) as two coupled rst-order equations. Giving the velocity its standard symbol v (t), we have x (t) = v (t) 1 F [x(t), v (t), t]. v (t) = m 1

(2)

## 1.5 1 0.5 0 -0.5 -1 -1.5 0 10 20 30 40 50 0.5 0 10 20 30 40 50 0.8

t=0.01 t=0.001

E(t) t

x(t)

0.7

0.6

Figure 1: Time dependence of the position x and the energy E = (1/2)kx2 +(1/2)mv 2 of a harmonic oscillator with k = m = 1 (which gives an oscillation period 2 ) integrated using the Euler method with two dierent time steps; t = 102 and 103 . To integrate this set of equations, we discretize the time-axis as t = t0 , t1 , . . ., with a constant time step tn+1 tn = t . The initial values x0 = x(t0 ) and v0 = v (t0 ) are used to start the integration process.

2.1

Euler algorithm

The simplest way to advance the time from tn to tn+1 is to use the rst-order approximation; xn+1 = xn + t vn , vn+1 = vn + t an , (3)

where an = F (xn , vn , tn )/m is the acceleration. Clearly, the error made in each step of this algorithm is O(2 t ). The method, which is called Eulers forward method, is in general not very useful in practice. For example, in systems with no damping or driving force, the energy should be conserved. However, with the Euler method the energy typically diverges with time, whereas in most higher-order methods the energy errors are bounded. Fig. 1 shows results obtained with the Euler method for the energy and the position of a harmonic oscillator with k = m = 1 (F = kx). Even for a time step as as small as 103 , the energy error is as large as 5% at t = 50 (corresponding to less than 9 oscillations); it grows faster than linearly with t. There are several ways to proceed to derive more accurate, higher-order integration schemes; we will here discuss manipulations leading to a few practically useful formulas.

2.2

## Leapfrog and Verlet algorithms

To simplify the discussion initially, we will here rst assume that there is no damping, i.e., the force and the acceleration are velocity independent. We begin by writing the Taylor expansion of 2

initial values v(-1/2) x(0) v(1/2) x(1) v(3/2) x(2) v(5/2) x(3)

Figure 2: Position and velocity grids used in the leapfrog method. The position is calculated at integer multiples of the time step, tn = nt , n = 0, 1, . . ., and the velocity is evaluated at the times tn1/2 = (n 1/2)t between these points. The integration starts with given n = 0 values. xn+1 = x(tn + t ) to second order in t ; 1 a(xn , tn ) + O(3 x(tn+1 ) = x(tn ) + t v (tn ) + 2 t ). 2 t Noting that v (tn ) + (t /2)a(xn , tn ) = v (tn+1/2 ) to order 2 t , we can rewrite this as x(tn+1 ) = x(tn ) + t v (tn + t /2) + O(3 t ). (5) (4)

We thus need a formula that propagates the velocity on a time grid with points tn+1/2 = tn + t /2, i.e., between the integer-labeled time points tn = t0 + nt used for the position. If we use the rstorder expansion of the velocity, v (tn+1/2 ) = v (tn1/2 ) + t a(tn1/2 ) + O(2 t ), the error remains O(3 ) in Eq. (5), but we have a problem since this requires the acceleration, and hence the position t x (on which the force depends), on the grid points tn+1/2 used only for the velocity. However, it appears intuitively clear, by symmetry, that we should actually use the acceleration at tn+1 to propagate the velocity from tn+1/2 to tn+3/2 , i.e., v (tn+3/2 ) = v (tn+1/2 ) + t a(xn+1 , tn+1 ). The scheme combining Eqs. (5) and (6) is often called the leapfrog method; vn+1/2 = vn1/2 + t an , xn+1 = xn + t vn+1/2 . The leapfrog grid is illustrated in Fig. 2. We normally have initial conditions in the form of x0 and v0 . In order to start the leapfrog method we need v1/2 , which we can get, up to an error 2 t , using v1/2 = v0 t a0 + O(2 t ). 2 (8) (7) (6)

3 3 The 2 t error may seem to spoil the O (t ) error scaling of the leapfrog method. However, the O (t ) error is a single-step error; we will see in Sec. 2.3 that the long-time error, i.e., after integrating a large number of steps, is O(2 t ), and hence the error in (8) is of appropriate order.

It turns out that the position x in the leapfrog method in fact has a single-step error of order 4 t, i.e., smaller than what might have been expected from the initial expansion (5). This can be seen most easily by deriving the method in another way, starting from two dierent Taylor expansions: 1 an + xn+1 = xn + t vn + 2 2 t 1 an x n 1 = x n t vn + 2 2 t 3 1 3 a n + O(4 t ), 6 t 1 3 a n + O(4 t ). 6 t

(9)

## Adding these two equations gives

4 xn+1 = 2xn xn1 + 2 t an + O (t ),

(10)

which is called the Verlet algorithm. It does not contain any velocities explicitly; the next x-value is obtained from two preceding x values. However, the Verlet algorithm is in fact completely equivalent to the leapfrog method (7). To see this, we use the fact that (xn xn1 )/t equals the velocity vn1/2 , up to a correction O(2 t ), and write xn+1 = xn + t (vn1/2 + t an ) + O(4 t ). (11)

Noting that vn1/2 + t an = vn+1/2 , we get exactly the leapfrog method (7). It is important to note that the substitution xn xn1 t vn1/2 done in order to get Eq. (11) from Eq. (10) does not introduce any further errors in this equation, as it can be considered as a denition of the velocity in terms of a dierence between the two x-values. Hence, the step error in the position is 2 O(4 t ) in the leapfrog method. For the velocity the error is O (t ), due to the way it is dened as a discretized derivative. As we will see below in Sec. 2.3, the accumulated errors in x and v over a nite integration time both scale as (2 t ). An important feature of the leapfrog algorithm is that it is time-reversible, i.e., if we calculate xn+1 from xn and vn+1/2 according to Eqs. (7), and then reverse the time direction to go back to xn , we use the same velocity vn+1/2 as in the forward step, up to the change of its sign. Therefore, in the absence of numerical round-o errors, we arrive back exactly at the original xn . This is not the case with the Euler metod, Eqs. (3), where we use vn to go from xn to xn+1 , but vn+1 to go backwards from xn+1 to xn . If we carry out several leapfrog integration steps backwards in time, t = tN , . . . , t0 , after having completed a calculation from t0 to tN , we should arrive back at the original initial conditions at t = t0 . This way the time-reversibility can in fact be used as a check of the sensitivity of a calculation to round-o errors (which will slightly break the symmetry). An important consequence of the time-reversibility is that the errors are bounded for a system with periodic motion. This follows because integrating backward or forward in time for a full period T involves the acceleration at exactly the same x points (assuming now that the period is exactly a multiple of the time step). Hence, the error in some quantity, e.g., the energy, has to be the same at t = T ; (T ) = (T ). Integrating forward from t = T to 0, starting from the point reached in a previous integration from t = 0 to T , we would expect an additional error similar to (T ), because the initial conditions of these two integrations only dier by the very small amount (T ), so that the total error after a backward and a forward integration should be 2(T ). To be consistent with the time-reversibility, according to which this second integration in fact must bring us back to the original starting point and hence zero error, we must have (T ) = 0. The error only vanishes completely at times t t0 = T, 2T, . . ., but the important point is that there can be no steady error increase, which does aect methods with no time-reversal symmetry (as shown for the Euler method in Fig. 1). The general form of the long-time error scaling is discussed in detail in Sec. 2.3. The bounded error of the leapfrog method when applied to systems with periodic motion is illustrated in Fig. 3, which shows results for the energy of the same harmonic oscillator as the one studied in Sec. 2.1 using the Euler method. The energy here oscillates with a period half of the periodicity 2 of the system. The vanishing of the error at every half period in this case is a consequence of the symmetric potential, which was not assumed in the discussion above. Note also 4

Figure 3: Time dependence of the energy of a harmonic oscillator with k = m = 1, integrated using the leapfrog method with t = 101 (thinner curve) and 102 (thicker curve). how small the deviations are from the true energy E = 1/2 (compare with Fig. 1), even for t as large as 0.1. There is yet a third equivalent formulation of the Verlet method, called the velocity Verlet algorithm. It is obtained from the original Verlet algorithm as follows: Adding xn+1 on both sides of Eq. (10),
4 2xn+1 = xn+1 + 2xn xn1 + 2 t an + O (t ).

(12)

we can dene the velocity using a two-step dierence; vn = 1 (xn+1 xn1 ), 2t (13)

which when used in (12) gives the position in terms of x and v at the previous step only; 1 xn+1 = xn + t vn + 2 an . 2 t (14)

In order to obtain an equation for the velocity, we rst write the original Verlet equation (10) for xn instead of xn+1 ; xn = 2xn1 xn2 + 2 (15) t an 1 , which we add to (10). Rearranging the result in the following way; xn+1 xn1 = xn xn2 + 2 t (an1 + an ), and using the velocity denition (13), we obtain an equation for vn ; 1 vn = vn1 + t (an1 + an ). 2 (17) (16)

f(t)

t
Figure 4: A function that is evaluated on a grid with spacing t (solid circles) can be approximated by an interpolating polynomial between those points (dashed curve), provided that the function is smooth on the scale t . Shifting the time step by one and again writing down Eq. (14) for the position, we arrive at the velocity Verlet algorithm: 1 an , xn+1 = xn + t vn + 2 2 t 1 vn+1 = vn + t (an+1 + an ). 2

(18)

This formulation of the Verlet algorithm is completely equivalent to (10) and (7) as far as the propagation of the position is concerned. It may preferrable to the leapfrog method in cases where there is some reason to use the same time grid for x and v , but note that more operations are required at each step of the velocity Verlet algorithm (however, the number of calls to the force/acceleration function is the same, which is typically the part dominating the processor time). It should also be noted that although the algorithms give identical results (up to round-o errors) for xn , the accuracy in the velocity is four times higher in the leapfrog method becuase it is there dened using x points separated by a single time step, instead of two steps in the velocity Verlet method. However, nothing of course prohibits us from also calculating vn+1/2 in terms of xn and xn+1 in the velocity Verlet method or vn in terms of xn1 and xn+1 in the leapfrog method. In view of this, all results are completely equivalent in the two formulations.

2.3

## Error build-up in the Verlet/leapfrog method

As we have seen above, the single-step error in the Verlet/leapfrog algorithm is O(4 t ) for the position and O(2 ) for the velocity. We are clearly also very interested in the accumulated errors t after a large number of steps have been carried out. To study the accumulated error in the position x calculated using the Verlet/leapfrog method, we use the original Verlet form of the algorithm, Eq. (10). We introduce a symbol for the deviation n of the calculated value xn from the exact solution xex n , xn = xex n + n , 6 (19)

## and insert this in Eq. (10). After rearranging, this gives

ex ex 2 4 n1 2n + n+1 = (xex n1 2xn + xn+1 ) + t an + O (t ).

(20)

Under the assumption that the true solution, the acceleration (i.e., the force), and the deviation are all smooth functions on a time scale set by t (which is true for t smaller than some max if t the true solution and the force are well behaved functions), we can imagine constructing high-order polynomials that go through all the (tn , n ) and (tn , xex n ) points, as illustrated in Fig. 4. We can then use these interpolating polynomials to work in continuum time and use the second derivative of a function [f (t) = (t) or xex (t)] in place of the discretized versions appearing in Eq. (20), 1 1 d d2 f (tn ) [f (tn+1/2 ) f (tn1/2 )] = 2 (fn1 2fn + fn+1 ), = dt2 t dt t to obtain from (20) (tn ) = x ex (tn ) + a(tn ) + O(2 t ), t 0, (21)

(22)

even though t is xed and nite. The error introduced by this quasi-continuum approach is given by the order of the imagined interpolating polynomials and will clearly be much smaller than the O(4 t ) error of the Verlet formula we are working with (under the assumption of smoothness on the scale set by t ). By denition of the exact solution, x ex (tn ) = a(tn ), and hence in Eq. (20) we are left with (tn ) 2 . (23) t We can write an expression for the error (T ) after a large number of integration steps in terms of the second derivative (T = tN t0 ; we will set t0 = 0):
T

(T ) =
0

(t) + (0) = dt
0

dt
0

(t ) (0) + (0). dt

(24)

By denition of the initial conditions, the error (0) = 0. Since the Verlet formula is completely symmetric with respect to reversal of the time direction (n n), the Taylor expansion of the (0) = 0. error around t0 can have no contribution of rst order (or any odd order), and hence also With the second derivative bounded by the quadratic form (23), the integral (24) clearly can scale no worse than as T 2 2 t . Hence, the error in x after a nite number of steps scales in the worst case in the time discretization and as T 2 in the total integration time. Since there is an scenario as 2 t unknown prefactor in the t dependence of the error in Eq. (23)the sign can even be mixedthe T -scaling can be much better in practice, as we have seen explicitly in Fig. 3 in the case of the harmonic oscillator. Note again that the above arguments hold only when t is suciently small for the solution to be smooth on this scale. The way the velocity is dened as a discrete derivative of the coordinate, its error will clearly also scale as 2 t for any t.

2.4

## Verlet/leapfrog methods for damped systems

In the derivation of the Verlet algorithm, and its leapfrog formulation, we assumed that the forcefunction contains no explicit velocity dependence. With a velocity dependent force, F (x, v, t), the problem is that the acceleration an = a(xn , vn , tn ) in (7) should be evaluated at the time-points

corresponding to the position xn , on which we do not have the velocity. To circumvent this problem, we rst separate the damping term from the rest of the force terms and write 1 a(x, v, t) = [F (x, t) G(v )]. (25) m If we simply approximate a(xn , vn , tn ) by [F (xn , tn ) G(vn1/2 )]/m, we make an error of order t . Since a is multiplied by 2 t to obtain xn+1 , the resulting algorithm acquires a contribution to the x error scaling as 3 ; an order lower than the Verlet algorithm without damping term. To remedy t this, we can use a scheme based on an intermediate estimation x n+1 of the position, obtained using the above approximation G(vn1/2 ) in place of G(vn ). We can subsequently correct for the O(3 n+1 by doing a second step, where we use a better approximation of the velocity; t ) error in x vn = ( xn+1 xn1 )/(2t ). The error in this velocity, and hence in G(vn ), is O(2 t ), i.e., one order smaller than in the rst approximation and sucient to render xn+1 calculated with it accurate to the same order as in the original leapfrog algorithm. To summarice this modied algorithm, these are the steps to be performed in the leapfrog version of the Verlet algorithm with a damping term: v n+1/2 = vn1/2 + t [F (xn , tn ) G(vn1/2 )]/m, x n+1 = xn + t v n+1/2 , vn = ( xn+1 xn1 )/(2t ), vn+1/2 = vn1/2 + t [F (xn , tn ) G(vn )]/m, xn+1 = xn + t vn+1/2 . This algorithm requires more work than the simple leapfrog method without damping. However, in cases where the processor time is dominated by evaluating the function F , the dierences are in practice only minor. In the important special case of linear damping, i.e., xn+1 xn1 + O(2 G(v ) = v = t ), 2t the algorithm can be simplied. Starting from the Verlet form (10) we can write xn+1 = 2xn xn1 + which we can rearrange as xn+1 (1 + t /2m) = 2xn xn1 (1 t /2m) + or xn+1 = xn + (1 t /2m)(xn xn1 ) + 2 t Fn /m , 1 + t /2m t (1 t /2)vn1/2 + 2 t Fn = xn + . 1 + t /2 2 t Fn + O(4 t ), m (29) 2 t [Fn (xn+1 xn1 )/(2t )] + O(4 t ), m (28) (26)

(27)

(30)

In this expression we can identify the velocity vn+1/2 and obtain the leapfrog algorithm in the presence of linear damping; vn+1/2 = xn+1 (1 t /2m)vn1/2 + t Fn /m , 1 + t /2m = xn + t vn+1/2 . 8

(31)

Using the velocity Verlet algorithm (18) with damping, we have a problem analogous to that in the leapfrog method; to evaluate vn+1 we need an+1 , which itself depends explicitly on vn+1 . We can proceed in a way similar to what we did above, rst obtaining an estimate v n+1 based on approximating an+1 = a(xn+1 , vn+1 , tn+1 ) by a(xn+1 , vn + an t , tn+1 ) and then rening; 1 xn+1 = xn + t vn + 2 an , 2 t 1 v n+1 = vn + t [an + a(xn+1 , vn + t an , tn+1 )], 2 1 vn+1 = vn + t [an + a(xn+1 , v n+1 , tn+1 )]. 2

(32)

Also in this case the algorithm can be further simplied in the case of linear damping. The expression corresponding to Eq. (16) in the presence of linear damping is xn+1 xn1 = xn xn2 + 2 t (Fn1 + Fn ) t (xn+1 xn1 + xn xn2 ), m 2 (33)

which leads to the following algorithm [the formula for the position in the original velocity Verlet algorithm with no damping, Eq. (18), remains unchanged]: 1 xn+1 = xn + t vn + 2 an , 2 t t 1 [vn (1 t /2m) + (Fn1 + Fn )]. vn+1 = 1 + t /2m 2m

(34)

2.5

## The Runge-Kutta method

The Runge-Kutta (RK) method can be considered the most classic of all high-order schemes. There is a whole range of methods of this type of dierent orders, but the RK name is most commonly associated with the fourth-order variant [discretization error O(5 t )]. We will here outline one way of constructing this algorithm; however, without proving the error scaling. We will rst give a complete proof of the much simpler second-order RK method. For simplicity, before considering the equations of motion, we will consider the slightly easier case of a single rst-order dierential equation, x (t) = f [x(t), t]. (35)

The second-order RK method for this equation corresponds to the mid-point rule of function integration: If we know the value of the function f [x(t), t] at the mid-point of an interval [tn , tn+1 ], its integral is known up to an error O(3 t ):
tn+1 tn

## f [x(t), t]dt = t f [x(tn+1/2 ), tn+1/2 ] + O(3 t ).

(36)

In the case of a function f (t), we could simply evaluate it at tn+1/2 , but in the present case we do not know the value of xn+1/2 = x(tn+1/2 ) and hence we need to approximate it. We only need an

approximant linear in t to keep the error of the integral cubic, so we can simply use the Euler formula for (35): t f (xn , tn ) + O(2 (37) x n+1/2 = xn + t ), 2 where is again used to indicate an intermediate step of the calculation. We then arrive at the second-order RK formula, xn+1 = xn + t f ( xn+1/2 , tn+1/2 ) + O(3 t ), which is typically written in the algorithm form k1 = t f (xn , tn ), k2 = t f (xn + k1 /2, tn+1/2 ), xn+1 = xn + k2 . In the case of the equations of motion for x and v , we use the Euler formula for both quantities in the rst step to obtain approximants for xn+1/2 and vn+1/2 and then use these to obtain better estimates in the same way as we did above, in Sec. 2.4, giving k1 = t a(xn , vn , tn ), l1 = t vn , k2 = t a(xn + l1 /2, vn + k1 /2, tn+1/2 ), l2 = t (vn + k1 /2), vn+1 = vn + k2 , xn+1 = xn + l2 . This algorithm is rarely used in practice, since the Verlet/leapfrog method is both simpler and more accurate. The algorith does, however, serve as a good warm-up for studying the signicantly more complicated fourth-order RK algorithm. We again rst consider the single equation (35). Again using results from simple function integration, we know that if the values of the function f [x(t), t] are evaluated at the time points tn , tn+1/2 , and tn , we can apply Simpsons formula to obtain the integral over the range [tn , tn+1 ] up to an error O(5 t ), i.e., t (fn + 4fn+1/2 + fn+1 ) + O(5 (41) xn+1 = xn + t ). 6 We hence need to nd approximations for fn+1/2 and fn+1 up to errors O(4 t ). The somewhat obscure scheme used for this purpose is illustrated in Fig. 5. As in the second-order derivation above, we begin by using an Euler formula to obtain an approximation for xn+1/2 . This time there will be two such estimates, and we use a superscript to distinguish between them; x 1 n+1/2 = xn + x 2 n+1/2 t f (xn , tn ) 2 t f ( x1 = xn + n+1/2 , tn+1/2 ). 2 (40) (39) (38)

(42)

Since f is the derivative of x, what we have done is to rst extrapolate xn+1/2 from xn using derivative information at the initial point (xn , tn ), denoted z0 in Fig. 5. The second extrapolation 10

Figure 5: Illustration of the scheme used to estimate the values of the function f at the true (t, x) points (tn+1/2 , xn+1/2 ) and (tn+1 , xn+1 ) (here shown as open circles) based on the known point z0 = (tn , xn ) (large solid circle). A point z1 is rst generated using f (z0 ) as the derivative. The points z2 and z3 are obtained using f (z1 ) as the derivative. is based on derivative information at the estimated point ( x1 n+1/2 , tn+2 ), denoted z1 in the gure, but still using (xn , tn ) as the point of departure. This second point at tn+1/2 is denoted z2 in the gure. The idea of this procedure is, roughly, that if the function f is smooth on the scale t , the actual point xn+1/2 should be between z1 and z2 . In the RK method the desired function value f1+n/2 is approximated by the average of the function at these two estimated points; 1 2 fn+1/2 [f (x1 n+1/2 , tn+1/2 ) + f (xn+1/2 , tn+1/2 )]. 2 (43)

We will not prove here the remarkable fact that this approximation has the required accuracy, i.e., x2 an O(4 t ) error. In order to estimate xn+1 , the function value at ( n+1/2 , tn+1/2 ) (i.e., z2 ) is used to obtain an approximation of the derivative, x n+1 = xn + t f (x2 n+1/2 , tn+1/2 ). (44)

This is in the spirit that the derivative should be evaluated at the mid-point between (xn , tn ) and (xn+1 , tn+1 ). Not knowing the actual point, we use the estimate z2 for this. Again, it can be proven that this leads to an error O(4 t ) in xn+1 . The fourt-order RK algorithm resulting from the above procedures is traditionally written as k1 = t f (xn , tn ), k2 = t f (xn + k1 /2, tn+1/2 ), k3 = t f (xn + k2 /2, tn+1/2 ), k4 = t f (xn + k3 , tn+1 ), 1 xn+1 = xn + (k1 + 2k2 + 2k3 + k4 ). 6 (45)

The correctness of this scheme up to an error of order O(5 t ) can be proven using a Taylor expansion. Before adapting the RK algorithm to equations of motion, we write down a more general form of

11

the RK scheme for two coupled equations; x (t) = f (x, y, t), y (t) = g(x, y, t), for which the RK algorithm generalizes to k1 = t f (xn , yn , tn ), l1 = t g(xn , yn , tn ), k2 = t f (xn + k1 /2, yn + l1 /2, tn+1/2 ), j2 = t g(xn + k1 /2, yn + l1 /2, tn+1/2 ), k3 = t f (xn + k2 /2, yn + l2 /2, tn+1/2 ), l3 = t g(xn + k2 /2, yn + l2 /2, tn+1/2 ), k4 = t f (xn + k3 , yn + l3 , tn+1 ), l4 = t g(xn + k3 , yn + l3 , tn+1 ), 1 xn+1 = xn + (k1 + 2k2 + 2k3 + k4 ), 6 1 yn+1 = yn + (l1 + 2l2 + 2l3 + l4 ), 6 which in turn easily generalizes to any number of coupled equations. In the case of the equations of motion, we make the identication x v , y x, f a, g v (note that v corresponds to a function g(x, y, t) = x in the notation used above). We then obtain k1 = t a(xn , vn , tn ), l1 = t vn , k2 = t a(xn + l1 /2, vn + k1 /2, tn+1/2 ), l2 = t (vn + k1 /2), k3 = t a(xn + l2 /2, vn + k2 /2, tn+1/2 ), l3 = t (vn + k2 /2), k4 = t a(xn + l3 , vn + k3 , tn+1 ), l4 = t (vn + k3 ), 1 vn+1 = yn + (k1 + 2k2 + 2k3 + k4 ), 6 1 xn+1 = xn + (l1 + 2l2 + 2l3 + l4 ). 6 (46)

(47)

(48)

The RK method is clearly not time-reversal symmetric, and hence the errors are unbounded even for periodic motion. Fig. 6 shows the time dependence of the energy obtained using the RK method for the same harmonic oscillator that was previously studied with the Euler (Fig. 1) and leapfrog (Fig. 3) methods. Comparing the RK and leapfrog methods at = 0.1, it can be seen that the errors are considerably smaller in the case of the RK method within the time range shown. However, the error grows linearly with time, and hence for long times the leapfrog method will in fact perform better in this case although its single-step error scaling is worse. However, for aperiodic motion, or motion with a very long period, there may be no practical advantage of the bounded error of the 12

t=0.2 t=0.1

E(t)

## 0.49994 0.49992 0.49990 0.49988 0

10

20

30

40

50

Figure 6: Time dependence of the energy of a harmonic oscillator with k = m = 1 integrated using the 4th-order Runge-Kutta method with time steps t = 0.2 and 0.1. leapfrog methodthe error can still grow large within a periodand then the RK method is often preferrable. An advantage of the RK method is that it can be used with a variable time step. This is not possible in the Verlet/leapfrog method, because data at tn and tn+1/2 are needed to advance the time to tn+1 . Using dierent time-steps tn+1 tn and tn tn1 would clearly ruin the symmetry of the scheme and lead to a worse error scaling. In the RK method, only data at tn are needed to advance to tn+1 , and hence t can be chosen arbitrarily in each step. There are adaptive methods that use this property of the Runge-Kutta scheme to automatically adjust the time step to keep the error within specied limits.

2.6

## Motion in more than one dimension

The methids we have discussed above generalize very easily to motion in more than one dimension. The equations of motion become vector equations; x (t) = v (t) (t) = 1 F [x(t), v (t), t], v m

(49)

which separate into equations for the components x , v ( = 1, 2, 3 in three dimensions). These equations are coupled only through the force function, the components F of which typically depend on all components of x and, in the case of damped motion, v . As a simple example, we consider planetary motion, which takes place in a single plane (in the case of a single planet considered here), i.e., we study the equations of motion in two dimensions. For a planet of mass m moving in the gravitational eld of a star with mass M , the force experienced by m is given by GM m r, (50) F (r ) = r3 13

where r is the position vector of the planet with respect to the star and r = |r |. In principle this is a two-body problem, but a result of elementary mechanics is that it can be reduced to a one-body problem for a reduced mass = M m/(m + M ). We will here for simplicity assume that M m and treat the star as a stationary body, using the original m for the mass of the planet. We then have the equations of motion for the x and y coordinates of the planet: x = vx , v x = GM x/r 3 , y = vy , v y = GM y/r , where r = x2 + y 2 . The leapfrog algorithm (7) for these equations is vx (n + 1/2) = vx (n 1/2) t GM x(n)[x2 (n) + y 2 (n)]3/2 , x(n + 1) = x(n) + t vx (n + 1/2), vy (n + 1/2) = vy (n 1/2) t GM y (n)[x2 (n) + y 2 (n)]3/2 , y (n + 1) = y (n) + t vy (n + 1/2). This scheme of course generalizes to other types of forces. The Runge-Kutta method can also easily be used in more than one dimension. (52)
3

(51)

2.7

## Program example; taking down a satellite

As an example of a numerical solution of equations of motion, we will here show how to calculate the orbit of a satellite that is subject to earths gravitational force, atmospheric drag, as well as the thrust of its own internal rocket motor. Specically, we will study the case of a satellite that is to be brought back to earth after having served its mission in a circular orbit, at a specied altitude. At the start of the calculation, the rocket motor, which produces a force directed opposite to the satellites velocity vector, is turned on for a specied duration of time. As the satellite slows down, it exits the orbit and begins to accelerate towards earth. As it enters the atmosphere, air drag will slow down the fall, and eventually the satellite crashes onto the ground. The program, which is available on the course web-site in the le crash.f90, produces a le with the time dependence of the altitude, the angular position, and the velocity. 2.7.1 Modeling the forces acting on a satellite

The problem is illustrated schematcally in Fig. 7. The satellite has mass m and its position and velocity vectors are r and v , respectively. Its dynamics is governed by a force F = Fgravity (r )er + Frocket (t)ev + Fdrag (r, v )ev , (53)

which is a sum of forces due to gravitation, the internal rocket motor, and air drag. All the forces in this problem keep the motion in the plane dened by the original circular orbit, and hence we consider a two dimensional space; r = xex + yey . The explicit dependencies of the forces on r , 14

Figure 7: Schematic illustration of the satellite landing problem. Earth, seen from above the north pole, is shown in dark gray; the atmosphere is indicated by the lighter shade of gray. The satellite is initially in a circular orbit around the equator. At the start of the calculation it is at the position indicated, at angle = 0 with the x axis. The rocket motor is turned on and the satellite begins to descend, eventually crashing at a location to be determined by the calculation. v , and t have been indicated above; the unit vectors are in the direction of r and v , respectively; er = r/r and ev = v/v . We will now specify their exact forms of the force terms. We already discussed motion in earths gravitational eld above in Sec. 2.6; the force Fgravity is given by Eq. (50), where M is earths mass, Me = 5.976 1024 kg, and the constant of gravitation G = 6.672 1011 m3 /kgs2 . To relate the initial altitude h0 to r , and to determine when the satellite crashes, we also need earths radius; re = 6378 km. For the braking process used to bring the satellite out of its orbit, we assume that the rocket engine produces a thrust corresponding to a constant decelerating force opposite to the direction of motion. We neglect the reduction of the satellites mass as fuel is consumed, i.e., we assume a contribution to the total acceleration that is constant in magnitude; Frocket = (Tbreak t)Bev , m (54)

where Tbrake is the time the motor stays on, and (t) is the theta function; (t) = 1 for t 0, (t) = 0 for t < 0. In the program, we will take B = 5m/s2 for the deceleration. We also need a model of atmospheric drag. We will assume that the damping force due to the air depends quadratically on the velocity; Fdrag (h) = Cd (h)v 2 ev (55)

where Cd is the drag coecient and (h) is the density of the atmosphere at altitude h and h = r re . The quadratic form is appropriate at high velocity, whereas a linear velocity dependence is closer to reality at low velocities. Here we will not incorporate a cross-over to such a low-velocity linear form, because the velocity will be high most of the time and the details would anyway be highly dependent on the shape of the satellite. We will adjust the drag coecient Cd , which has the unit m2 , so that a reasonable terminal velocity for a free fall is achieved. 15

(kg/m )

0.5

0 0

10

20

30

40

r (km)

Figure 8: Density of the atmosphere as a function of altitude above sea level. The circles are values for the so-called 1976 standard atmosphere, and the curve was calculated using the model (57) with k1 = 1.2 104 m and k2 = 2.2 104 m. The terminal velocity is determined by Fdrag (h = 0, v )/m = g, or g Cd = . m (0)v 2 (56)

Our values for GM and earths radius re correspond to a gravitational acceleration g = 9.77355 m/s2 . Since we are just interested in a rough estimate, we will neglect the fact that the true terminal velocity will be aected by the non-constant air denity, and, to a lesser extent, the altitude dependence of g. We can later use our program to check the true terminal velocity of a body of given Cd /m that drops from a high altitude. Let us here assume that the shape and mass of the satellite are such that the terminal velocity (56) is 100m/s if and g are taken as constants equal to their sea-level values. With the density of air at sea level taken as 1.225 kg/m3 , we then get Cd /m 8 104 m2 /kg, which we will use in the program. We will use the following simple model for the atmospheric density (h) = exp h + k1 h k3/2
3/2

1.225 kg/m3 ,

(57)

with k1 and k2 adjusted to give an approximate agreement with tabulated values for the density. This form was arrived at simply by trial and error, using tabulated data up to r = 40km. With k1 = 1.2 104 m and k2 = 2.2 104 m the agreement with the data is good, as shown in Fig. 8. For altitudes > 40 km it is not clear how well this model reproduces actual conditions. In fact, the atmospheric density at very high altitudes varies considerably with time and is inuenced by, e.g., the solar wind (and, of course, even at low altitudes the actual air density depends on the weather!). The model density (57) suces for semi-realistic illustration purposes. 2.7.2 Computer program

The computer program carrying out the calculations is available on the course web site, in the le crash.f90. We here discuss its main subroutines and functions. 16

The program uses a module systemparam, which appears at the beginning of the program le and contains some of the parameters and variables used in several of the program units:
module systemparam real(8), real(8), real(8), real(8), real(8), parameter parameter parameter parameter parameter :: :: :: :: :: pi=3.141592653589793d0 gm=3.987d14 ! Earths mass times constant of gravity arocket=5.d0 ! Retardation due to rocket motor dragc=8.d-4 ! Air drag coefficient re=6.378d6 ! Earths radius

## end module systemparam

The parameters, all of which were discussed above, are explained in the module by comments; SI units are used throughout the program. The time step dt and the run-time tbrake of the rocket motor are read from the keyboard by the main program. A variable dt2 with half the time step is also dened. The statement use systemparam appears in the beginning of most of the program units, which thus have access to all the declarations of the modeule. The main program satellitecrash reads the following input data from the user; the initial altitude, the run-time of the rocket motor, the time step to be used in the integration of the equations of motion, the fraction of the time-points to be recorded, and the maximim time to integrate (used in case inappropriate input data are used, for which the satellite may never crash, or may do so only after an exceedingly long time). Initial conditions for the position and the velocity are set according to the altitude entered (to which the earths radius is added to give the distance r0 from earths center). The initial position is x = r0 , y = 0, in accordance with Fig. 7. The initial velocity is that corresponding to a circular orbit of radius r0 ; vx = 0, vy = GM/r0 . The main program repeatedly calls a subroutine rkstep, which uses the Runge-Kutta algorithm to advance the time by a step t =dt. Before each call to this routine, the following information is written to a le sat.dat: the time, the normalized angular position (from 0 to 1, with 0 corresponding to the initial angle 0, as calculated using the x and y coordinates in the subroutine polarposition), the altitude, and the velocity. The program stops when a negative altitude h = r re has been reached (the satellite has crashed), or when the maximum integration time is exceeded. In the integration subroutine rkstep(t0,x0,y0,vx0,vy0), the argument list contains the present time, coordinates, and velocities. On exit, the coordinates and velocities have been replaced by their values at the next time step. When the subroutine is executed, the time value at the next step t1 as well as the half-step th between t0 and t1 are rst set; t1=t0+dt; th=t0+dt2. The implementation of the Runge-Kutta algorithm (48) for the x and y components follows:
call accel(x0,y0,vx0,vy0,t0,ax,ay) kx1=dt2*ax ky1=dt2*ay

17

lx1=dt2*vx0 ly1=dt2*vy0 call accel(x0+lx1,y0+ly1,vx0+kx1,vy0+ky1,th,ax,ay) kx2=dt2*ax; ky2=dt2*ay lx2=dt2*(vx0+kx1) ly2=dt2*(vy0+ky1) call accel(x0+lx2,y0+ly2,vx0+kx2,vy0+ky2,th,ax,ay) kx3=dt*ax ky3=dt*ay lx3=dt*(vx0+kx2) ly3=dt*(vy0+ky2) call accel(x0+lx3,y0+ly3,vx0+kx3,vy0+ky3,t1,ax,ay) kx4=dt2*ax ky4=dt2*ay lx4=dt2*(vx0+kx3) ly4=dt2*(vy0+ky3) x1=x0+(lx1+2.d0*lx2+lx3+lx4)/3.d0 y1=y0+(ly1+2.d0*ly2+ly3+ly4)/3.d0 vx1=vx0+(kx1+2.d0*kx2+kx3+kx4)/3.d0 vy1=vy0+(ky1+2.d0*ky2+ky3+ky4)/3.d0

In these lines of code, note that the variables corresponding to k1 , k2 , k4 , l1 , l2 , l4 (but not k3 , l3 ) in Eq. (48) have been dened with a factor 1/2 (using the half time step dt2) in order to avoid unnecessary divisions by 2. At the end of the subroutine the new coordinates and velocities are copied to the variables of the argument list; x0=x1; y0=y1; vx0=vx1; vy0=vy1. The accelerations needed in the Runge-Kutta algorithm are obtained by calls to the subroutine accel(x,y,vx,vy,t,ax,ay), where ax and ay are the accelerations in the x and y directions. This subroutine implements the various F/m terms discussed above, in the following way:
r=sqrt(x**2+y**2) v2=vx**2+vy**2 v1=sqrt(v2) !*** evaluates the acceleration due to gravitation r3=1.d0/r**3 ax=-gm*x*r3 ay=-gm*y*r3

!*** evaluates the acceleration due to air drag if (v1 > 1.d-12) then ad=dragc*airdens(r)*v2 ax=ax-ad*vx/v1 ay=ay-ad*vy/v1 end if !*** evaluates the acceleration due to rocket motor thrust if (t < tbrake .and. v1 > 1.d-12) then ax=ax-arocket*vx/v1 ay=ay-arocket*vy/v1 end if

18

1500

1500

v (m/s)

1000

1000

500

500

50

## 100 150 200 250

50

100

150

200

t (s)

h (km)

Figure 9: Time and altitude dependence of the velocity a body dropped from an altitude of 200 km in the presence of the atmosphere (solid curves) and without atmosphere (dashed curves). The operations here should be self-explanatory; the air drag calculation uses a function airdens(r) which implements the athmospheric density model (57). The condition v1 > 1.d-12 used in the air drag and rocket motor contributions is included for safety, in order to avoid an unlikely division by zero. In the case of the air drag, the force at zero velocity is zero so there should be no contribution anyway. In the case of the rocket motor contribution, a zero velocity would happen only if the motor runs long enough for the satellite to come to a halt, which should not happen in situations of interest (the way the brake force is dened, it always acts in the direction opposite to the velocity, and hence one cannot make the satellite y o into space with this program). If the velocity happens to be exactly zero at some time step, the motor will be neglected for that step. 2.7.3 Some results

In order to test the program, we can give a high initial altitude, > 200 km, for which there is no signicant air drag, and run the motor for 0 seconds. The satellite should then continue forever in its circular orbit. We can also set the atmospheric drag coecient dragc=0.d0 and the initial velocity vy=0 (instead of the value of the circular orbit used). The satellite should then experience free fall. If we use a very low altutude, e.g., 1 km, the gravitational force is constant to a good approximation, and the time to reach the ground should be t = 2r0 /g (with g 9.774 m/s2 ). The program passed these tests. Nest, we will calculate the actual terminal velocity of a body falling fram a high altitude in the presence of the atmosphere. We have chosen the drag coecient Cd /m = 8 104 m2 /kg, so that the terminal velocity according to the relation (56) is 100 m/s (our values correspond more precisely to 99.865 m/s). The actual terminal velocity will be dierent, however, because of the altitude dependence of g and . Since the atmospheric eect is larger, we would expect a slightly higher terminal velocity. Running the program with an initial velocity 0 and an initial altitude of 200 km produced the time and altitude dependencies of the velocity shown in Fig. 9, where results are also

19

100 80

8000

6000

h (km)

60 40 20 0

v (m/s)
0 0.2 0.4 0.8 1

4000

2000

/2

0.6

0.2

0.4

/2

0.6

0.8

Figure 10: Altitude and velocity versus angular distance traveled for a satellite initially in a circular orbit at 100 km and then subjected to atmospheric drag. shown for the case of no atmosphere (obtained by setting dragc=0.d0). In these graphs the eects of the atmosphere are seen to become strong quite sharply at an altitude of approximately 40 km. The terminal velocity is 102.6 m/s. These calculations were run with a time step of 0.1 s, however, 1 s produces almost indistinguishable results (but gives a less precise crash time). It is also interesting to investigate the slow-down solely due to atmospheric drag of a satellite in a low orbit, i.e., without igniting its rocket motor. With the atmospheric density model used, the density at 200 km altitude is only 1019 of that at sea-level, and the satellite stays in orbit almost indenitely (in reality a satellite at that altitude must be regularly propelled to stay in orbit). At 100 km, the relative density in our model is 1.5 108 , which actually has a signicant eect on the orbit; by coincidence, the satellite in this case crashes almost exactly after one revolution around earth, as shown in Fig. 10. Notice that the magnitude of the velocity is aected very little until the satellite has dropped to 50 km. It is hence primarily the direction of the velocity that is aected at the initial stages of the decent. The calculation was run with t = 1 s, but even 10 s gives an acceptable convergence in this case. Fig. 11 show results for an initial altitude of 120 km. In this case the satellite stays up for tens of orbits, eventually crashing after almost 88 hours. The points represent the altitude versus the angular position sampled every 30 s (the integration was done using t = 1 s). Finally, we turn on the rocket engine. We investigate a satellite in a 200 km orbit, keeping the rocket engine running for various lengths of time. Results for the altitude versus the angular position are shown in Fig. 12. Running the motor for 5 s does not suciently slow down the satellite and it instead assumes an almost elliptical orbit; it will eventually crash, but only after a very large number of orbits. After 10 s and 20 s motor runs the satellite crashes without completing another full orbit.

20

Figure 11: Altitude versus angular position of a satellite initially in a circular orbit at 120 km and then, at the beginning of the simulation, subjected to atmospheric drag. The points represent sampling of the path at intervals of 30 s.

200

150

h (km)

100
Tbreak = 5 s Tbreak = 10 s Tbreak = 20 s

50

0 0

0.2

0.4

/2

0.6

0.8

Figure 12: Altitude versus the angular position of a satellite initially in a circular orbit at 200 km altitude. The rocket engine, which causes the satellite to decelerate at 5 m/s2 , is turned on at angular position = 0 for dierent durations of time. In the case of the 5 s burn, the satellite stays in orbite for a long time in an almost elliptical orbit; only the rst orbit is shown here.

21

Figure 13: A pendulum of mass m, suspended by a weightless rod of length l. The denition of the angle x is indicated.

Chaotic motion

Many mechanical systems exhibit chaotic motion in some regions of their parameter spaces. Essentially, the term chaotic motion, or chaos, refers to aperiodic motion and sensitivity of the time evolution to the initial conditions. A chaotic system is in practice unpredictable on long time scales, althoght the motion is in principle deterministic, because minute changes in the initial conditions can lead to large changes in the behavior after some time (the time-scale on which this unpredictability sets in is governed by the so-called Lyapunov exponent). Although a chaotic system is unpredictable, its motion is not completely random. In particular, the way the system approaches chaos often exhibits universality, i.e., seemingly dierent systems make the transition from regular, periodic motion to chaotic motion in very similar ways, often through a series of quantitatively universal period dublings (bifurcations). We will here only touch briey on this subjecy, by studying one of the simplest model systems exhibiting chaotic motion; a periodically driven damped pendulum. We will limit the discussion to how chaotic motion can be characterized and quantied graphically based on numerical simulation data; methods that can be directly applied to other systems as well.

3.1

## The driven damped pendulum

For a pendulum of mass m suspended by a rigid rod of length l (assumed to be massless), the potential energy is V (x) = mgl[1 cos (x)], (58) where x is the displacement angle, as shown in Fig. 13. For small angles, this system can be approximated by a harmonic oscillator with spring constant k = mgl, but here we will keep the full form of the potential. In addition, we include a damping x , which could arise from the suspension mechanism at the upper end of the rod. We will also consider a periodic driving force Fdr = Q sin (t/), which could be due to, e.g., a time-dependent electrical eld if the pendulum is charged. The equation of motion is thus x = k sin (x) v + Q sin (t/), (59)

where v = x is the angular velocity. This equation cannot be solved analytically when = 0 and Q = 0. We will discuss the trajectory this system follows in its phase space x(t), v (t), using numerical results obtained with the Runge-Kutta method. 22

1 0.5 0.5

x(t)

v(t)

0.0

-0.5

-0.5

20

40

60

80

100

-0.5

0.0

0.5

1.0

x(t)

Figure 14: Damped motion of a pendulum with k = 1, = 0.1, Q = 0, and initial condition x = 1, v = 0. Left: displacement angle versus time; Right: phase space trajectory.

1 0.5

1.0 0.5

x(t)

0 -0.5 -1 0 20 40

v(t)

60

80

100

## -1.5 -1.0 -0.5 0.0

0.5

1.0

x(t)

Figure 15: Motion of a damped driven pendulum with k = 1, = 0.1, Q = 1, = 2, and initial condition x = 1, v = 0. Left: displacement angle versus time; Right: phase space trajectory.

23

Figure 16: Limit cycles of period 1, 2, and 4 times the driving period 2/ for a pendulum with k = 1, = 0.1, = 2/3, and three dierent values of the driving amplitude Q. Note that the angle x is periodic.

Figure 17: Limit cycles for a pendulum with k = 1, = 0.1, = 2/3, Q = 1.08, for two dierent initial conditions; x0 = 0.7 (left) and x0 = 0.8 (right); v0 = 0 in both cases. For = 0 and Q = 0, the pendulum exhibits simple periodic motion. In the presence of damping, > 0, it eventually comes to rest; this is illustrated by a plot of x(t) in the left panel of Fig. 14. Looking at the phase space, shown in the right panel of the same gure, the trajectory is seen to spiral toward the point (x = 0, v = 0). This point is said to be an attractor of this dynamical system. In Fig. 15 a driving force of amplitude Q = 1 and = 2 (period 2/ = ) has been added. In this case, the system is seen to initially exhibit damped, irregular motion (transient behavior), but eventually it settles into a periodic motion with the same frequency as the driving force. There is a phase-space attractor in this case as well; a closed loop. In general, an attractor is a region of phase space that a system approaches at long times for a set of boundary conditions. There can be more than one attractor, in which case the space of initial conditions subdivides into basins of attraction. The basins of attraction can be very complex (fractal), reecting a sensitive dependence on the initial conditions. An attractor is also often referred to as a limit set. A periodic attractor is a limit-cycle in this terminology. Periodic attractors more complicated than the single closed loop are also possible. Turning up the driving force Q, a damped pendulum ca undergo a series of period doublings, i.e., the period of the motion changes from 2/ to /, /2, etc., at certain values of Q. The corresponding limit-cycle 24

Figure 18: Bifurcation diagram for the k = 1, = 0.1, = 2/3, Q = 2/3 pendulum. in phase space then develops a more complex loop structure. Examples of such doublings are shown in Fig. 16. Here the graphing of the trajectory was started only after a time long enough for the transient behavior to have died away, so that only the attractor itself is seen. It can also be noted in these graphs that the attractor is not symmetric with respect to changing the signs of x and v . This is an example of symmetry breaking. There are also corresponding attractors which are obtained by (x, v ) (x, v ) in the gure (under which the equations of motion are symmetric). There are in fact yet other attractors in addition to those shown in Fig. 16 in this region of the parameter space. The attractors shown in Fig. 16 were obtained with the initial condition x0 = 1, v0 = 0. Fig. 17 shows two attractors for a pendulum with Q = 1.08 and the other paremeters equal to those in Fig. 16. In both ceses v0 = 0, but the initial angle x0 was dierent. For x0 = 0.7, shown in the left panel, the limit cycle is a symmetry-broken one with period four, as the one for Q = 1.081 in Fig. 16, but for x0 .8, shown in the right panel, a symmetric attractor is obtained. There is also of course an x x, v v counterpart to the symmetry-broken attractor. Which one of the three attractors that is realized (i.e., the basin of attraction to which the initial conditions belong) depends in a complex way on (x0 , v0 ). In order to systematically investigate the behavior versus a parameter of the system, e.g., the driving force Q, it is useful to not consider the full phase space, but only a cut through it (or, more generally, for a system with a higher-dimensional phase space, a surface in that space). For example, we can graph the velocity at the point where the angle x passes the value 0 from above. Such a construct is called a Poincar e section. For periodic motion, the Poincar e section is a discrete set of points; for the limit cycles shown in Fig. 16 there are 1, 2, and 4 points in the Poincar e section through x = 0. By graphing the set of points versus Q, we can investigate how the periodicity of the motion evolves (although in general it may be dicult, or even be impossible, to know exactly what Poincar e section to use in order to capture the actual periodicity this way, since a limit cycle may have a looping portion that does not cross the chosen cut). Such a plot is shown

25

Figure 19: The parameter region Q = 1.83 1.87 of the bifurcation diagram in Fig. 18 shown in greater detail. in Fig. 18, for the velocity as x crosses zero from above (hence all the velocities are negative). A very complex behavior can be observed, with regions of period doublings and regions where there does not appear to be a discrete set of points at allsuch regions correspond to chaotic behavior.1 Plots like these, often referred to as bifurcation diagrams, have been investigated in great detail for many systems. It turns out that the bifurcationsthe points where the periodicity doubleshave a universal structure, i.e., many dierent systems show a quantitatively similar behavior. Chaotic behavior occurs after an innite number of increasingly closely spaced period-doubling bifurcations, and this series of bifurcations exhibits universality. in Fig 19, the region Q = 1.83 1.87 of Fig. 18 is shown in greater detail (in order to get a sucient density of points, a new run with a ner Q-spacing had to be carried out). Here several period-doubling bifurcations can be seen preceding the chaotic regime for Q > 1.863. The universality of period-doubling bifurcations is contained in the so-called logistic map. This is a discrete recursion relation, where the next value in a sequence of numbers is given by xn+1 = xn (1 xn ). (60)

One can liken the discrete index n to the successive points on a Poincar e section. A graph similar to Fig. 18 is obtained by iterating the logistic map numerous times for each of a sequence of values of the control parameter , plotting the points (, xn ), n = 0, 1, 2, . . . (after some initial number of steps, to avoid transients). Such a graph is shown in Fig. 20. Clearly, there are great similarities with the bifurcations seen in Fig. 19, and in fact the two plots could be scaled in such a way that they coincide perfectly. This is true in general for systems that approach chaos through a
1 The same initial conditions, x0 = 1, v0 = 0 were used for all Q values in this plot. In regions where there are several basins of attraction, the plot could hence look dierent for other initial conditions.

26

Figure 20: Bifurcation diagram for the logistic map. In the right panel, a region of the left panel is shown on a more detailed scale, to reveal other bifurcation sequences. sequence of period-doubling bifurcations, as was rst noted by Mitchell Feigenbaum in 1974. The period-doublings are self-similar, i.e., they look exactly the same on dierent scales; the smaller bubbles: seen in Fig. 20 can be magnied and then look exactly as the bigger ones, and they themselves bifurcate into smaller bubbles, etc., until chaos sets in in the limit of an innite number of bifurcations. There is not just one period-doubling sequence in the logistic map; they appear in an innite number of places in a fractal fashion. An example is shown in the right panel of Fig. 20, which is a magnication of a region of the left panel where only some dim structures can be seen. Upon magnication, these features can be seen to constitute exactly the same type of bifurcation as in the left graph (but upside down). In addition to its famous period doubling bifurcations, the logistic map also exhibits other universal features, e.g., associated with the window structure (the almost empty spaces between the chaotic regimes, which correspond to periodic motion). As we have seen, period-doubling bifurcations do occur in the damped pendulum, and the similarity with the logistic map can be seen clearly, e.g., in Fig. 19. However, not all transitions to chaos in the pendulum are of this form; several routes tho chaos are known, but the period-doubling is the most well studied and understood. The phase space attractor for a pendulum in a chaotic state is shown in Fig. 21. In this case there is no limit-cycle, and instead the attractor lls a nite region of the two-dimensional phase space. There is clearly some structure in the attractor, with some regions more likely to be occupied than others. One can also represent the phase space behavior of a periodically driven systems by stroboscopic sampling. This amounts to plotting x and v at time intervals equal to the period 2/ of the driving force. For a system with periodicity equal to that of the driving force, this plot will clearly have a single point, and if the plot has n points the periodicity is n2/. For a chaotic systems the points fall on an extended region of the phase space, and this region has a fractal dimensionality for chaotic motion. Examples of stroboscopic sampling are shown in Fig. 22. Plots like these are also oftent referred to as Poincar e sections, and the sets of occupied points are also attractors. The right panel in Fig. 22 is an example of a strange attractor; one with a fractal dimensionality (i.e., a Hausdor dimensionality between 1 and 2 in this case). Stroboscopic sampling is often the easiest 27

Figure 21: Phase space trajectory (attractor) for the k = 1, = 0.1, = 2/3, Q = 1.2 pendulum; this system is chaotic.

Figure 22: Stroboscopic sampling plots for the k = 1, = 0.1, = 2/3 pendulum at two dierent driving apliitudes. The left graph shows a case where the period is 3, and the right graph shows a chaotic case with a strange attractor. way to examine whether a periodically driven system is chaotic.

28
{"ft_lang_label":"__label__en","ft_lang_prob":0.9019523,"math_prob":0.9972213,"size":57370,"snap":"2020-10-2020-16","text_gpt3_token_len":15551,"char_repetition_ratio":0.14536485,"word_repetition_ratio":0.05750389,"special_character_ratio":0.28175005,"punctuation_ratio":0.11786491,"nsfw_num_words":1,"has_unicode_error":false,"math_prob_llama3":0.998054,"pos_list":[0],"im_url_duplicate_count":[null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2020-04-05T08:02:26Z\",\"WARC-Record-ID\":\"<urn:uuid:4089561b-4f18-4e65-a7b1-b144b0246913>\",\"Content-Length\":\"421661\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:3ddf6d12-91ee-4a33-9863-aa2196bb1de2>\",\"WARC-Concurrent-To\":\"<urn:uuid:5f1b1c5e-e4d0-4fa3-ae90-99990b728624>\",\"WARC-IP-Address\":\"151.101.250.152\",\"WARC-Target-URI\":\"https://pt.scribd.com/document/184122466/c-Motion\",\"WARC-Payload-Digest\":\"sha1:3JFRSP74PLQUECQA44ORSKGQYDQMSY6X\",\"WARC-Block-Digest\":\"sha1:C4FJPECTLWYCHO7QIVCARK4MEQSNSX2D\",\"WARC-Identified-Payload-Type\":\"application/xhtml+xml\",\"warc_filename\":\"/cc_download/warc_2020/CC-MAIN-2020-16/CC-MAIN-2020-16_segments_1585370529375.49_warc_CC-MAIN-20200405053120-20200405083120-00186.warc.gz\"}"}