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

Authors: P. Pokorny
Journal: Mathematical Problems in Engineering
 (2009)

---

## 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.


---

## Body

## 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)dxdt=f(x,p),x∈Rn,p∈R,
where f is a smooth function of both x and p. A point xE is called an equilibrium point of the system (2.1) if f(xE,p)=0, that is, x(t)=xE is a constant solution. For a system with one parameter p we want to investigate the dependence xE(p). The brute force approach would be to choose a fixed value p̃ of the parameter p and then to solve a system of n equations for n unknowns
(2.2)f(x,p̃)=0.
This could be done for various values p̃, 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 off even more. First, we introduce
(2.3)y=(x,p)∈Rn+1,F(y)=f(x,p).
Then the equation to solve is
(2.4)F(y)=0,
where F:Rn+1→Rn. We assume that the solutions form a curve in Rn+1.Continuation is a predictor—corrector method. Given one pointyo on the curve, we first find the prediction yp, and then we find its correction yc. Then we iterate this procedure. As a result we find a series of points approximating the curve in Rn+1. As the corrector, we can use Newton’s method after we add one more equation
(2.5)(y-yp)·(yp-yo)=0
to get a system of n+1 equations for n+1 unknowns. The geometrical meaning of this additional equation is that we want to find a point y on the curve such that the correction y-yp is perpendicular to the predicted displacement yp-yo.


## 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)dxdt=f(x,p),x∈Rn,p∈R
has a periodic solution with the least period T>0, then
(3.2)dxdt=g(x,p,T),
where
(3.3)g(x,p,T)=Tf(x,p),
has a periodic solution with the least period 1. For other than periodic solutions T may be any fixed positive constant.Let us denote the flow of (3.2) by ϕ(t,x,p,T), that is, the solution x(t) of (3.2), satisfying the initial condition
(3.4)x(0)=xic,
can be written as
(3.5)x(t)=ϕ(t,xic,p,T).
To find a periodic solution, that is, a solution satisfying
(3.6)ϕ(1,x,p,T)-x=0,
we have n+2 unknowns. To use Newton's method for the correction we need to add two more equations:
(3.7)(x-xo)·g(xo,p,T)=0
(which is called the phase condition) and, as stated before,
(3.8)(y-yp)·(yp-yo)=0,
where y=(x,p,T)∈Rn+2.The geometrical meaning of (3.7) is that we want the change x-xo 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 ofn+2 algebraic equations (3.6), (3.7), and (3.8) for n+2 unknowns x,p,T 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 parameterp 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)E(x,ẋ)=E0,
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(3.10)ẍ+x=0,x∈R.
The general solution is (up to a time shift)
(3.11)x=acost.
The constant a (which is arbitrary) is related to the total energy
(3.12)E(x,ẋ)=x22+ẋ22=a22,
which is conserved. In polar coordinates r,φ the system can be expressed as
(3.13)ṙ=0,φ̇=1.
The flow of this system is (omitting p and T)
(3.14)ϕ(t,r0,φ0)=(r0,φ0+t),
and the condition for periodic solution
(3.15)ϕ(T,r,φ)-(r,φ)=(0,2π)
is
(3.16)0=0,T=2π.
The condition of the total energy being E0 gives the equation
(3.17)r=2E0,
and the phase condition results in
(3.18)φ=φ0.
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
(3.19)ℱ(T,r,φ)=(0,0,0,0),
where
(3.20)ℱ(T,r,φ)=(0,T-2π,r-2E0,φ-φ0),
then the Jacobi matrix of partial derivatives of ℱ, used in Newton's method, is
(3.21)J=(000100010001)
(here we differentiate with respect to T,r,φ, resp.). When evaluated numerically, J 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

