Mathematical Problems in Engineering
Volume 2009, Article ID 104547, 15 pages
http://dx.doi.org/10.1155/2009/104547
Research Article

## Continuation of Periodic Solutions of Dissipative and Conservative Systems: Application to Elastic Pendulum

Department of Mathematics, Prague Institute of Chemical Technology, Technicka 5, 166 28 Prague, Czech Republic

Received 2 December 2008; Revised 26 January 2009; Accepted 21 February 2009

Academic Editor: José Roberto Castilho Piqueira

Copyright © 2009 P. Pokorny. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Continuation is an efficient algorithm for finding solutions of systems of nonlinear algebraic equations where the solutions form a one-dimensional continuum. Such systems arise naturally when investigating equilibrium points and periodic solutions of ordinary differential equations with one parameter. Continuation of isolated periodic solutions of dissipative systems is a well-established technique. Less attention has been devoted to continuation of periodic solutions of conservative systems, where periodic solutions typically form a one-parameter family. To specify a single periodic solution, additional condition must be considered. However, this gives an over-determined system, which has no solution when working with approximate numerical values. We propose a simple algorithm which solves this difficulty by using singular value decomposition of the Jacobian matrix. This algorithm is applied to the conservative model of elastic pendulum. A branch of periodic solutions with constant energy is found which is born by the period doubling bifurcation of vertical oscillations.

#### 1. Motivation

When investigating the model of elastic pendulum [1], we have found a period doubling bifurcation when vertical oscillation loses its stability and a new branch of periodic solution is generated. To follow this branch of periodic solutions for varying parameter of the system and for constant energy, we developed a simple algorithm for continuation of periodic solutions in conservative systems. As we believe this algorithm will be useful for investigation of other systems as well we present the problems and solutions related to continuation of periodic solutions of conservative systems in the first half of this paper. The second half describes the results of this algorithm when applied to the model of the elastic pendulum.

Basic introduction to dynamical systems is [2]; principles of continuation can be found in [3]. To follow the second part, we invite the reader to check [1], although detailed knowledge of the second part is not necessary for understanding the algorithm.

#### 2. Continuation of Equilibrium Points

For clarity we sketch the main idea behind continuation of equilibrium points of ODE’s. A nice introduction can be found in [3]. Consider a system of ordinary differential equations with one parameter where is a smooth function of both and . A point is called an equilibrium point of the system (2.1) if , that is, is a constant solution. For a system with one parameter we want to investigate the dependence . The brute force approach would be to choose a fixed value of the parameter and then to solve a system of equations for unknowns This could be done for various values , possibly from an equidistant grid. A good candidate for the solver is Newton’s method.

Instead of solving these problems independently, it would be convenient to use the solution for the previous value of the parameter as the initial approximation for the next (close) value of the parameter.

Continuation makes use of the smoothness of even more. First, we introduce Then the equation to solve is where . We assume that the solutions form a curve in .

Continuation is a predictor—corrector method. Given one point on the curve, we first find the prediction , and then we find its correction . Then we iterate this procedure. As a result we find a series of points approximating the curve in . As the corrector, we can use Newton’s method after we add one more equation to get a system of equations for unknowns. The geometrical meaning of this additional equation is that we want to find a point on the curve such that the correction is perpendicular to the predicted displacement .

#### 3. Continuation of Periodic Solutions

##### 3.1. Previous Results

Continuation of periodic orbits of dissipative systems is a well-established technique [35]. Less attention has been devoted to continuation of periodic orbits of conservative systems. In [6] a conservative system is embedded in a one-parameter family of dissipative systems by adding a small gradient perturbation term to the vector field in such a way that a periodic orbit can only exist when the perturbation is zero. This has been generalized to systems having more than one first integrals in [7].

We use a different approach. It can be summarized as follows. To identify one periodic solution, additional equations are needed. This gives an overdetermined system, which has no solution when working with approximate values. Singular value decomposition is used to find the correction that minimizes the residue. This paper describes the principles, problems, and remedy of continuation of periodic solutions.

##### 3.2. Dissipative Systems

Before we go to conservative systems, we describe continuation of periodic solutions of dissipative systems.

First, we rescale the right-hand side of the ODE. If the equation has a periodic solution with the least period then where has a periodic solution with the least period . For other than periodic solutions may be any fixed positive constant.

