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

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 𝑑𝑥𝑑𝑡=𝑓(𝑥,𝑝),𝑥∈𝑅^𝑛,𝑝∈𝑅, (2.1) where 𝑓 is a smooth function of both 𝑥 and 𝑝. A point 𝑥_𝐸 is called an equilibrium point of the system (2.1) if 𝑓(𝑥_𝐸,𝑝)=0, 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 𝑓(𝑥,𝑝)=0. (2.2) 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 𝑦=(𝑥,𝑝)∈𝑅^𝑛+1,𝐹(𝑦)=𝑓(𝑥,𝑝). (2.3) Then the equation to solve is 𝐹(𝑦)=0, (2.4) where 𝐹∶𝑅^𝑛+1→𝑅^𝑛. We assume that the solutions form a curve in 𝑅^𝑛+1.

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 𝑅^𝑛+1. As the corrector, we can use Newton’s method after we add one more equation (𝑦−𝑦_𝑝)⋅(𝑦_𝑝−𝑦_𝑜)=0 (2.5) to get a system of 𝑛+1 equations for 𝑛+1 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 [3–5]. 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 𝑑𝑥𝑑𝑡=𝑓(𝑥,𝑝),𝑥∈𝑅^𝑛,𝑝∈𝑅 (3.1) has a periodic solution with the least period 𝑇>0, then 𝑑𝑥𝑑𝑡=𝑔(𝑥,𝑝,𝑇), (3.2) where 𝑔(𝑥,𝑝,𝑇)=𝑇𝑓(𝑥,𝑝), (3.3) has a periodic solution with the least period 1. 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 𝑥(0)=𝑥_ic, (3.4) can be written as 𝑥(𝑡)=𝜙(𝑡,𝑥_ic,𝑝,𝑇). (3.5) To find a periodic solution, that is, a solution satisfying 𝜙(1,𝑥,𝑝,𝑇)−𝑥=0, (3.6) we have 𝑛+2 unknowns. To use Newton's method for the correction we need to add two more equations: (𝑥−𝑥_𝑜)⋅𝑔(𝑥_𝑜,𝑝,𝑇)=0 (3.7) (which is called the phase condition) and, as stated before, (𝑦−𝑦_𝑝)⋅(𝑦_𝑝−𝑦_𝑜)=0, (3.8) where 𝑦=(𝑥,𝑝,𝑇)∈𝑅^𝑛+2.

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 𝑛+2 algebraic equations (3.6), (3.7), and (3.8) for 𝑛+2 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 𝐸(𝑥,𝑥̇)=𝐸₀, (3.9) 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 𝑥̈+𝑥=0,𝑥∈𝑅. (3.10) The general solution is (up to a time shift) 𝑥=𝑎cos𝑡. (3.11) The constant 𝑎 (which is arbitrary) is related to the total energy ̇𝑥𝐸(𝑥,𝑥)=²2+𝑥̇²2=𝑎²2, (3.12) which is conserved. In polar coordinates 𝑟,𝜑 the system can be expressed as 𝑟̇=0,𝜑̇=1. (3.13) The flow of this system is (omitting 𝑝 and 𝑇) 𝜙(𝑡,𝑟₀,𝜑₀)=(𝑟₀,𝜑₀+𝑡), (3.14) and the condition for periodic solution 𝜙(𝑇,𝑟,𝜑)−(𝑟,𝜑)=(0,2𝜋) (3.15) is 0=0,𝑇=2𝜋. (3.16) The condition of the total energy being 𝐸₀ gives the equation 𝑟=2𝐸₀, (3.17) and the phase condition results in 𝜑=𝜑₀. (3.18) 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 ℱ(𝑇,𝑟,𝜑)=(0,0,0,0), (3.19) where ℱ(𝑇,𝑟,𝜑)=(0,𝑇−2𝜋,𝑟−2𝐸₀,𝜑−𝜑₀), (3.20) then the Jacobi matrix of partial derivatives of ℱ, used in Newton's method, is 𝐽=000100010001 (3.21) (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 𝑉: 𝐴=𝑈⋅𝑊⋅𝑉^𝑇. (3.22) If the system of equations 𝐴⋅𝑥=𝑏 (3.23) has a unique solution, then it can be found as 𝑥=𝑉⋅𝑊⁻¹⋅𝑈^𝑇⋅𝑏. (3.24) 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: (𝑊⁻¹)_𝑖𝑖=1𝑤_𝑖𝑖,if𝑤_𝑖𝑖≠00,if𝑤_𝑖𝑖=0. (3.25)

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 𝐴⋅𝑥−𝑏. (3.26) 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 10⁻⁶) 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 𝑓(𝑥,𝑝)=0 (3.27) when searching equilibrium points, we appended the parameter 𝑝 to the state vector 𝑥 to get 𝑦=(𝑥,𝑝)∈𝑅^𝑛+1, 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 𝑦=(𝑥,𝑝,𝑇)∈𝑅^𝑛+2. Then we had to add not only the orthogonality condition (3.8) but also the phase condition (3.7) to get the system of 𝑛+2 equations for 𝑛+2 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 𝜙(1,𝑥,𝑝,𝑇)−𝑥=0 (3.28) in a typical case. To pick out a single point we need three additional conditions: the orthogonality condition for 𝑦=(𝑥,𝑝,𝑇)(𝑦−𝑦_𝑝)⋅(𝑦_𝑝−𝑦₀)=0, (3.29) the phase condition (𝑥−𝑥₀)⋅𝑔(𝑥₀,𝑝,𝑇)=0, (3.30) and one more extra condition. A possible choice for this extra condition is the constant energy condition 𝐸(𝑦)−𝐸₀=0, (3.31) where 𝐸₀ is a fixed constant, possibly depending on the parameter 𝑝. A general form of this condition can be written as 𝜓(𝐸(𝑦),𝑝)=0, (3.32) 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 ℱ(𝑦)=0, (3.33) where ℱ∶𝑅^𝑛+2→𝑅^𝑛+3, which stands for 𝑛+3 equations for 𝑛+2 unknowns, then in each corrector step we can solve (3.33) by Newton-like iterations 𝑦_𝑘+1=𝑦_𝑘−Δ𝑦_𝑘, (3.34) where Δ𝑦_𝑘 is the solution of 𝑛+3 equations for 𝑛+2 unknowns 𝐽(𝑦_𝑘)⋅Δ𝑦_𝑘=ℱ(𝑦_𝑘) (3.35) with 𝐽_𝑖𝑗=𝜕ℱ_𝑖𝜕𝑦_𝑗. (3.36) 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 𝐽=𝑈⋅𝑊⋅𝑉^𝑇 (3.37) and then we find Δ𝑦_𝑘=𝑉⋅𝑊⁻¹⋅𝑈^𝑇⋅ℱ(𝑦_𝑘), (3.38) where, as described above, we compute 𝑊⁻¹ by replacing large values 𝑤_𝑖𝑖 on the diagonal with 𝑤⁻¹_𝑖𝑖 while replacing low values with 0.