Anym×n matrix A, whose number of rows m is greater than or equal to the number of columns n, can be written (see [8]) as the product of three matrices: an m×n column-orthogonal matrix U, an n×n diagonal matrix W with positive or zero elements on the diagonal (the singular values), and the transpose of an n×n orthogonal matrix V:
(3.22)A=U·W·VT.
If the system of equations
(3.23)A·x=b
has a unique solution, then it can be found as
(3.24)x=V·W-1·UT·b.
Here W-1 is the inverse to W if it exists (which is the case when all the entries wii on the diagonal of the matrix W are nonzero). In general the matrix W-1 is a diagonal matrix with entries that are computed by the following formula:
(3.25)(W-1)ii={1wii,ifwii≠00,ifwii=0.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)A·x-b.
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 matrixA, and then to find W-1 we treat values wii in the matrix W below a certain threshold (say 10-6) 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, U of type m×m, W of type m×n and V of type n×n as before. This choice is more convenient if m<n. In our case, when m>n, this choice means that we add m-n zero rows to the matrix W and we add m-n columns to the matrix U, 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(3.27)f(x,p)=0
when searching equilibrium points, we appended the parameter p to the state vector x to get y=(x,p)∈Rn+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 parameterp but also the unknown period T to get y=(x,p,T)∈Rn+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 n+2 equations for n+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(3.28)ϕ(1,x,p,T)-x=0
in a typical case. To pick out a single point we need three additional conditions: the orthogonality condition for y=(x,p,T)(3.29)(y-yp)·(yp-y0)=0,
the phase condition
(3.30)(x-x0)·g(x0,p,T)=0,
and one more extra condition. A possible choice for this extra condition is the constant energy condition
(3.31)E(y)-E0=0,
where E0 is a fixed constant, possibly depending on the parameter p. A general form of this condition can be written as
(3.32)ψ(E(y),p)=0,
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
(3.33)ℱ(y)=0,
where ℱ:Rn+2→Rn+3, which stands for n+3 equations for n+2 unknowns, then in each corrector step we can solve (3.33) by Newton-like iterations
(3.34)yk+1=yk-Δyk,
where Δyk is the solution of n+3 equations for n+2 unknowns
(3.35)J(yk)·Δyk=ℱ(yk)
with
(3.36)Jij=∂ℱi∂yj.
As we have shown above, J is not a regular square matrix, so we solve (3.35) by the singular value decomposition: first, we decompose J into
(3.37)J=U·W·VT
and then we find
(3.38)Δyk=V·W-1·UT·ℱ(yk),
where, as described above, we compute W-1 by replacing large values wii on the diagonal with wii-1 while replacing low values with 0.


#### 3.3.4. Damping

It is safer to replace (3.34) with
(3.39)yk+1=yk-qΔyk,
where the damping factor q is found in the following way. We start with q=1; we compute yk+1 and ℱ(yk+1). Then we compare the norms of residues; if |ℱ(yk+1)|<|F(yk)|, we are done. If not, we repeatedly decrease q (by a factor, say, 2) and we recompute yk+1 and ℱ(yk+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 computingW-1, large values wii are replaced with wii-1 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 wii in decreasing order. First we replace the lowest singular value with zero, we compute Δyk and then yk+1 and |ℱ(yk+1)|. Then we replace the two lowest singular values and compute |ℱ(yk+1)| again. This is done for all possible choices and the case giving the lowest residue |ℱ(yk+1)| is used.


#### 3.3.6. How Many Iterations

The other question is how many iterations (3.34) should be done. Initially the residue |ℱ(yk+1)| decreases with k 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.


## 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)dxdt=f(x,p),x∈Rn,p∈R
has a periodic solution with the least period T>0, then
(3.2)dxdt=g(x,p,T),
where
(3.3)g(x,p,T)=Tf(x,p),
has a periodic solution with the least period 1. For other than periodic solutions T may be any fixed positive constant.Let us denote the flow of (3.2) by ϕ(t,x,p,T), that is, the solution x(t) of (3.2), satisfying the initial condition
(3.4)x(0)=xic,
can be written as
(3.5)x(t)=ϕ(t,xic,p,T).
To find a periodic solution, that is, a solution satisfying
(3.6)ϕ(1,x,p,T)-x=0,
we have n+2 unknowns. To use Newton's method for the correction we need to add two more equations:
(3.7)(x-xo)·g(xo,p,T)=0
(which is called the phase condition) and, as stated before,
(3.8)(y-yp)·(yp-yo)=0,
where y=(x,p,T)∈Rn+2.The geometrical meaning of (3.7) is that we want the change x-xo 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 ofn+2 algebraic equations (3.6), (3.7), and (3.8) for n+2 unknowns x,p,T 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 parameterp 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)E(x,ẋ)=E0,
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(3.10)ẍ+x=0,x∈R.
The general solution is (up to a time shift)
(3.11)x=acost.
The constant a (which is arbitrary) is related to the total energy
(3.12)E(x,ẋ)=x22+ẋ22=a22,
which is conserved. In polar coordinates r,φ the system can be expressed as
(3.13)ṙ=0,φ̇=1.
The flow of this system is (omitting p and T)
(3.14)ϕ(t,r0,φ0)=(r0,φ0+t),
and the condition for periodic solution
(3.15)ϕ(T,r,φ)-(r,φ)=(0,2π)
is
(3.16)0=0,T=2π.
The condition of the total energy being E0 gives the equation
(3.17)r=2E0,
and the phase condition results in
(3.18)φ=φ0.
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
(3.19)ℱ(T,r,φ)=(0,0,0,0),
where
(3.20)ℱ(T,r,φ)=(0,T-2π,r-2E0,φ-φ0),
then the Jacobi matrix of partial derivatives of ℱ, used in Newton's method, is
(3.21)J=(000100010001)
(here we differentiate with respect to T,r,φ, resp.). When evaluated numerically, J 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