Let us denote the flow of (3.2) by , that is, the solution of (3.2), satisfying the initial condition can be written as To find a periodic solution, that is, a solution satisfying we have unknowns. To use Newton's method for the correction we need to add two more equations: (which is called the phase condition) and, as stated before, where .

The geometrical meaning of (3.7) is that we want the change to be perpendicular to the tangent vector of the trajectory. This is a natural condition because any point on the periodic trajectory satisfies (3.6) and we want to remove this degeneracy.

For an isolated periodic solution, which is a typical case for dissipative ODE’s, the system of algebraic equations (3.6), (3.7), and (3.8) for unknowns can be solved by Newton’s method. Unfortunately this is not a typical case for conservative systems which we discuss in the next section.

##### 3.3. Conservative Systems
###### 3.3.1. Constant Parameter

In this subsection we assume that the parameter is fixed and we leave it out.

In a typical case, periodic solutions to a conservative system form a one parameter family. We can choose one solution by adding one more equation namely, the condition of constant total energy. However, then we get a numerically ill-conditioned system.

This may be illustrated by the simplest possible case—the harmonic oscillator The general solution is (up to a time shift) The constant (which is arbitrary) is related to the total energy which is conserved. In polar coordinates the system can be expressed as The flow of this system is (omitting and ) and the condition for periodic solution is The condition of the total energy being gives the equation and the phase condition results in Now, the system (3.16)–(3.18) has a unique solution, and there seems to be no problem. There is none indeed when we can solve it analytically and exactly. However, in most cases we have to solve it only approximately and numerically using Newton’s method.

When we rewrite the system (3.16)–(3.18) in single line where then the Jacobi matrix of partial derivatives of , used in Newton's method, is (here we differentiate with respect to resp.). When evaluated numerically, will have all entries nonzero and approximate values. Then we get an overdetermined system which has no solution in a typical case. Singular value decomposition is a good approach to such systems. It is explained in the next section.

###### 3.3.2. Singular Value Decomposition

Any matrix , whose number of rows is greater than or equal to the number of columns , can be written (see [8]) as the product of three matrices: an column-orthogonal matrix , an diagonal matrix with positive or zero elements on the diagonal (the singular values), and the transpose of an orthogonal matrix : If the system of equations has a unique solution, then it can be found as Here is the inverse to if it exists (which is the case when all the entries on the diagonal of the matrix are nonzero). In general the matrix is a diagonal matrix with entries that are computed by the following formula:

If there is no solution of (3.23) which is a typical case for an overdetermined system, then (3.24) gives the least squares solution, that is, the solution that minimizes the norm of the residue Even in cases when (3.23) has infinitely many solutions, then (3.24) gives a reasonable result: the solution with the least norm.

In numerical computation we find the singular value decomposition of the given matrix and then to find we treat values in the matrix below a certain threshold (say ) as zero. This gives a robust algorithm.

Note on dimensions of the matrices: there is another possible choice of the dimensions of the matrices in the decomposition (3.22), namely, of type , of type and of type as before. This choice is more convenient if . In our case, when , this choice means that we add zero rows to the matrix and we add columns to the matrix , but these extra columns are insignificant, because they are multiplied by zero. Up to this formal difference, the two choices are equivalent.

###### 3.3.3. Varying Parameter

Now, when we described the difficulty and the remedy with periodic solutions of conservative systems with a fixed parameter, we want to build the final algorithm for continuation of periodic solutions of conservative systems with one varying parameter.

Remember that adding one parameter in the system when searching equilibrium points, we appended the parameter to the state vector to get and we had to add the orthogonality condition (2.5), saying that the correction should by orthogonal to the prediction.

For continuation of periodic solutions of dissipative systems, we append not only the parameter but also the unknown period to get . Then we had to add not only the orthogonality condition (3.8) but also the phase condition (3.7) to get the system of equations for unknowns to solve in the corrector step (by Newton’s method).

For continuation of periodic solutions of conservative systems, we have a three-dimensional continuum of solutions of in a typical case. To pick out a single point we need three additional conditions: the orthogonality condition for the phase condition and one more extra condition. A possible choice for this extra condition is the constant energy condition where is a fixed constant, possibly depending on the parameter . A general form of this condition can be written as where is a smooth function.

