# On the application of the Lindstedt-Poincaré method to the Lotka-Volterra system

###### Abstract

We apply the Lindstedt-Poincaré method to the Lotka-Volterra model and discuss alternative implementations of the approach. By means of an efficient systematic algorithm we obtain an unprecedented number of perturbation corrections for the two dynamical variables and the frequency. They enable us to estimate the radius of convergence of the perturbation series for the frequency as a function of the only model parameter. The method is suitable for the treatment of systems with any number of dynamical variables.

## 1 Introduction

In the last three decades there has been some interest in the application of perturbation theory to nonlinear dynamical systems such as the Lotka-Volterra model. Murty et al[1] applied perturbation theory to a three-species ecological system and obtained the first perturbation correction to the population of each species. However, they did not take into account that the secular terms spoil the approximate result that will not exhibit the expected periodic behaviour. Grozdanovski and Shepherd[2] applied the well-known Lindstedt-Poincaré method to remove secular terms and obtained the first two perturbation corrections to a two-species system. Consequently, their approximate results exhibit the expected periodic behaviour. Navarro[3] also applied the Lindstedt-Poincaré method to the same two-species system and obtained periodic approximate expressions of second order. This author proposed a symbolic algorithm for the computation of periodic orbits but surprisingly did not show results of larger order. Navarro and Poveda[4] applied Navarro’s approach to a three-species system and derived perturbation corrections of first and second order for the populations of the three species for some particular values of the model parameters. All these studies lead to the conclusion that the Lindstedt-Poincaré method gives reasonable results for some model parameters and suggest that the technique may be useful for the analysis of more realistic and related nonlinear dynamical problems.

Unfortunately Grozdanovski and Shepherd[2] did not explicitly indicate the initial conditions chosen for the solution of the first-order differential equations that provide the corrections at every perturbation order. Since their strategy is not clearly delineated it is difficult to derive a systematic approach for the calculation of perturbation corrections of greater order. On the other hand, Navarro[3] and Navarro and Poveda[4] put forward a systematic symbolic algorithm but they did not appear to exploit it to obtain perturbation corrections of large order. Besides, their presentation of the algorithm appears to be rather obscure for anybody who is not familiar with such technique.

The aim of this paper is the analysis of the approach proposed by Grozdanovski and Shepherd[2] with the purpose of deriving a systematic method for the calculation of perturbation corrections of any order to the Lotka-Volterra model. Such results may give us some clue about the convergence properties of the perturbation series. In addition to it we investigate the possibility of generalizing the method for the treatment of more realistic systems with more than two degrees of freedom.

In section 2 we outline the Lotka-Volterra model and compare alternative implementations of the Lindstedt-Poincaré approach. In section 3 we put forward a generalization of the method that enables one to treat all the previously discussed cases. In section 4 we carry out a large order calculation of the perturbation corrections and estimate the radius of convergence of the perturbation series for the frequency. Finally, in section 5 we outline the application of the method to more general and realistic dynamical systems and draw conclusions.

## 2 The Lindstedt-Poincaré method

In this section we briefly discuss the Lotka-Volterra model and delineate the application of the Lindsted-Poincaré technique. For concreteness we will follow Grozdanovski and Shepherd[2] because their treatment of the dynamical equations is clear and straightforward.

### 2.1 The model

The dynamical equations for the model are

(1) |

where , , and are positive parameters and the point indicates derivative with respect to time. The physical meaning of the parameters is not relevant for present purposes because we are mainly interested in the success of the perturbation approach. Besides, most probably nobody will apply the Lotka-Volterra model to an actual ecological system today because it is quite unrealistic. The interested reader may resort to the papers cited above for more information[1, 2, 3, 4] (and the references cited therein).

We can get rid of some model parameters by means of the following transformations of the independent and dependent variables

(2) |

The resulting equations

(3) |

depend on just one parameter .

There is a stationary point at and . Therefore, if we define

(4) |

the resulting dynamical equations will depend on the perturbation parameter

(5) |

In order to apply the Lindstedt-Poincaré method we define the dimensionless time

(6) |

where is the unknown frequency of oscillation. In this way we have

(7) |

Following Grozdanovski and Shepherd[2] we are using the same symbols and for the solutions of equations (5) and (7). Besides, we have chosen a dot to indicate the derivative with respect to either or . Although this practice may be unwise when one is studying a practical problem and wants to reconstruct and from and it is harmless in the present case because our aim is to show how to obtain perturbation corrections of large order and study the convergence of the perturbation series.

### 2.2 Perturbation equations

### 2.3 Systematic approach

The purpose of this subsection is to derive general expressions for the solutions to the perturbation equations (10) that enable us to develop a systematic algorithm for the calculation of corrections of sufficiently large order.

Neither Grozdanovski and Shepherd[2] nor Navarro[3] consider the initial conditions of the perturbation corrections and explicitly. Here we choose

(11) |

because they greatly facilitate the calculation of and from and :

(12) |

Note that, given and we obtain the product and . Later on we will show why always appears associated to the perturbation parameter in this particular way.

In order to solve equations (10) we rewrite them in matrix form

(13) |