3.3.4. Damping

It is safer to replace (3.34) with 𝑦_𝑘+1=𝑦_𝑘−𝑞Δ𝑦_𝑘, (3.39) where the damping factor 𝑞 is found in the following way. We start with 𝑞=1; we compute 𝑦_𝑘+1 and ℱ(𝑦_𝑘+1). Then we compare the norms of residues; if |ℱ(𝑦_𝑘+1)|<|𝐹(𝑦_𝑘)|, we are done. If not, we repeatedly decrease 𝑞 (by a factor, say, 2) and we recompute 𝑦_𝑘+1 and ℱ(𝑦_𝑘+1) until we get a lower residue. If this is not achieved in a small number of iterations (say, 20) 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 0. 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 𝑦_𝑘+1 and |ℱ(𝑦_𝑘+1)|. Then we replace the two lowest singular values and compute |ℱ(𝑦_𝑘+1)| again. This is done for all possible choices and the case giving the lowest residue |ℱ(𝑦_𝑘+1)| is used.

3.3.6. How Many Iterations

The other question is how many iterations (3.34) should be done. Initially the residue |ℱ(𝑦_𝑘+1)| 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 [9–15]. 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 [17–19]. 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 [21–24].

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 (0,−𝑔).

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 𝑥̈₁1=(𝑥²₁+𝑥²₂−1)𝑥₁,𝑥̈₂1=(𝑥²₁+𝑥²₂−1)𝑥₂−𝑝. (4.1) This system has the only parameter 𝑝=(𝑚_𝐵+𝑚_𝑆/2)𝑔𝑘ℓ₀, (4.2) where ℓ₀ is the length of the unloaded spring.