With these three conditions however we get an overdetermined system which can be conveniently solved by the singular value decomposition in each corrector step.

If we write (3.28)–(3.31) in a single line as where , which stands for equations for unknowns, then in each corrector step we can solve (3.33) by Newton-like iterations where is the solution of equations for unknowns with As we have shown above, is not a regular square matrix, so we solve (3.35) by the singular value decomposition: first, we decompose into and then we find where, as described above, we compute by replacing large values on the diagonal with while replacing low values with .

###### 3.3.4. Damping

It is safer to replace (3.34) with where the damping factor is found in the following way. We start with ; we compute and . Then we compare the norms of residues; if , we are done. If not, we repeatedly decrease (by a factor, say, ) and we recompute and until we get a lower residue. If this is not achieved in a small number of iterations (say, ) then a better prediction must be done with a smaller step size.

###### 3.3.5. No Threshold Needed

Two more technical questions remain open. First, when computing large values are replaced with and low values are replaced with . How to choose the level to decide whether a value is large or low? Instead of choosing one fixed level, we use a different approach. We sort the singular values in decreasing order. First we replace the lowest singular value with zero, we compute and then and . Then we replace the two lowest singular values and compute again. This is done for all possible choices and the case giving the lowest residue is used.

###### 3.3.6. How Many Iterations

The other question is how many iterations (3.34) should be done. Initially the residue decreases with until the noise level (given by the round-off error and by the error of the methods used) is reached. Then the residue attains similar values. This suggests to compute further iterations until the residue fails to decrease.

To illustrate this algorithm, we apply it to the model of elastic pendulum in the next section.

#### 4. Elastic Pendulum

Elastic pendulum is a simple model that exhibits a wide and surprising range of dynamic phenomena [915]. The first known publication about the elastic pendulum is the paper [16] by Vitt and Gorelik. Deterministic chaos has been observed in numerical simulation and presented in the form of Poincaré section, autocorrelation function, Lyapunov exponent and power spectrum in [1719]. The bifurcation diagram of a plane elastic pendulum is sketched in [20]. The most thorough treatment of small amplitude oscillation of both plane and space elastic pendulum can be found in the works by Peter Lynch [2124].

Consider a pendulum consisting of a point bob of mass suspended on a homogeneous elastic spring of mass with the elasticity constant ; see Figure 1. Assume a homogeneous gravitational field .

Figure 1: Heavy spring elastic pendulum consists of a point mass attached to a spring with mass which is fixed at the other end point.

In [1] we derive the equations of motion This system has the only parameter where is the length of the unloaded spring.

In [1] we also derive the stability condition for the vertical oscillations where is the amplitude of the vertical oscillation given by initial conditions. The vertical oscillation of the heavy spring elastic pendulum with the relative amplitude is stable (see [1]) for where The error of this estimate of is less than for .

Figure 2 shows the regions in the - plane with stable (light gray, or blue in color version) and unstable (dark gray, or red in color version) vertical oscillation with amplitude and the parameter . For and sufficiently small, the vertical oscillations are stable. For a fixed let us observe what happens when we increase . When crosses the critical value the vertical oscillations (with the period ) lose their stability, one eigenvalue of the Jacobi matrix leaves the unit circle through the point in the complex plane, and a period doubling bifurcation gives rise to a family of periodic oscillations (with the period close to ). Another period doubling bifurcation exists when decreasing the parameter for the critical value . We want to find the branch of these periodic oscillations.

Figure 2: The stable region (light gray, or blue in color version) and the unstable region (dark gray, or red in color version) for the vertical oscillation with amplitude and the external field .

As there is a continuum of periodic solutions for one value of the parameter, we want to pick out one of these periodic solutions by specifying the energy. The total energy of the system consists of the potential and of the kinetic energy where We are interested in periodic solutions with the total energy where is the energy of the vertical oscillations (4.3):

To find the initial approximation of one periodic solution, we choose a parameter value and Then we can use numerical integration and Poincaré section given by the horizontal plane going through the only stable equilibrium point of the system; see Figure 3. For various initial conditions we get various closed invariant curves densely filled by the intersections of the trajectory with the Poincaré section. From this picture we can estimate the coordinates and . From the position of the Poincaré section we know that and we compute from the known total energy. We estimate the period to be because the period of vertical oscillations before the period doubling bifurcation is . Finally we choose the initial prediction vector (determining the direction in which the curve will be traveled) meaning that we want to start with decreasing parameter . Now we can start the continuation. The results are discussed in the next section.