so that the solution is simply given by

(14) |

where

(15) |

Note that equation (14) is consistent with the initial conditions (11).

In order to identify the resonant terms that would give rise to secular terms we rewrite the first-order differential equations as second order ones; for example

(16) |

Therefore, we set so that

(17) |

These equations are a generalization of the Lemma 1 in the paper by Grozdanovski and Shepherd[2] and the proposal of Navarro[3]. It is worth noting that the same value of satisfies both equations (17). We are not aware of a rigorous proof of this result but we can test it by means of our calculations of large order. This point was not discussed in the earlier papers on the application of the Lindstedt-Poincaré method to multidimensional systems[2, 3, 4] probably because they did not try to obtain a set of explicit equations for a systematic application of the approach.

Unfortunately, the perturbation corrections obtained in this way are considerably more complicated than those derived by Grozdanovski and Shepherd[2]. For example, at first order we obtain

(18) | |||||

The coefficients of and agree with the ones derived earlier by those authors and the remaining terms are necessary to satisfy the initial conditions (11). We also find that removes the resonant terms. The perturbation corrections derived by Navarro[3] with somewhat different initial conditions appear to be simpler but they are restricted to .

The perturbation corrections of second order are so complicated that we do not show them here. Besides, is nonzero and a rather cumbersome function of and :

(19) | |||||

The occurrence of rather too complicated perturbation corrections appears to be the price that one has to pay for obtaining the simpler equations (12) for the calculation of and .

At first sight the dependence of on the phase may appear to be the consequence of a wrong calculation. However, we have verified that present solutions already satisfy the perturbation equations and comparison with numerical results reveals a good agreement. For example, figure 1 compares the curve for , and calculated by perturbation theory of zeroth and second order and an accurate numerical result. We appreciate that the addition of the perturbation corrections shown above already improves the analytical results. In the next subsection we will explain the reason for the discrepancy between our expressions and those of Grozdanovski and Shepherd[2] in a more transparent way.

### 2.4 The straightforward Fourier expansion

The purpose of this subsection is merely to show why it is possible to obtain many different solutions at every order of perturbation theory. We may solve the differential perturbation equations (10) by inserting Fourier expansions of the form

(20) |

For the first order we obtain (we omit the superscript for simplicity)

(21) |

We appreciate that if we choose we obtain exactly the results of Grozdanovski and Shepherd[2]. However, there is an infinite number of perfectly valid solutions that emerge from arbitrary choices of , , and provided that they satisfy the ratios . One of such possible solutions is that shown above that satisfies the boundary conditions (11). The solutions derived by the authors just mentioned seem to be the simplest ones and are therefore most convenient for large-order calculations. We should find suitable general conditions to produce such simple results at every order of perturbation theory.

## 3 Generalization of the systematic approach

In the preceding sections we discussed two possible solutions: those that lead to the simple initial conditions (12) and those that are considerably simpler but lead to somewhat complicated initial conditions. The problem at hand is that we have not yet specified the initial conditions for the perturbation equations that lead to the latter. Simpler solutions are obviously most convenient for the calculation of analytic perturbation corrections of large order because they will render the computation algorithm more efficient and less time and memory consuming.

Fortunately, it is not difficult to make the general approach of subsection 2.3 more flexible so that it yields results that are as simple as those of Grozdanovski and Shepherd[2]. We simply choose general initial conditions of the form

(22) |

Now the solution to the matrix perturbation equations (13) is given by

(23) |

and we can choose the arbitrary real numbers and so that the pair of solutions at order is as simple as possible. In what follows we simply set them so that the coefficients of and in vanish (we can, of course, choose instead). More precisely, and are solutions to the equations

(24) |

It is obvious that in this way the solutions and are completely determined.

To second order we have

(26) | |||||

and

(27) | |||||

These solutions are different from those of Grozdanovski and Shepherd[2] but all of them satisfy the perturbation equations. We would have obtained exactly their results if we had chosen and that make the coefficients of and in vanish. We just did it in this way to stress the ambiguity of the results already outlined above in subsection 2.4. Note that if we substitute for in equations (24) we modify the, in principle arbitrary, initial conditions for the solutions to the perturbation equations (10). In either case we have .

The perturbation corrections of third order are given by the coefficients (we again omit the superscripts)

(28) |

From all these results we obtain

(29) |

that was not calculated by earlier authors as far as we know.

Assisted by available computer algebra software we have calculated , , , , , , , interactively and our analytical results suggest that , for the boundary conditions (22) given by (24). Here we just show the next two perturbation corrections to the frequency:

We want to point out that up to this point we have carried out the calculation order by order interactively (that is to say: without programming the equations for the calculation of the perturbation corrections). Obviously, this strategy is unsuitable for the calculation of large order we are interested in. However, even in this rather inefficient way we derived perturbation corrections of order larger than those shown by Navarro[3] who proposed a symbolic algorithm for this purpose.

In closing this section we want to make a couple of considerations about the perturbative solution of this model. To begin with note that we can rewrite equation (7) as

(31) |