Anym×n matrix A, whose number of rows m is greater than or equal to the number of columns n, can be written (see [8]) as the product of three matrices: an m×n column-orthogonal matrix U, an n×n diagonal matrix W with positive or zero elements on the diagonal (the singular values), and the transpose of an n×n orthogonal matrix V:
(3.22)A=U·W·VT.
If the system of equations
(3.23)A·x=b
has a unique solution, then it can be found as
(3.24)x=V·W-1·UT·b.
Here W-1 is the inverse to W if it exists (which is the case when all the entries wii on the diagonal of the matrix W are nonzero). In general the matrix W-1 is a diagonal matrix with entries that are computed by the following formula:
(3.25)(W-1)ii={1wii,ifwii≠00,ifwii=0.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)A·x-b.
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 matrixA, and then to find W-1 we treat values wii in the matrix W below a certain threshold (say 10-6) 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, U of type m×m, W of type m×n and V of type n×n as before. This choice is more convenient if m<n. In our case, when m>n, this choice means that we add m-n zero rows to the matrix W and we add m-n columns to the matrix U, 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(3.27)f(x,p)=0
when searching equilibrium points, we appended the parameter p to the state vector x to get y=(x,p)∈Rn+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 parameterp but also the unknown period T to get y=(x,p,T)∈Rn+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 n+2 equations for n+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(3.28)ϕ(1,x,p,T)-x=0
in a typical case. To pick out a single point we need three additional conditions: the orthogonality condition for y=(x,p,T)(3.29)(y-yp)·(yp-y0)=0,
the phase condition
(3.30)(x-x0)·g(x0,p,T)=0,
and one more extra condition. A possible choice for this extra condition is the constant energy condition
(3.31)E(y)-E0=0,
where E0 is a fixed constant, possibly depending on the parameter p. A general form of this condition can be written as
(3.32)ψ(E(y),p)=0,
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
(3.33)ℱ(y)=0,
where ℱ:Rn+2→Rn+3, which stands for n+3 equations for n+2 unknowns, then in each corrector step we can solve (3.33) by Newton-like iterations
(3.34)yk+1=yk-Δyk,
where Δyk is the solution of n+3 equations for n+2 unknowns
(3.35)J(yk)·Δyk=ℱ(yk)
with
(3.36)Jij=∂ℱi∂yj.
As we have shown above, J is not a regular square matrix, so we solve (3.35) by the singular value decomposition: first, we decompose J into
(3.37)J=U·W·VT
and then we find
(3.38)Δyk=V·W-1·UT·ℱ(yk),
where, as described above, we compute W-1 by replacing large values wii on the diagonal with wii-1 while replacing low values with 0.


### 3.3.4. Damping