Figure 3: Poincaré section can be used to find the initial approximation of the periodic orbit, which is represented by an isolated point while surrounding closed curves correspond to invariant tori.

#### 5. Results

In this section we illustrate the above algorithm for continuation of periodic solutions by the Hamiltonian model of elastic pendulum. We choose the fixed energy corresponding to the amplitude of vertical oscillation with amplitude . When decreasing the parameter , the vertical oscillation becomes unstable for the critical value . Period doubling bifurcation gives rise to a one parameter family of periodic solutions that exist down to the parameter value . We use the above algorithm for continuation of this branch of periodic solutions.

Figure 4 shows how the period of the periodic solution grows with decreasing parameter from the value just after the period doubling bifurcation up to the value .

Figure 4: The period of periodic oscillation (as a function of the parameter ) changes from for high up to for low .

Figure 5 shows for the same parameter values how the range of coordinates (the horizontal displacement of the bob) grows and then decreases down to zero.

Figure 5: The range of periodic oscillation as a function of the parameter .

When decreasing the parameter the shape of the trajectory in the - plane changes in an interesting way. This can be illustrated by plotting the trajectory for a few parameter values see Figure 6. First, as the vertical oscillation with period loses its stability a new periodic solution is born with twice the period and a small range of . As the parameter goes down, the range of grows but the cup-like shape of the oscillation remains. For a still decreasing parameter the top parts of the trajectory begin to bend and then they turn down. The period grows. For parameter values close to the trajectory is four times folded and the range of decreases to zero. The period is now close to which is four times the period of vertical oscillation. Finally our periodic solution merges with the vertical oscillation.

Figure 6: The projection of the trajectory to the - plane for four values of the decreasing parameter . Top left: For slightly below the critical value the period doubling bifurcation gives rise to the new oscillation with twice the period and with a small range. Top right: For decreasing the cup-like shape of the trajectory widens and the top parts start to bend. Bottom left: For still lower values the outer parts of the trajectory bend even more and the period grows. Bottom right: Finally the period is close to which is four times the period of the vertical oscillation. For our periodic oscillation merges with the vertical oscillation.

#### 6. Conclusion

The difference between continuation of periodic solutions of dissipative and conservative systems was described. For dissipative systems periodic solutions are generically isolated. To specify a single point on the trajectory we need to add the orthogonality condition and the phase condition. Then continuation of periodic solutions can be done by the predictor—corrector algorithm. For conservative systems, however, periodic solutions form one-dimensional continuum parametrized by the energy. To specify a single solution we must add one more equation. When working with approximate numerical values, this yields an overdetermined system in the corrector step. To solve this problem we suggest to use singular value decomposition. This gives a simple algorithm that was tested by the Hamiltonian model of elastic pendulum. Using stability analysis from [1] we knew that a period doubling bifurcation exists for this model for a known parameter value. In this bifurcation the vertical oscillation loses its stability and a branch of new periodic solutions is generated. In the present paper we used the suggested algorithm to compute this new branch of periodic solutions for varying parameter . It turns out that the new periodic solution starts with twice the period of the vertical oscillation, it gets wider, it bends and folds, its period grows to which is four times the period of the vertical oscillation, and finally the new oscillation merges with the vertical oscillation again.

#### Acknowledgments

This work is supported by the project MSM 6046137306. The access to the METACentrum supercomputing facilities MSM 6383917201 is highly appreciated.

#### References