Therefore, we can obtain and from the solutions to the perturbation equations for and perturbation parameter . This transformation is convenient because we will not have in the analytic solutions to the perturbation corrections which results in the use of less computer memory. Note that Grozdanovski and Shepherd[2] also defined the parameter to write their expressions for and in a more compact way. However, they did not appear to exploit this fact in a systematic way.

## 4 Large-order calculations

In this section we will show that the algorithm discussed in section 3 is actually useful for the calculation of perturbation corrections of sufficiently large order and will exploit the fact that the perturbation equations (13) and (23), supplemented by (17) and (24), are suitable for programming in available computer algebra systems. For concreteness and simplicity we will focus on the perturbation series for the frequency

(32) |

and will try to determine its radius of convergence.

In general, the radius of convergence of the power-series expansion

(33) |

is determined by the singularity of the function closest to the origin: . There are many ways of estimating the singularities of an unknown function from its known power-series expansion. One of them is given by the Padé approximants[5]

(34) |

where the coefficients and are chosen so that

(35) |

It is commonly assumed that the stable zero of (as and increases) closest to the origin provides an estimate of .

In some cases it is more convenient to resort to quadratic Hermite-Padé approximants[5]

(36) |

where the coefficients , and are chosen so that

(37) |

In this case the singularity closest to the origin is a stable root of

(38) |

With the algorithm of section 3 we have been able to obtain for as analytical functions of . By means of diagonal Padé () and Hermite-Padé () approximants we estimated for the expansion variable . Figure 2 shows a good agreement between both types of approximants. We appreciate that the radius of convergence is a monotonously decreasing function of the model parameter . In a recent paper Amore et al[6] calculated the radius of convergence of the frequency of the van der Pol oscillator with unprecedented accuracy by means of Hermite-Padé approximants constructed from the Lindstedt-Poincaré series with an extremely large number of terms. We are therefore confident of the accuracy of present results.

As an additional verification of the accuracy of our results we have carried out a perturbation calculation of order for and obtained and with Hermite-Padé approximants and , respectively. On the other hand, the diagonal Padé approximants exhibit a stable pole close to the origin at . Based on these results we can safely conclude that .

As mentioned before Navarro[3] developed a symbolic algorithm for the computation of the Lindstedt-Poincaré perturbation corrections and applied it to the Lotka-Volterra model but did not show any results beyond second order. As far as we know there is no perturbation calculation in the literature of order as high as the one shown here. The usefulness of such calculation is obvious because it enables us to estimate the practical range of utility of the Lindstedt-Poincaré perturbation theory for the treatment of dynamical systems. In the case of the Lotka-Volterra model we clearly appreciate that this approximation is not valid if the initial populations and are such that .

## 5 Further comments and conclusions

In this paper we have developed a systematic method for the application of the Linstedt-Poincaré perturbation theory to the Lotka-Volterra model. In particular we discussed the initial conditions for the perturbation equations that were not taken into account explicitly in earlier papers[2, 3]. Present analysis reveals that one can obtain an infinite number of solutions to the perturbation equations and the choice of one of them depends solely on convenience. Here we weighted the possibility of simpler expressions for the calculation of the parameters and on the one side against the simplicity of the solutions on the other. In the latter case we obtained perturbation corrections of considerably larger order than those derived earlier[2, 3]. From them we could estimate the radius of convergence of the perturbation series for the frequency. This result is important for the estimation of the range of validity of the approximate perturbation solutions to the dynamical equations. It shows that the resulting analytical expressions are bounded to fail for some initial conditions.

Present approach can be easily generalized to periodic nonlinear systems of any number of dynamical variables. For example, if we can rewrite the perturbation corrections to the dynamical equations in the form

(39) |

where is an matrix and and are dimensional column vectors, then the solution to the perturbation equation of order is given by

(40) |

where is an dimensional column vector with arbitrary elements that we choose in order to obtain the simplest solutions. In order to carry out this calculation we just need but its construction is a textbook excercise[7].

## Acknowledgments

The research of P. Amore was supported by the Sistema Nacional de Investigadores (México)

## References

- [1] K. N. Murty, M. A. S. Srinivas, and K. R. Prasad, “Approximate analytical solutions to the three-species ecological system”, J. Math. Anal. Appl. 145, 89-99 (1990).
- [2] T. Grozdanovski and J. J. Shepherd, “Approximating the periodic solutions of the Lotka-Volterra system”, ANZIAM J. 49, C243-C257 (2008).
- [3] J. F. Navarro, “A symbolic algorithm for the computation of periodic orbits in non-linear differential systems”, J. Adv. Appl. Math. 1, 160-174 (2016).
- [4] J. F. Navarro and R. Poveda, “Computation of periodic orbits in three-dimensional Lotka-Volterra systems”, Math. Meth. Appl. Sci. 40, 1-16 (2017).
- [5] R. E. Shafer, “On quadratic approximation”, SIAM J. Numer. Anal. 11, 447-460 (1974).
- [6] P. Amore, J. P. Boyd, and F. M. Fernández, “High order analysis of the limit cycle of the van der Pol oscillator”, arXiv:1711.09978 [math.NA].
- [7] T. M. Apostol, Calculus, Second ed. (Blaisdell, Waltham, 1969).