It is safer to replace (3.34) with
(3.39)yk+1=yk-qΔyk,
where the damping factor q is found in the following way. We start with q=1; we compute yk+1 and ℱ(yk+1). Then we compare the norms of residues; if |ℱ(yk+1)|<|F(yk)|, we are done. If not, we repeatedly decrease q (by a factor, say, 2) and we recompute yk+1 and ℱ(yk+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 computingW-1, large values wii are replaced with wii-1 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 wii in decreasing order. First we replace the lowest singular value with zero, we compute Δyk and then yk+1 and |ℱ(yk+1)|. Then we replace the two lowest singular values and compute |ℱ(yk+1)| again. This is done for all possible choices and the case giving the lowest residue |ℱ(yk+1)| is used.


### 3.3.6. How Many Iterations

The other question is how many iterations (3.34) should be done. Initially the residue |ℱ(yk+1)| decreases with k 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.


## 3.3.1. Constant Parameter

In this subsection we assume that the parameterp 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)E(x,ẋ)=E0,
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(3.10)ẍ+x=0,x∈R.
The general solution is (up to a time shift)
(3.11)x=acost.
The constant a (which is arbitrary) is related to the total energy
(3.12)E(x,ẋ)=x22+ẋ22=a22,
which is conserved. In polar coordinates r,φ the system can be expressed as
(3.13)ṙ=0,φ̇=1.
The flow of this system is (omitting p and T)
(3.14)ϕ(t,r0,φ0)=(r0,φ0+t),
and the condition for periodic solution
(3.15)ϕ(T,r,φ)-(r,φ)=(0,2π)
is
(3.16)0=0,T=2π.
The condition of the total energy being E0 gives the equation
(3.17)r=2E0,
and the phase condition results in
(3.18)φ=φ0.
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
(3.19)ℱ(T,r,φ)=(0,0,0,0),
where
(3.20)ℱ(T,r,φ)=(0,T-2π,r-2E0,φ-φ0),
then the Jacobi matrix of partial derivatives of ℱ, used in Newton's method, is
(3.21)J=(000100010001)
(here we differentiate with respect to T,r,φ, resp.). When evaluated numerically, J 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

Anym×n matrix A, whose number of rows m is greater than or equal to the number of columns n, can be written (see [8]) as the product of three matrices: an m×n column-orthogonal matrix U, an n×n diagonal matrix W with positive or zero elements on the diagonal (the singular values), and the transpose of an n×n orthogonal matrix V:
(3.22)A=U·W·VT.
If the system of equations
(3.23)A·x=b
has a unique solution, then it can be found as
(3.24)x=V·W-1·UT·b.
Here W-1 is the inverse to W if it exists (which is the case when all the entries wii on the diagonal of the matrix W are nonzero). In general the matrix W-1 is a diagonal matrix with entries that are computed by the following formula:
(3.25)(W-1)ii={1wii,ifwii≠00,ifwii=0.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)A·x-b.
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 matrixA, and then to find W-1 we treat values wii in the matrix W below a certain threshold (say 10-6) 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, U of type m×m, W of type m×n and V of type n×n as before. This choice is more convenient if m<n. In our case, when m>n, this choice means that we add m-n zero rows to the matrix W and we add m-n columns to the matrix U, 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(3.27)f(x,p)=0
when searching equilibrium points, we appended the parameter p to the state vector x to get y=(x,p)∈Rn+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 parameterp but also the unknown period T to get y=(x,p,T)∈Rn+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 n+2 equations for n+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(3.28)ϕ(1,x,p,T)-x=0
in a typical case. To pick out a single point we need three additional conditions: the orthogonality condition for y=(x,p,T)(3.29)(y-yp)·(yp-y0)=0,
the phase condition
(3.30)(x-x0)·g(x0,p,T)=0,
and one more extra condition. A possible choice for this extra condition is the constant energy condition
(3.31)E(y)-E0=0,
where E0 is a fixed constant, possibly depending on the parameter p. A general form of this condition can be written as
(3.32)ψ(E(y),p)=0,
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
(3.33)ℱ(y)=0,
where ℱ:Rn+2→Rn+3, which stands for n+3 equations for n+2 unknowns, then in each corrector step we can solve (3.33) by Newton-like iterations
(3.34)yk+1=yk-Δyk,
where Δyk is the solution of n+3 equations for n+2 unknowns
(3.35)J(yk)·Δyk=ℱ(yk)
with
(3.36)Jij=∂ℱi∂yj.
As we have shown above, J is not a regular square matrix, so we solve (3.35) by the singular value decomposition: first, we decompose J into
(3.37)J=U·W·VT
and then we find
(3.38)Δyk=V·W-1·UT·ℱ(yk),
where, as described above, we compute W-1 by replacing large values wii on the diagonal with wii-1 while replacing low values with 0.


## 3.3.4. Damping