1. P. Pokorny, “Stability condition for vertical oscillation of 3-dim heavy spring elastic pendulum,” Regular and Chaotic Dynamics, vol. 13, no. 3, pp. 155–165, 2008.
2. M. W. Hirsch, S. Smale, and R. L. Devaney, Differential Equations, Dynamical Systems, and an Introduction to Chaos, vol. 60 of Pure and Applied Mathematics, Elsevier/Academic Press, Amsterdam, The Netherlands, 2nd edition, 2004.
3. R. Seydel, “Tutorial on continuation,” International Journal of Bifurcation and Chaos, vol. 1, no. 1, pp. 3–11, 1991.
4. E. J. Doedel, R. C. Paffenroth, A. R. Champneys et al., “AUTO2000: continuation and bifurcation software for ordinary differential equations,” Tech. Rep., Department of Computer Science, Concordia University, Montreal, Canada, 2000.
5. M. Marek and I. Schreiber, Chaotic Behaviour of Deterministic Dissipative Systems, Cambridge University Press, Cambridge, UK, 1991.
6. J.-A. Sepulchre and R. S. MacKay, “Localized oscillations in conservative or dissipative networks of weakly coupled autonomous oscillators,” Nonlinearity, vol. 10, no. 3, pp. 679–713, 1997.
7. F. J. Muñoz-Almaraz, E. Freire, J. Galán, E. Doedel, and A. Vanderbauwhede, “Continuation of periodic orbits in conservative and Hamiltonian systems,” Physica D, vol. 181, no. 1-2, pp. 1–38, 2003.
8. W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, Cambridge University Press, Cambridge, UK, 2nd edition, 1992.
9. M. G. Olsson, “Why does a mass on a spring sometimes misbehave?” American Journal of Physics, vol. 44, no. 12, pp. 1211–1212, 1976.
10. T. E. Cayton, “The laboratory spring-mass oscillator: an example of parametric instability,” American Journal of Physics, vol. 45, no. 8, pp. 723–732, 1977.
11. M. G. Rusbridge, “Motion of the sprung pendulum,” American Journal of Physics, vol. 48, no. 2, pp. 146–151, 1980.
12. E. Breitenberger and R. D. Mueller, “The elastic pendulum: a nonlinear paradigm,” Journal of Mathematical Physics, vol. 22, no. 6, pp. 1196–1210, 1981.
13. H. M. Lai, “On the recurrence phenomenon of a resonant spring pendulum,” American Journal of Physics, vol. 52, no. 3, pp. 219–223, 1984.
14. B. A. Anicin, D. M. Davidovic, and V. M. Babovic, “On the linear theory of the elastic pendulum,” European Journal of Physics, vol. 14, no. 3, pp. 132–135, 1993.
15. D. M. Davidović, B. A. Aničin, and V. M. Babović, “The libration limits of the elastic pendulum,” American Journal of Physics, vol. 64, no. 3, pp. 338–342, 1996.
16. A. Vitt and G. Gorelik, “Oscillations of an elastic pendulum as an example of the oscillations of two parametrically coupled linear systems,” Journal of Technical Physics, vol. 3, no. 2-3, pp. 294–307, 1933, translated from Russian by Lisa Shields with an Introduction by Peter Lynch.
17. H. N. Núñez-Yépez, A. L. Salas-Brito, C. A. Vargas, and L. Vicente, “Onset of chaos in an extensible pendulum,” Physics Letters A, vol. 145, no. 2-3, pp. 101–105, 1990.
18. R. Cuerno, A. F. Rañada, and J. J. Ruiz-Lorenzo, “Deterministic chaos in the elastic pendulum: a simple laboratory for nonlinear dynamics,” American Journal of Physics, vol. 60, no. 1, pp. 73–79, 1992.
19. R. Carretero-Gonzalez, H. N. Nunez-Yepez, and A. L. Salas-Brito, “Regular and chaotic behaviour in an extensible pendulum,” European Journal of Physics, vol. 15, no. 3, pp. 139–148, 1994.
20. S. V. Kuznetsov, “The motion of the elastic pendulum,” Regular and Chaotic Dynamics, vol. 4, no. 3, pp. 3–12, 1999.
21. P. Lynch, “Resonant motions of the three-dimensional elastic pendulum,” International Journal of Non-Linear Mechanics, vol. 37, no. 2, pp. 345–367, 2002.
22. P. Lynch, “The swinging spring: a simple model of atmospheric balance,” in Large-Scale Atmosphere-Ocean Dynamics, Vol. II, J. Norbury and I. Roulstone, Eds., pp. 64–108, Cambridge University Press, Cambridge, UK, 2002.
23. D. D. Holm and P. Lynch, “Stepwise precession of the resonant swinging spring,” SIAM Journal on Applied Dynamical Systems, vol. 1, no. 1, pp. 44–64, 2002.
24. P. Lynch and C. Houghton, “Pulsation and precession of the resonant swinging spring,” Physica D, vol. 190, no. 1-2, pp. 38–62, 2004.