In [1] we also derive the stability condition for the vertical oscillations 𝑥₁=0,𝑥₂=−1−𝑝+𝑎sin𝑡, (4.3) 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 |𝑎|<|𝑎_𝐶(𝑝)|, (4.4) where 𝑎_𝐶(𝑝)≐(3𝑝)^2/3−1. (4.5) The error of this estimate of 𝑎_𝐶(𝑝) is less than 0.005 for 0<𝑝<0.6.

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 𝑎=0.1 let us observe what happens when we increase 𝑝. When 𝑝 crosses the critical value 𝑝₁≐(1/3)(1−𝑎)^3/2≐0.2846 the vertical oscillations (with the period 2𝜋) lose their stability, one eigenvalue of the Jacobi matrix leaves the unit circle through the point −1 in the complex plane, and a period doubling bifurcation gives rise to a family of periodic oscillations (with the period close to 4𝜋). Another period doubling bifurcation exists when decreasing the parameter 𝑝 for the critical value 𝑝₂≐(1/3)(1+𝑎)^3/2≐0.3846. 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 𝐸(𝑥₁,𝑥₂,𝑥̇₁,𝑥̇₂)=𝐸_pot(𝑥₁,𝑥₂)+𝐸_kin(𝑥̇₁,𝑥̇₂), (4.6) where 𝐸_pot=12(𝑥²₁+𝑥²₂−1)²+𝑝𝑥₂,𝐸_kin=12(𝑥̇²₁+𝑥̇²₂). (4.7) We are interested in periodic solutions with the total energy 𝐸(𝑥₁,𝑥₂,𝑥̇₁,𝑥̇₂)=𝐸₀, (4.8) where 𝐸₀ is the energy of the vertical oscillations (4.3): 𝐸₀1=𝐸(0,−1−𝑝+𝑎sin𝑡,0,𝑎cos𝑡)=2(𝑎²−𝑝²−2𝑝). (4.9)

To find the initial approximation of one periodic solution, we choose a parameter value 𝑝<𝑝₂ and 𝑎=0.1 Then we can use numerical integration and Poincaré section given by the horizontal plane 𝑥₂=−1−𝑝 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 𝑥₂=−1−𝑝 and we compute 𝑥̇₂ from the known total energy. We estimate the period to be 𝑇=4𝜋 because the period of vertical oscillations before the period doubling bifurcation is 2𝜋. Finally we choose the initial prediction vector (determining the direction in which the curve will be traveled) 𝑣₀=(0,0,0,0,−1,0) 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 𝑎=0.1. When decreasing the parameter 𝑝, the vertical oscillation becomes unstable for the critical value 𝑝₂≐0.3846. Period doubling bifurcation gives rise to a one parameter family of periodic solutions that exist down to the parameter value 𝑝₃≐0.0655. 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 𝑇=4𝜋 just after the period doubling bifurcation up to the value 𝑇=8𝜋.

Figure 4 

The period 𝑇 of periodic oscillation (as a function of the parameter 𝑝) changes from 4𝜋 for high 𝑝 up to 8𝜋 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 Δ𝑥1=𝑥1max−𝑥1min 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 2𝜋 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 𝑝₃≐0.0655 the trajectory is four times folded and the range of 𝑥₁ decreases to zero. The period is now close to 8𝜋 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 𝑥1-𝑥2 plane for four values of the decreasing parameter 𝑝. Top left: For 𝑝, slightly below the critical value 𝑝2≐0.3846 the period doubling bifurcation gives rise to the new oscillation with twice the period and with a small 𝑥1 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 8𝜋 which is four times the period of the vertical oscillation. For 𝑝=𝑝3≐0.0655 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 8𝜋 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.