It is safer to replace (3.34) with
(3.39)yk+1=yk-qΔyk,
where the damping factor q is found in the following way. We start with q=1; we compute yk+1 and ℱ(yk+1). Then we compare the norms of residues; if |ℱ(yk+1)|<|F(yk)|, we are done. If not, we repeatedly decrease q (by a factor, say, 2) and we recompute yk+1 and ℱ(yk+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 computingW-1, large values wii are replaced with wii-1 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 wii in decreasing order. First we replace the lowest singular value with zero, we compute Δyk and then yk+1 and |ℱ(yk+1)|. Then we replace the two lowest singular values and compute |ℱ(yk+1)| again. This is done for all possible choices and the case giving the lowest residue |ℱ(yk+1)| is used.


## 3.3.6. How Many Iterations

The other question is how many iterations (3.34) should be done. Initially the residue |ℱ(yk+1)| decreases with k 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 massmB suspended on a homogeneous elastic spring of mass mS with the elasticity constant k; see Figure 1. Assume a homogeneous gravitational field (0,-g).Figure 1
Heavy spring elastic pendulum consists of a point massmB attached to a spring with mass mS which is fixed at the other end point.In [1] we derive the equations of motion
(4.1)ẍ1=(1x12+x22-1)x1,ẍ2=(1x12+x22-1)x2-p.
This system has the only parameter
(4.2)p=(mB+mS/2)gkℓ0,
where ℓ0 is the length of the unloaded spring.In [1] we also derive the stability condition for the vertical oscillations
(4.3)x1=0,x2=-1-p+asint,
where a is the amplitude of the vertical oscillation given by initial conditions. The vertical oscillation of the heavy spring elastic pendulum with the relative amplitude a is stable (see [1]) for
(4.4)|a|<|aC(p)|,
where
(4.5)aC(p)≐(3p)2/3-1.
The error of this estimate of aC(p) is less than 0.005 for 0<p<0.6.Figure2 shows the regions in the p-a plane with stable (light gray, or blue in color version) and unstable (dark gray, or red in color version) vertical oscillation with amplitude a and the parameter p. For p and a sufficiently small, the vertical oscillations are stable. For a fixed a=0.1 let us observe what happens when we increase p. When p crosses the critical value p1≐(1/3)(1-a)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 p for the critical value p2≐(1/3)(1+a)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 amplitudea and the external field p.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(4.6)E(x1,x2,ẋ1,ẋ2)=Epot(x1,x2)+Ekin(ẋ1,ẋ2),
where
(4.7)Epot=12(x12+x22-1)2+px2,Ekin=12(ẋ12+ẋ22).
We are interested in periodic solutions with the total energy
(4.8)E(x1,x2,ẋ1,ẋ2)=E0,
where E0 is the energy of the vertical oscillations (4.3):
(4.9)E0=E(0,-1-p+asint,0,acost)=12(a2-p2-2p).To find the initial approximation of one periodic solution, we choose a parameter valuep<p2 and a=0.1 Then we can use numerical integration and Poincaré section given by the horizontal plane x2=-1-p 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 x1 and ẋ1. From the position of the Poincaré section we know that x2=-1-p and we compute ẋ2 from the known total energy. We estimate the period to be T=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) v0=(0,0,0,0,-1,0) meaning that we want to start with decreasing parameter p. 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 amplitudea=0.1. When decreasing the parameter p, the vertical oscillation becomes unstable for the critical value p2≐0.3846. Period doubling bifurcation gives rise to a one parameter family of periodic solutions that exist down to the parameter value p3≐0.0655. We use the above algorithm for continuation of this branch of periodic solutions.Figure4 shows how the period T of the periodic solution grows with decreasing parameter p from the value T=4π just after the period doubling bifurcation up to the value T=8π.Figure 4
The periodT of periodic oscillation (as a function of the parameter p) changes from 4π for high p up to 8π for low p.Figure5 shows for the same parameter values how the range of x1 coordinates (the horizontal displacement of the bob) grows and then decreases down to zero.Figure 5
The rangeΔx1=x1max-x1min of periodic oscillation as a function of the parameter p.When decreasing the parameterp the shape of the trajectory in the x1-x2 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 x1. As the parameter p goes down, the range of x1 grows but the cup-like shape of the oscillation remains. For a still decreasing parameter p the top parts of the trajectory begin to bend and then they turn down. The period T grows. For parameter values close to p3≐0.0655 the trajectory is four times folded and the range of x1 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.The projection of the trajectory to thex1-x2 plane for four values of the decreasing parameter p. Top left: For p, slightly below the critical value p2≐0.3846 the period doubling bifurcation gives rise to the new oscillation with twice the period and with a small x1 range. Top right: For decreasing p the cup-like shape of the trajectory widens and the top parts start to bend. Bottom left: For still lower p 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 p=p3≐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 p. 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.



---

*Source: 104547-2009-06-02.xml*