J. Symbolic Computation (1995) 11, 1–000

Normal Form Method and Seminumerical Solutions of Nonlinear ODEs ^{†}

Victor F. Edneral

Institute for Nuclear Physics

of Moscow University,

Leninskie Gory,

Moscow, 119899, Russia

e-mail: edneral@thp.npi.msu.su

(Received )

1. Introduction

A normal form method is based on a transformation of ODEs system to a simpler set called the normal form. The importance of this method for an investigation of ODEs near a stationary point has been recognized for a long time. For the history of this subject see, for instance, (Arnold, Anosov, 1988) or (Guckenheimer, Holmes, 1986). Definitions of a normal form and a normalizing transformation can be formulated in different approaches for general and for special cases. So there are very developed approaches for Hamiltonian systems, see, for example, (Deprit, 1969), (Hori, 1966), (Bruno, 1972), (Bruno, 1990, chapters 1,2). There are also many algorithms (and their realizations) for creating of normal forms and corresponding transformations. For the Hamiltonian case see an improved algorithm of Deprit and Hori in (Mersman, 1970) and its symbolic algebra realization under REDUCE system in (Shevchenko, Sokolsky, 1993), a matter of a numerical creation of normal forms of Hamiltonians is described in (Godziewski, Maciejewski, 1990). Questions of integrability and convergence of the normalizing transformation in Hamiltonian case are discussed in (Bruno, 1972, 1990), (Ito, 1989, 1992). Concerning algorithms for the creation of a normal form in a general case we mention here (in an addition to Bruno’s books) papers (Walcher 91,93, Vallier, 1993). A connection between symmetry properties and convergence of a normalizing transformation have been investigated in (Bruno, Walcher, to be printed).

In this paper we will use the algorithm based on the approach which was developed by A.D.Bruno (Bruno, 1971, 1972, 1989) for Poincare–Dulac’s normal form. The important advantage of this approach is a possibility to investigate a wide class of equations in the united, easy algorithmized frame. In particular it provides a constructive way for a creation of approximations of families of periodic and conditionally periodic solutions in a form of power/Fourier series for real families or in a form of power series in time dependent exponents for complex families. For this paper it is especially important that a problem of analyticity of used transformations is investigated and we are guaranteed

† Supported by CNPq Brazil, project #300894/93-7.

0747–7171/90/000000 + 00 $03.00/0 \begin{matrix}c \\ ◯ \end{matrix} 1995 Academic Press Limited

that a domain of a convergency of corresponding (infinite as a rule) series is a nonzero one. The latest circumstance allows us to calculate periodic solutions and corresponding frequencies near stationary points with acceptable precision by finite formulaes and to be sure that with guaranty we can increase this precision by a calculation of highest terms. Except solutions themselves we can get origin conditions, which initiate periodic solutions in the form of series in parameters or/and (for Hamiltonian systems) in the terms of the mechanical energy. I.e. we can produce some elements of a phase analyses of the systems.

Another advantage of the using approach is an algorithmic simplicity of a creation of the normal form and the corresponding transformation. We have for this procedure a direct recurrent formulae. Its usage does not demand keeping of some large intermediate results as it is in other algorithms. The approach is free from a necessity to solve intermediate systems of linear equations and from restrictions on low resonance cases.

Below we describe a creation of the normal form and its application to a creation of periodic solutions of well known equations. We assist our consideration with concrete examples produced by the NORT program (Edneral, Khrustalev, 1985, 1992) or by NORT and REDUCE system. The intention of the author is a demonstration of the method on so transparent as it is possible examples. Specularity of this paper is also in an attempt to use a common approach to Hamiltonian and non Hamiltonian systems.

2. Applications of normal forms

Let us consider a system of ordinary autonomic nonlinear differential equations:

\frac{d\overline{X}}{dt} = \overline{Φ}(\overline{X}, \overline{a})  (2.1)

here \overline{Φ} = \{Φ_{1}, …, Φ_{n}\} - a vector of right sides of the system which are functions in \overline{X} = \{X_{1}, …, X_{n}\} — variables of the system and in \overline{a} = \{a_{1}, …, a_{m}\} — its parameters. It means that \overline{X} are unknown functions in \overline{a} and time — t.

Such a type of equations originates from many scientific and engineering problems where oscillations, vibrations or wave processes take place. The standard way of an investigation of such systems is:

1 a bifurcation analysis, i.e. an investigation of a picture of a behavior of system solutions in dependence on parameters. It is especially important if sharp change of this behavior occurs at some their values;

2 a phase picture, i.e. an investigation of a behavior of system solutions in dependence

on initial conditions; 3 a calculation of system solutions.

If we have analytic solutions of a system, we have a clear picture of behavior of the system, but we can have such solutions very rarely. As a rule we have numerical solutions, but a numerical investigation of the above items 1 and 2 is sometimes a complicate problem. It is not simple also to create numerical solutions in unstable cases. As a sequence we may be interested in a some ”intermediate” approach which would be between analytic and numerical methods.

The normal form method widely used for a bifurcation analysis. About methods of such an investigation see, for example, (Hassard, Kazarinoff, Wan, 1981) and (Guckenheimer, Holmes, 1986). At these books you can see that a numerical bifurcation analyses is based really on the normal form method. The main idea here is in replacing of system (2.1) by some ”model” system with polynomial right sides of a finite order and in recasting then low orders of them to the canonical (normal) form. From the lowest coefficients of the normal form we can make conclusions about a qualitative behavior of the origin system. It is enough to know only lowest orders of the normal form for such an analysis. Sometimes this job is produced by hands, sometimes by computer algebra systems, see, for example, (Rand, Armbruster, 1987).

Here we try to demonstrate, that calculation of high orders of the normal form can be useful also for targets 2 and 3 mentioned above. We restrict our consideration here mainly by a creation of approximations to periodic solutions but we have a positive material of using this approach in a case of conditionally periodic solution.

3. The normal form method

At the investigation of equations of type (2.1) there three preliminary steps are usually produced. The first step is a shift of variables of the system \overline{X} in such a way that we have in new variables: \overline{Φ}(0, \overline{a}) = 0. The constant zero vector will be a solution of the system. A corresponding point (zero here) is called an ”immovable”, ”stationary” or ”equilibrium” point. A system can have more then one such a point. Each of them is to be investigated separately. The second step is a reduction of the system to the ”model” form which approximates vector \overline{Φ}(\overline{X}, \overline{a}) by a set of polynomials of a finite order. If \overline{Φ} is an analytical function of \overline{X} in some domain of the stationary point it can be realised with any desirable precision. Often this step is made simultaneously with a reduction of the system to the central manifold. After that right sides of the system will be polynomials without constant (absolute) terms. By a linear complex change of its variables it is possible to reduce a linear parts of the system to the canonical form. After this (the third) step we have:

\frac{d\overline{Y}_{i}}{dt} = Â\overline{Y} + \overline{F} (\overline{Y} , \overline{a})  (3.1)

where Â matrix has a canonical (Jordan’s) form and \overline{Y} – new variables of the system after the above transformations.

In this paper we suppose also that matrix Â has a diagonal form. Really the described method is formulated without such a restriction. It has rather a technical character and can be eliminated with a development of the NORT program. We suppose also that the dependence of the right side in parameters has a polynomial character also. It has a more principal character.

Finally we will have the system which is considered in this paper (below for a simplification of a notation we will hint the dependence on parameters):

\frac{dY_{i}}{dt} = λ_{i}Y_{i} + F_{i}(Y_{1}, Y_{2}, …, Y_{n}), i = 1, …, n  (3.2)

here {λ_{i}} is a non-zero vector of eigenvalues of the linear parts of the right sides. {F_{i}} is a vector of polynomials of finite degrees without constant and linear terms. In other words we impose further restrictions:

• we take into account autonomic † systems with a polynomial nonlinearity

• systems are investigated near a stationary point, i.e. the vector {(0, …, 0)} is a solution of the system;

• we suppose that a linear part of the right hand side can be diagonalized and this operation has been done (actually the algorithm can be generalized for a case of a triangular matrix (Bruno, 1988). But anyway the eigenvalues must be know). We also suppose that its eigenvalues λ_{i} have numerical values†.

Let us now rewrite the equation (3.2) in the form:

\frac{dY_{i}}{dt} = λ_{i}Y_{i} + Y_{i} ∑_{\overline{Q}_{i}} F_{i,\overline{Q}_{i}} \overline{Y}^{\overline{Q}^{i}},  (3.3)

where \overline{Y}^{\overline{Q}^{i}} is defined as a product:

∏_{j=1}^{n} Y_{j}^{Q_{ij}}, Q_{ij} = \{ \begin{matrix}0, 1, 2, … & if i ≠ j, \\ −1, 0, 1, … & if i = j. \end{matrix}

Q_{ij} here are components of a vector \overline{Q}_{i} = \{Q_{ij}\}. Q_{ij} are non-negative integers with the exception of {Q_{ii}} ≥ −1, because factor Y_{i} has been carried out from sum in (3.3). Summation in (3.3) is over all possible \overline{Q}_{i}-vectors.

In such terms a normalizing transformation may be written as:

Y_{i} = Z_{i} + Z_{i} ∑_{\overline{Q}_{i}} H_{i,\overline{Q}_{i}} \overline{Z}^{\overline{Q}^{i}}.  (3.4)

Coefficients H in transformation (3.4) can be chosen so that equation (3.3) will have the form:

\frac{dZ_{i}}{dt} = λ_{i}Z_{i} + Z_{i} ∑_{(\overline{Q}_{i},\overline{λ})=0} G_{i,\overline{Q}_{i}} \overline{Z}^{\overline{Q}^{i}}.  (3.5)

The great difference between (3.5) and (3.3) is a restriction on the domain of summation, which is defined by the equation:

(\overline{Q}_{i}, \overline{λ}) =^{def} ∑_{j=1}^{n} Q_{ij}λ_{j} = 0.  (3.6)

The form (3.5) is called the normal form for (3.3). The H and G coefficients in (3.4) and (3.5) may be found from the recurrent formulae (Bruno, 1971, 1989):

† Actually periodic dependence on time can be eliminated by changing time dependent functions for the proper solutions of additional ODEs. For instance cos(ωt) can be simulated by a couple of additional autonomic equations with proper origin conditions.

† This is not a very dramatic requirement, see, e.q., a treatment of van der Pol’s equation below.

\begin{matrix}G_{i,\overline{Q}_{i}} + (\overline{Q}_{i}, \overline{λ})H_{i,\overline{Q}_{i}} = \\ − ∑_{j=1}^{n} ∑_{\overline{P} +\overline{R}=\overline{Q}_{i}}(P_{j} + δ_{ij})H_{i,\overline{P}} G_{j,\overline{R}} \\ +Φ_{i,\overline{Q}_{i}}. \end{matrix}  (3.7)

Here the second summation in the right-hand side is over all vectors \overline{P} and \overline{R} for which \overline{P} +\overline{R} = \overline{Q}_{i}; P_{j} in parentheses is the j- component of vector \overline{P}; Φ_{i,\overline{Q}_{i}} is a coefficient at the

factor Z_{i}\overline{Z}^{\overline{Q}^{i}} in the polynomial F_{i} from (3.3), arguments of which have been transformed by (3.4) †.

The uncertainty in (3.7) is fixed below by conventions:

\begin{matrix}H_{i,\overline{Q}_{i}} = 0  if (\overline{Q}_{i}, \overline{λ}) = 0, \\ G_{i,\overline{Q}_{i}} = 0  if (\overline{Q}_{i}, \overline{λ}) ≠ 0 \end{matrix}  (3.8)

Sometimes the normalizing transformation fixed by (3.8) is called a ”basic transformation”.

3.1. the main algorithm

The algorithm of the calculation of G and H is based on (3.7),(3.8). Let us choose the representation of sets of coefficients G_{i,\overline{Q}_{i}} and H_{i,\overline{Q}_{i}} so that they would be collected in such groups where a group of order n contains only elements with such vector-indexes \overline{Q_{i}} = \{q_{i,j}\} that the sum of their components is equal to the order, i.e. ∑_{j} q_{i,j} = n for each i. You can calculate the group of the next order for G and H by using groups of G and H with smaller order only, i.e. (3.7) is a recurrent formulae.

The algorithm:

Let n is a dimension of the system, for its normalization till order m we are to do:

(i) for i = 1, 2, …, n do:

Calculate all square in Y_{i} elements in right sides of F_{i}(\overline{Y}), i.e. calculate the group of the first order of the set Φ_{i,\overline{Q}_{i}} and sort it into two sets in dependence on a value of scalar product (3.6). The first set where this product is zero will be the first order group of G_{i} and the second set after a division on the value of the corresponding scalar product will be the first order group of H_{i}

(ii) for k = 2, 3, …, m do:

(a) for i = 1, 2, …, n do: calculate the group of order k of the set Φ_{i} for which substitute the Y_{i} evaluated by (3.4) till order k − 1 in the F_{i}(\overline{Y} (\overline{Z})) and choose coefficients at monomials

Z_{i}\overline{Z}^{\overline{Q}^{i}} as Φ_{i,\overline{Q}_{i}};

(b) for i = 1, 2, …, n do:

Calculate the groups of G_{i} and H_{i} of order k by a division the set Φ_{i} into two

† Take into account that the scalar production (\overline{Q}_{i}, \overline{λ}) will appeare in denominators of H coefficients with any possible values of \overline{Q}. If one of eigenvalue has a symbolic value this scalar production will provides a series of (complex in general) poles in this symbolic value.

sets as in item (a). After that you will have the full group of order k for G_{i} and preliminary values of such a group for H_{i}, because those values of the vector \overline{Q}_{i} which takes place at non-zero elements into the summation in (3.7) can not zero the scalar product, as a sum of a vector which zeroes this product (\overline{R}) and a vector which does not do it (\overline{P}).

(c) for i = 1, 2, …, n do:

for j = 1, 2, …, n do:

supplement the group of order k of H_{i} with properly sorted multiplications of all elements of such groups of H_{i} and G_{j} so that their summary order (i.e. ∑_{s=1}^{n}(P_{s} +R_{s})) would be equal to k. These multiplications would be multiplied also by a factor (P_{j} + δ_{i,j}), which can be equal to zero and therefore should be calculated before the other factors. Before the supplement all elements are to be divided on the corresponding scalar products.

It is not so heavy to estimate a cost of this algorithm but here such estimation has no a deep sense because of relative high complexity of a calculation of right sides of the nonlinear system. Under such circumstances it is very important to calculate these parts very economic, using as possible the fact that we need at the each step of (ii) a single homogeneous ”slice” of \overline{F} of order k only and its terms of lower orders will be not changing during the following operations. The problems with the optimization of the right sides evaluation is one of the main reason of an absence of closed interface for the NORT package.

3.2. periodic solutions, a convergence of series and adequacy of the normal form to the origin system

The sums in (3.4) and (3.5) may include an infinite number of terms for physically interesting cases even though the sum in (3.3) has a finite number of ones. The question about properties of a convergence and a smoothness of transformation arises. For adequacy of the normal form to the origin system in a some domain of a stationary point we need to have analytical normalizing transformation. The problem of analyticity of (3.4) have been investigated in (Bruno, 1971, 1989). The result is shortly that the transformation (3.4) is divergent in a general case, but there can be found some subsets of the formal invariant sets of the normalized system (3.5) where the normalizing transformation is convergent and which are therefore adequate to some invariant subsets of the original system. Bruno calls them analytical sets.

We will see below the subspace of a coordinates space of normalized system which includes only variables Z_{i} for which Re(λ_{i}) = 0. Actually, if this is a null-space we have no deep problems with the convergence, because in this case there is a result about topological equivalence of the original system to the system which is normalized till a finite order †. The variables for which corresponding eigenvalues are not pure imaginary increase to infinity (for Re(λ_{i}) > 0) or vanish (for Re(λ_{i}) < 0) with time. Therefore we have for periodic solutions only the above mentioned subspace. Let us extract in this subspace the set A:

† We will not discuss here this circle of questions and send a reader for the references to Bruno’s books.

\begin{aligned}A = \{\overline{Z} : & ψ_{i} = λ_{i}Z_{i}a, & if Re(λ_{i}) = 0; \\ & Z_{i} = 0, & if Re(λ_{i}) ≠ 0; i = 1, …, n\} \end{aligned}  (3.9)

where ψ_{i} is an element of a right side of the normalized system and a is a free parameter, the same for all values of i.

Theorem 3.1 The set A contains periodic and conditionally periodic solutions of the normalized system including stationary points.

Further let us collect all pure imaginary eigenvalues in groups with the property: each group includes all λ’s which are co-measured in pairs †. Following Bruno we will speak about coordinate subspace which contains only the coordinates (variables) with numbers corresponding one of these groups as about the ”rational” subspace. Let us designate the intersection of the set A with a union of all such rational subspaces by A^{∼}.

Theorem 3.2 The set A^{∼} are filled with periodic solutions of the normalized system including stationary points and there exists an analytical invertible transformation which transfers the original system to the system which is normalized on the set A^{∼}.

The important sequence from the above theorems is that if we find the set A^{∼} we find the periodical solutions of the normalized system, which contracts to the stationary point. Along these solutions the normalizing transformation will be analytical. Thus the normalized system will adequate to the origin one along these solutions.

It is easy to see that a fulfillment of (3.9) is enough for integration of the system, because of a will be a constant. It provides also the way for building integrals of motion.

A set exists also which contains conditionally periodic solutions such that the normalizing transformation is analytic along them. The definition of such a set see in (Bruno, 1989).

3.3. the general scheme of the ODEs investigation

General scheme of an investigation of a nonlinear ODEs system by the normal form method near each stationary point is:

1 recasting the system to a (model) polynomial form;

2 the linear normalization of the system, i.e. the reduction of a linear part of the right sides of the system to Jordan’s or to a diagonal form and an investigation of the corresponding eigenvalues;

3 the non-linear normalization of the system,—the creation of the normal form and the corresponding normalizing transformation; (sometimes it is called recasting the system in a ”resonance form”).

4 a bifurcation analysis of the system in parameters by the investigation of lowest

non-disappeared orders of the normal form; 5 building the periodic and conditionally periodic solutions which contract to the

stationary point; 6 decreasing an order of the normalized system if it necessary and a repetition of

the present investigation from item 2 near each stationary point of the system decreased.

† I.e. a ratio of any non-zero elements within the group is a rational number.

4. Second order systems

All questions about convergence and integrability of normal forms for all second order systems of type 3.2 near a stationary point have been investigated by (Bruno, 1971,1989) for the case when the stationary point is an elementary especial point, i.e. both eigenvalues of the system are not equal zero simultaneously. The result of this investigation shortly is: the normal forms for equations defined as the above are integrable. For each case there are written integrals of the normalized system. The result regarding the convergence is the following. Let λ_{2} ≠ 0 and λ =^{def} λ_{1}/λ_{2}. If Im(λ) ≠ 0 or if λ > 0 then the transformation is convergent and we can produce an approximate solution of the origin system from known integrals of the normal form by transformation (3.4) with any desirable precision. In other words cases of ”focus” and ”knot” can be treated without any additional demands. At a negative irrational λ the convergence will take place if for all non-zero vectors \overline{Q} with integer elements there exist positive ε, ν such that |(\overline{Q}, \overline{λ})| > ε(|q_{1}| + |q_{2}|)^{−ν}. This condition can be checked before creation of the normal form. The latest case is a one partial case of a saddle point. But at real non-positive rational λ = −m/n ≤ 0 we will have convergence under some additional requirements on the normal form. It is interesting case. It includes the cases of a center and a limit circle. Let us look at a couple of examples.

4.1. duffing’s equation

This is an equation of the second order, which origins from a problem of a mathematical pendulum:

\frac{d^{2}ϕ}{dt^{2}} + ω_{0}^{2} sin(ϕ) = 0  (4.1)

where ϕ is an angle of deviation of the pendulum and which can be approximated by the form (units below have been chosen so that ω_{0} = 1):

\frac{d^{2}ϕ}{dt^{2}} = −ϕ + \frac{1}{6}ϕ^{3}  (4.2)

By considering ϕ = \sqrt{6}X we get Duffing’s equation:

\frac{d^{2}X}{dt^{2}} = −X + X^{3}  (4.3)

This is a Hamiltonian equation with the full energy H:

H = \frac{1}{2}(\frac{dX}{dt})^{2} + \frac{1}{2}X^{2} − \frac{1}{4}X^{4}  (4.4)

By linear complex changing for variables:

\begin{aligned}X & = & Y_{1} + Y_{2} \\ \frac{dX}{dt} & = & i (Y_{1} − Y_{2}) \end{aligned}  (4.5)

we can rewrite it in the diagonalized form:

\begin{aligned}\frac{dY_{1}}{dt} & = & i & Y_{1} − \frac{i}{2} (Y_{1} + Y_{2})^{3} \\ \frac{dY_{2}}{dt} & = & −i & Y_{2} + \frac{i}{2} (Y_{1} + Y_{2})^{3} \end{aligned}  (4.6)

Note that a couple of equations the above mentioned have complex conjugated coefficients at terms with exchange Y_{1} ↔ Y_{2} It will take place always if the origin system has real coefficients.

The vector of eigenvalues of a linear part of the system is \overline{λ} = \{i, −i\}.

In accordance with the definition of the normal form:

\frac{dZ_{i}}{dt} = λ_{i}Z_{i} + Z_{i} ∑ ; \begin{matrix}p_{i} ≥ −1 \\ p_{1}, …, p_{i−1}, p_{i+1}, …, p_{n} ≥ 0 \\ (\overline{λ}, \overline{p}) = 0 \end{matrix} g_{i,p_{1},...,p_{n}} ∏_{k=1}^{n} Z_{k}^{p_{k}}, i = 1, …, n

we will have for the case n = 2 at sums in the right sides of equations only terms, where (\overline{λ}, \overline{p}) = i (p_{1} − p_{2}) = 0. I.e. sums at the right sides of the normal form of Duffing’s equation will include only terms where p_{1} = p_{2}:

\begin{aligned}\frac{dZ_{1}}{dt} & = & i & Z_{1} + Z_{1} (g_{1,1,1}Z_{1}Z_{2} + g_{1,2,2}Z_{1}^{2}Z_{2}^{2} + …) \\ \frac{dZ_{2}}{dt} & = & −i & Z_{2} + Z_{2} (g_{2,1,1}Z_{1}Z_{2} + g_{2,2,2}Z_{1}^{2}Z_{2}^{2} + …) \end{aligned}  (4.7)

The condition (3.9) for the second order equation with eigenvalues λ_{1} = −λ_{2} has the form:

∑_{k=1,...} g_{1,k,k} (Z_{1}Z_{2})^{k} = − ∑_{k=1,...} g_{2,k,k} (Z_{1}Z_{2})^{k}

As for any real (originally) equation of the second order the equalities:

\begin{matrix}g_{1,i,i} = g_{2,i,i}^{∗} \\ h_{1,i,j} = h_{2,j,i}^{∗} \end{matrix}  (4.8)

takes place †, this condition for real equations of the second order has the form:

∑_{k=1,...} Re(g_{1,k,k}) (Z_{1}Z_{2})^{k} = 0  (4.9)

After a calculation of the normal form for Duffing’s equation, it can be seen that (4.9) is fulfilled identically because all g_{1,i,i} and g_{2,i,i} are pure imaginary here. It means that in the case of Duffing’s equation we have the case which is usually called the ”center”, when periodic solutions exist for any initial conditions.

Really, by multipling the first from equations (4.7) by Z_{2} and the second one by Z_{1} we will have after their addition \frac{d(Z_{1}Z_{2})}{dt} = 0 and so (4.7) has a trivial form:

† The asterisk here notes a complex conjugation

\begin{aligned}\frac{dZ_{1}}{dt} & = & & i ω(Z_{1}Z_{2}) Z_{1} \\ \frac{dZ_{2}}{dt} & = & − & i ω(Z_{1}Z_{2}) Z_{2} \end{aligned}

because ω(Z_{1}Z_{2}) =^{def} 1 + g_{1,1,1}/i + g_{1,2,2}/i + … is a real constant. Therefore the family of solutions of Duffing’s equation in the normal form is:

Z_{1}(t) = c_{1}e^{+iω(c_{1}c_{2})t}, Z_{2}(t) = c_{2}e^{−iω(c_{1}c_{2})t}

where c_{1}, c_{2} are voluntary constants.

Now we can create the approximation to origin equation (4.3) by substitution the above found Z_{i} into Y_{i} with (3.4) and after that into X. If we choose complex conjugate values for c_{1}, c_{2} we will have an approximation of the real solution in the form of truncated Fourier series.

The calculation by the NORT package till 15^{th} order in C = \sqrt{c_{1}c_{2}} took 53 seconds on 33MHz 4Mbyte RAM PC/AT-486DS computer in the symbolic mode of REDUCE 3.4.2. Due to equalities (4.8) we do not need to calculate series for both of equations (4.6). So we have 16 terms for each sums in normal form (4.7) and 255 terms for the each normalizing transformation from Z_{i} to Y_{i}. The final expression for a frequency has 16 terms and the expression for X(t) has 141 terms.

We got these series as series in the C variable. It is not convenient and for the final representation we calculated the H as a series in C by substituting found X and \frac{dX}{dt} in (4.4). By inverting this series and substituting C as a series of H into expressions for the ω and X we will have the end result as series in H. For an economy of the place we show here this result till the 5^{th} order only †:

ω = 1 − \frac{3}{4}H − \frac{69}{64}H^{2} − \frac{633}{256}H^{3} − \frac{110421}{16384}H^{4} − \frac{1318329}{65536}H^{5} + …
X = \sqrt{2H}×
[ cos(ωt)(1 + \frac{9}{16}H + \frac{271}{256}H^{2} + \frac{10779}{4096}H^{3} + \frac{243613}{32768}H^{4} + \frac{2963587}{131072}H^{5})
− cos(3ωt)H(\frac{1}{16} + \frac{3}{16}H + \frac{1209}{2048}H^{2} + \frac{127233}{65536}H^{3} + \frac{6907221}{1048576}H^{4})
+ cos(5ωt)H^{2}(\frac{1}{256} + \frac{11}{512}H + \frac{3107}{32768}H^{2} + \frac{25567}{65536}H^{3})
− cos(7ωt)H^{3}(\frac{1}{4096} + \frac{1}{512}H + \frac{5805}{524288}H^{2})
+ cos(9ωt)H^{4}(\frac{1}{65536} + \frac{21}{131072}H)
− cos(11ωt)H^{5}(\frac{1}{1048576}) + …]

(4.10)

The result of this calculation was verified by two ways. The first one was a direct substitution of series (4.10) in the original equation (4.3). After that we had only terms with negligible orders of H. The second way was a comparison of the numerical solutions of Duffing’s equation evaluated by Runge–Kutta method (by NAG’s d02baf procedure with a ’tolerance’ = 10^{−12}) and the values tabulated of the series the above mentioned at different values of H. In view of (4.4) we have H = H(X = ϕ/\sqrt{6}, \frac{dX}{dt} = \frac{dϕ}{dt}/\sqrt{6}). On

† If somebody is interested in the full series in the REDUCE readable or FORTRAN form, he can reach them via the name service by anonymous ftp from euler.mat.ufrgs.br (143.54.2.220) in the directory pub/norm

the other hand H has its maximum physically reasonable value when at a zero velocity a maximum deviation of a pendulum takes place, i.e. when ϕ = π/2 and \frac{dϕ}{dt} = 0. In this case we can choose H ≤ H_{max} ≃ 0.163. Let us now introduce a function of a maximum relative error during one period f_{err}(H):

f_{err} = sup_{t∈[0,2π/ω]} \sqrt{\frac{(X_{ser} − X_{num})^{2} + (dX_{ser}/dt − dX_{num}/dt)^{2}}{X_{ser}^{2} + (dX_{ser}/dt)^{2}}}  (4.11)

This function indicates a maximum relative deviation in phase space between series (4.10) X_{ser} and its symbolically evaluated derivation on the one hand and numerical solutions of (4.3) X_{num}, dX_{num}/dt on the another hand as a function of full energy H. We have:

\begin{aligned}f_{err}(H = 0.1) & ≃ 1.8 × 10^{−8} \\ f_{err}(H = 0.125) & ≃ 7.4 × 10^{−7} \\ f_{err}(H_{max} = 0.163) & ≃ 7.4 × 10^{−5} \end{aligned}

You can see that riched precision can be used for practical goals in a physical diapason of energy. A use of a ready finite formulae can be sometimes more preferable at a real time calculations than a numerical solution of differential equations. Of course form (4.10) must be improved for a real usage by standard methods of prepearing of series for numerical tabulation.

4.2. van der pol’s equation

This equation origins from a problem of vibrations in electronic circuits:

\frac{d^{2}X}{dt^{2}} = −X + (ε^{2} − X^{2}) \frac{dX}{dt}  (4.12)

By linear complex changing variables (4.5) it can be rewritten in the diagonalized form:

\frac{dY_{1}}{dt} = i Y_{1} + \frac{1}{2} (Y_{1} − Y_{2}) [ε^{2} − (Y_{1} + Y_{2})^{2}]
\frac{dY_{2}}{dt} = −i Y_{2} + \frac{1}{2} (Y_{2} − Y_{1}) [ε^{2} − (Y_{1} + Y_{2})^{2}]
\frac{dε}{dt} = 0

As for Duffings‘ case these equations have complex conjugated coefficients at terms with exchange Y_{1} ↔ Y_{2}.

We consider here the parameter ε as a new variable. Such trick allows us to free eigenvalues from a parameter dependence. It is necessary to be careful under this operation. We can produce the trick only if the expansion of the expression 1/(\overline{Q}_{i}, \overline{λ}(ε)) is correct at the point ε = 0 for all \overline{Q}_{i} used in the summation in (3.4). (See the footnote on the correspondent page).

And as for Duffing’s equation sums at the right sides of the normal form will include only terms where p_{1} = p_{2}. The third ”additional” equation will be not changed. So a difference between (4.7) and the system presented below is only in polynomial dependence of the right sides on ε:

\frac{dZ_{1}}{dt} = i Z_{1} + Z_{1} ∑_{k=0,1,...} (g_{1,1,1,2k}(Z_{1}Z_{2})ε^{2k} + g_{1,2,2,2k} (Z_{1}Z_{2})^{2}ε^{2k} + …)
\frac{dZ_{2}}{dt} = −i Z_{2} + Z_{2} ∑_{k=0,1,...} (g_{2,1,1,2k}(Z_{1}Z_{2})ε^{2k} + g_{2,2,2,2k}(ε)(Z_{1}Z_{2})^{2}ε^{2k} + …)
\frac{dε}{dt} = 0

(4.13)

The convergent condition (4.9) will have now the form:

Re[ ∑_{j,k=0,1,...} g_{1,j,j,2k} ε^{2k} (Z_{1}Z_{2})^{j}] = 0,  g_{1,0,0,0} =^{def} 0  (4.14)

Contrary to Duffing’s equation the above is not satisfied automatically, but it may be rewrite in the form:

aU + bV + cU^{2} + dUV + eV^{2} + … = 0

with U = Z_{1}Z_{2}, V = ε^{2} and a ≠ 0. This is an experimental fact and it allows us to get the unique (because of the implicit function theorem) solution of (4.14) in the form:

Z_{1}Z_{2} = ∑_{k=1,2,...} q_{k}ε^{2k}

It is easy to see that if the above is satisfied then Z_{1}Z_{2} is a constant in time. So we can continue the evaluations as for Duffing’s case, but now constancys c_{1}, c_{2} are not free (this is a case of a”limit circle”):

\begin{matrix}Z_{1}(t) = c_{1}e^{+iω(c_{1}c_{2})t}, Z_{2}(t) = c_{2}e^{−iω(c_{1}c_{2})t} \\ c_{1}c_{2} = ∑_{k=1,2,...} q_{k}ε^{2k} \end{matrix}

After we create the approximation to the origin equation (4.12) by substitution of the found above Z_{i} into Y_{i} with the help of (3.4) and next in X. If we choose complex conjugate values for c_{1}, c_{2} we also for Duffing’s equation will have an approximation of the real solution in the form of truncated Fourier series.

But form (4.12) is not well for a graphic representation because † a vibration amplitude X is proportional to the value of ε and in a picture we will have different scales for different values of the parameter. We eliminate this problem by introducing for the end result a new variable χ = X/ε. The equation for this variable will have the form:

\frac{d^{2}χ}{dt^{2}} = −χ + ε^{2} (1 − χ^{2}) \frac{dχ}{dt}  (4.15)

We did not choose this form for the investigation because in form (4.12) we choose variables X ∼ ε as having the same weight near zero. Actually we could use form (4.15) but the result would have smaller precision at the same order of the corresponding normal form. For (4.15) we will have (we display it here only till the 16^{th} order):

† Really, you can see from the above that Z_{1}Z_{2} ∼ ε^{2} and X ∼ Y_{i} ∼ Z_{i}.

ω = 1 − \frac{1}{16}ε^{4} + \frac{17}{3072}ε^{8} + \frac{35}{884736}ε^{12} − \frac{678899}{5096079360}ε^{16} + …
χ = cos(ωt)(2 + \frac{1}{64}ε^{4} − \frac{23}{49152}ε^{8} − \frac{51619}{169869312}ε^{12} + \frac{948555443}{19568944742400}ε^{16})

+ cos(3ωt)ε^{4}(−\frac{3}{32} + \frac{101}{12288}ε^{4} + \frac{24061}{28311552}ε^{8} − \frac{279818087}{815372697600}ε^{12}) + cos(5ωt)ε^{4}(−\frac{5}{96} + \frac{1865}{110592}ε^{4} − \frac{328835}{254803968}ε^{8} − \frac{111998015}{293534171136}ε^{12}) + cos(7ωt)ε^{8}(\frac{1379}{110592} − \frac{10923199}{3185049600}ε^{4} + \frac{21049213549}{183458856960000}ε^{8}) + cos(9ωt)ε^{8}(\frac{61}{20480} − \frac{1769369}{589824000}ε^{4} + \frac{161113663733}{237817036800000}ε^{8}) + cos(11ωt)ε^{12}(−\frac{409871}{331776000} + \frac{1359229760383}{1872809164800000}ε^{4}) + cos(13ωt)ε^{12}(−\frac{715247}{3715891200} + \frac{2076538440769}{5243865661440000}ε^{4}) + cos(15ωt)ε^{16}(\frac{526426361}{4661213921280}) + cos(17ωt)ε^{16}(\frac{392636471}{29964946636800}) + sin(3ωt)ε^{2}(−\frac{1}{4} + \frac{15}{512}ε^{4} − \frac{779}{1179648}ε^{8} − \frac{4538017}{6794772480}ε^{12}) + sin(5ωt)ε^{6}(\frac{85}{2304} − \frac{8095}{1327104}ε^{4} − \frac{1252495}{6115295232}ε^{8}) + sin(7ωt)ε^{6}(\frac{7}{576} − \frac{99967}{13271040}ε^{4} + \frac{415949513}{382205952000}ε^{8}) + sin(9ωt)ε^{10}(−\frac{9791}{2457600} + \frac{117258703}{70778880000}ε^{4}) + sin(11ωt)ε^{10}(−\frac{5533}{7372800} + \frac{1657839733}{1486356480000}ε^{4}) + sin(13ωt)ε^{14}(\frac{21731177}{57802752000}) + sin(15ωt)ε^{14}(\frac{138697}{2774532096}) + …

(4.16)

The calculation by the NORT package till 28^{th} order in ε took 21 minutes on the same computer as for Duffing’s case. We had 113 terms for each sums in the normal form (4.13) and 1232 terms for the each normalizing transformation from Z_{i} to Y_{i}. The expression for a frequency has 9 terms and the expression for X(t) has 114 terms. Note, that the power series for a frequency of van der Pol’s equation has been calculated till 164^{th} order in ε in (Andersen, Geer, 1983).

The comparison of this result with a numerical one in the terms (4.11) gives:

\begin{aligned}f_{err}(ε^{2} = 0.5) & ≃ 8.9 × 10^{−9} \\ f_{err}(ε^{2} = 0.75) & ≃ 4.1 × 10^{−6} \\ f_{err}(ε^{2} = 1.0) & ≃ 2.9 × 10^{−4} \end{aligned}

Results of the tabulation of a full (in 28 order) variant of (4.16) for the above values of ε are presented in Fig.4.2.

Besides a solution of equation we can create also the expressions for the origin conditions which provide the periodic solution of (4.15) as a function (series) in ε. Really, we have periodic solution when the expression Z_{1}Z_{2} is a constant in time. Let us inverse series (3.4), (4.5) and change Z_{i} in this production for their expressions via χ, dχ/dt and ε. The result must be constant along the trajectory of the limit cycle too. Therefore we have equation of that trajectory (there is no a clear dependence on time below, but χ and dχ/dt depend on time by implicit way):

\frac{d}{dt}[Z_{1}(χ, dχ/dt, ε)Z_{2}(χ, dχ/dt, ε)] = 0

\frac{d^{2}χ}{dt^{2}} = −χ + ε^{2} (1 − χ^{2}) \frac{dχ}{dt}

-3 -2 -1 0 1 2 3-2.5 -2 -1.5 -1 -0.5 0  0.5  1  1.5  2  2.5^{χ̇}  ^{χ}  ^{ε^{2} = 0.50}^{ε^{2} = 0.75}^{ε^{2} = 1.00}

Figure 1. Solutions of van der Pol’s equation for different values of ε^{2}.

Differenting symbolically the above using the van der Pol’s equation for an evaluation of the second derivation d^{2}χ/dt^{2} we will have the corresponding equation in the form of the series:

∑_{i,j,k=0,1,...} p_{i,j,k} ε^{2k}χ^{i} {dχ \atop dt}^{j} = 0, p_{0,0,0} =^{def} 0

Suppose we are interested in the dependence of the origin condition for the derivation in (4.15) {dχ \atop dt}|_{t=t_{0}} on ε if χ|_{t=t_{0}} = 0. In such a case we substitute zero for χ in the sum above and solve the equation resulted. The result for van der Pol’s equation is:

\begin{matrix}χ = & 0 \\ {dχ \atop dt} = & 2 + {17 \atop 96}ε^{4} − {1577 \atop 552960}ε^{8} − {102956839 \atop 55738368000}ε^{12} + {48722480822161 \atop 157315969843200000}ε^{16} + … \end{matrix}  (4.17)

For the couple χ|_{t=t_{0}}, {dχ \atop dt}|_{t=t_{0}} = 0 we have:

\begin{matrix}χ = & 2 + {1 \atop 96}ε^{4} − {1033 \atop 552960}ε^{8} + {1019689 \atop 55738368000}ε^{12} + {9835512276689 \atop 157315969843200000}ε^{16} + … \\ {dχ \atop dt} = & 0 \end{matrix}  (4.18)

5. Henon–Heiles’ equations as an example of the forth order system

Henon–Heiles’ system originally appeared in the theory of motion of particle in an axial symmetric gravitation field, for example in the problem of motion of stars in the Galactic field (Henon, Heiles, 1964). It is a system of two second order differential equations:

\begin{aligned}\frac{d^{2}X}{dt^{2}} & = & −X − 2XY \\ \frac{d^{2}Y}{dt^{2}} & = & −Y − X^{2} + Y^{2} \end{aligned}  (5.1)

System (5.1) can be rewritten in form (3.2) by the linear complex substitution:

\begin{aligned}X & = & Y_{1} + Y_{2} \\ Y & = & Y_{3} + Y_{4} \\ \frac{dX}{dt} & = & i (Y_{1} − Y_{2}) \\ \frac{dY}{dt} & = & i (Y_{3} − Y_{4}) \end{aligned}

after which we will have instead of (5.1):

\begin{aligned}\frac{dY_{1}}{dt} & = & i & Y_{1} + i (Y_{1} + Y_{2}) (Y_{3} + Y_{4}) \\ \frac{dY_{2}}{dt} & = & −i & Y_{2} − i (Y_{1} + Y_{2}) (Y_{3} + Y_{4}) \\ \frac{dY_{3}}{dt} & = & i & Y_{3} + \frac{i}{2} [(Y_{1} + Y_{2})^{2} − (Y_{3} + Y_{4})^{2}] \\ \frac{dY_{4}}{dt} & = & −i & Y_{4} − \frac{i}{2} [(Y_{1} + Y_{2})^{2} − (Y_{3} + Y_{4})^{2}] \end{aligned}  (5.2)

The non-linear normalization for (5.1) is produced by the transformation:

Y_{i} = Z_{i} + Z_{i} ∑ ; \begin{matrix}p_{i} ≥ −1 \\ p_{1}, …, p_{i−1}, p_{i+1}, …, p_{4} ≥ 0 \end{matrix} h_{i,p_{1},p_{2},p_{3},p_{4}} Z_{1}^{p_{1}}Z_{2}^{p_{2}}Z_{3}^{p_{3}}Z_{4}^{p_{4}}, i = 1, …, 4  (5.3)

The vector of eigenvalues of the linear part will consist of two couples of complex conjugated imaginary units: \overline{λ} = \{i, −i, i, −i\}. I.e. the normal form for (5.2) will be:

\frac{dZ_{i}}{dt} = λ_{i}Z_{i} + Z_{i} ∑ ; \begin{matrix}p_{i} ≥ −1 \\ p_{1}, …, p_{i−1}, p_{i+1}, …, p_{4} ≥ 0 \\ p_{1} + p_{3} = p_{2} + p_{4} \end{matrix} g_{i,p_{1},p_{2},p_{3},p_{4}} Z_{1}^{p_{1}}Z_{2}^{p_{2}}Z_{3}^{p_{3}}Z_{4}^{p_{4}}, i = 1, …, 4  (5.4)

5.1. periodic solutions of Henon–Heiles’ system

The domain in phase space, where the normalizing transformation will be analytic for a system which has two couples of complex conjugated imaginary units as eigenvalues, is defined by equations (3.9):

λ_{i}Z_{i}a = λ_{i}Z_{i} + Z_{i} ∑ ; \begin{matrix}p_{i} ≥ −1 \\ p_{1}, …, p_{i−1}, p_{i+1}, …, p_{4} ≥ 0 \\ p_{1} + p_{3} = p_{2} + p_{4} \end{matrix} g_{i,p_{1},p_{2},p_{3},p_{4}} Z_{1}^{p_{1}}Z_{2}^{p_{2}}Z_{3}^{p_{3}}Z_{4}^{p_{4}}, i = 1, …, 4  (5.5)

where a is an expression, which does not depend on i index. It can be proved that a does not depend on time if (5.5) are satisfied.

The set in phase space of (5.4) for which (5.5) is fulfilled will produce by the inverse transformation to (5.3) the invariant periodic set of (5.2) that contracts to the stationary point (a zero point here) (Bruno, 1993).

Equations (5.5) are simplified using Groebner bases technique (Buchberger, 1985), (Boege, Gebauer, Kredel, 1986) with the REDUCE (Hearn, 1987) computer algebra system. We obtained 8 proper solutions of (5.5) and they produce 6 families of different periodic solutions of (5.1) with 4 different values of a frequency.

As 4 equations (5.5) are not independent we have as a solution a hyperplane rather than a point:

solution 1: Z_{1}Z_{4} = Z_{2}Z_{3}, Z_{3}^{2} = 3Z_{1}^{2};
solution 2: Z_{1}Z_{4} = Z_{2}Z_{3}, Z_{1}^{2} = 3Z_{3}^{2};
solution 3: Z_{1}Z_{4} = −Z_{2}Z_{3}, Z_{1}^{2} = Z_{3}^{2};
solution 4: Z_{3} = Z_{4} = 0;
solution 5: Z_{1} = Z_{2} = 0;
solution 6: Z_{2} = Z_{4} = 0;
solution 7: Z_{1} = Z_{3} = 0;

Corresponding frequences for the solutions of (5.1) are:

ω_{1,4} = 1 − \frac{5}{6}H + \frac{17}{48}H^{2} + \frac{127517}{38880}H^{3} + \frac{51952319}{3732480}H^{4} + \frac{1675438657}{111974400}H^{5}
ω_{2,5} = 1 − \frac{5}{6}H − \frac{95}{48}H^{2} − \frac{54935}{7776}H^{3} − \frac{22030445}{746496}H^{4} − \frac{200207485}{1492992}H^{5}
ω_{3} = 1 + \frac{1}{3}H − \frac{2}{3}H^{2} + \frac{5389}{2430}H^{3} − \frac{52393}{5832}H4 + \frac{29471957}{729000}H^{5}
ω_{6} = ω_{7} = 1

where H is the full energy:

H = \frac{1}{2}[(\frac{dX}{dt})^{2} + (\frac{dY}{dt} )^{2} + X^{2} + Y^{2}] + X^{2}Y − \frac{1}{3}Y^{3}

For each frequence the above mentioned the corresponding solution in analytical form (as a Fourier series in t and H) was obtained. This is the solution which corresponds to ω_{1}:

X =

\frac{\sqrt{3}}{2799360} cos(8ωt)ε^{8} +
\frac{71\sqrt{2}}{18662400} cos(7ωt)ε^{7} +
\frac{\sqrt{3}}{259200} cos(6ωt)ε^{6} (130ε^{2} + 9)+
\frac{17\sqrt{2}}{1555200} cos(5ωt)ε^{5} (341ε^{2} + 30)+
\frac{\sqrt{3}}{583200} cos(4ωt)ε^{4} (3451ε^{4} + 3204ε^{2} + 540)+
\frac{\sqrt{2}}{3110400} cos(3ωt)ε^{3} (1238051ε^{4} + 398160ε^{2} + 64800)+
\frac{\sqrt{3}}{1399680} cos(2ωt)ε^{2} (68930ε^{6} + 2152809ε^{4}
+933120ε^{2} + 233280)+
\frac{\sqrt{2}}{18662400} cos(ωt)ε (−51700693ε^{6} − 5862600ε^{4}
+4795200ε^{2} + 9331200)+
\frac{\sqrt{3}}{1555200}ε^{2} (8445851ε^{6} + 1072800ε^{4}
+21600ε^{2} − 777600)
Y =  −
\frac{1}{2799360} cos(8ωt)ε^{8} +
\frac{71\sqrt{6}}{18662400} cos(7ωt)ε^{7} −
\frac{1}{259200} cos(6ωt)ε^{6} (130ε^{2} + 9)+
\frac{17\sqrt{6}}{1555200} cos(5ωt)ε^{5} (341ε^{2} + 30)−
\frac{1}{583200} cos(4ωt)ε^{4} (3451ε^{4} + 3204ε^{2} + 540)+
\frac{\sqrt{6}}{3110400} cos(3ωt)ε^{3} (1238051ε^{4} + 398160ε^{2} + 64800)−
\frac{1}{1399680} cos(2ωt)ε^{2} (68930ε^{6} + 2152809ε^{4}
+933120ε^{2} + 233280)−
\frac{\sqrt{6}}{18662400} cos(ωt)ε (51700693ε^{6} + 5862600ε^{4}
−4795200ε^{2} − 9331200)−
\frac{1}{1555200}ε^{2} (8445851ε^{6} + 1072800ε^{4}
+21600ε^{2} − 777600)

Results of the calculation were verified by two ways. The first one was a direct substitution of the series in the original equations (5.1). All results of this procedure contained only terms with negligible orders. The second way was a comparison of the numerical solutions of equations (5.1) calculated by the Runge–Kutta method with the tabulated values of the approximate solutions at values of energy H = \frac{1}{24} and H = \frac{1}{12}. We evaluated the values of f_{err} which are maximum values of relative deviations of the numerical solutions from the approximate solutions in a phase space during one period. The results are:

\frac{d^{2}X}{dt^{2}} = −X − 2XY, \frac{d^{2}Y}{dt^{2}} = −Y − X^{2} + Y^{2}

-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5-0.4 -0.3 -0.2 -0.1  0  0.1 0.2 0.3 0.4 0.5^{Ẏ}  ^{Y}  ^{ω_{1}} ^{+}^{+}^{+}  ^{ω_{2}} ^{◇}^{◇}^{◇}  ^{ω_{3}} ^{□}^{□}  ^{□}  ^{ω_{4}} ^{×}^{×}  ^{ω_{5}}

Figure 2. Intersections of periodic solutions of Henon–Heiles’ system with plane X = 0, \frac{dX}{dt} = \frac{dX}{dt} (X = 0, Y, Y′, H) at the full energy \frac{1}{12} .

H = \frac{1}{24} :  H = \frac{1}{12} :
for solutions with ω_{1} :  f_{err} = 1.3 × 10^{−5} f_{err} = 7.6 × 10^{−4}
for solutions with ω_{2} :  f_{err} = 4.7 × 10^{−5} f_{err} = 4.4 × 10^{−3}
for solutions with ω_{3} :  f_{err} = 1.1 × 10^{−5} f_{err} = 5.0 × 10^{−4}
for solutions with ω_{4} :  f_{err} = 1.3 × 10^{−5} f_{err} = 7.6 × 10^{−4}
for solutions with ω_{5} :  f_{err} = 2.9 × 10^{−5} f_{err} = 1.8 × 10^{−3}

We did not check numerically the complex solutions with ω_{6}, ω_{7}.

In Fig.5.1 intersections of periodic solutions of Henon–Heiles’ system with plane X = 0, \frac{dX}{dt} = \frac{dX}{dt} (X = 0, Y, Y′, H) for the full energy \frac{1}{12} are displayed. Solution 5 lies in that plane and other solutions cross this plane in a couple points each except solution 4 which looks like a digit ”8” and crosses the plane X = 0 in a single point.

It is necessary to remark that for the value H = \frac{1}{6} we will have a chaotic regime of a behavior of the system (Henon, Heiles, 1964). So \frac{1}{24} and \frac{1}{12} are not physically small values.

The creation of the normal form and the normalizing transformation took about 30 secs and brought 110 terms to the form and 1250 terms to the transformation.

6. The program

The earlier attempts of the author to realize the creation of high enough orders of the normal form on algebraic level of REDUCE-2 system were not very successful as a sequence of small capacity of computers which the author used. And for realization of the described frame the NORT program (Edneral, Khrustalev, 1985, 1992), was created. The program NORT is written in Standard LISP and contains about 1500 LISP operators. It is a package of procedures to treat truncated multivariate power series of arbitrary dimensions. In addition to the procedures for arithmetic operations there are special functions for the creation of normal forms and functions for substitutions, for calculation of roots (when it is possible), for differentiating, for printing and for inverting of multivariate power series. Coefficients of the truncated power series may be treated in three different (precision rational, floating point or an approximate rational (Rodionov, in progress) ) arithmetics according to the user’s choice. Separately user may choose a representation of numbers (in a floating point or in a rational form) for output. This output has a REDUCE acceptable form, but REDUCE’s internal procedures are not called by the program. The program uses a recurrent polynomial representation separately in each order of the stored series. There is a just compact use of memory. Unfortunately at this moment the system has no friendly user interface and has a research character. More that for nonlinear systems of ODE there is a so wide variety that we have no clear ideas about proper realization of the interface as it was possible for second orders systems (pay your attention to the nice package ALKAHEST III (Fitch, Norman, Moore, 1986)).

The time of a garbage collection for above examples of second order equations was smaller that 3% of evaluation time. This can characterize the NORT system as a system with a good enough internal LISP organization.

Regarding a speed of calculations under different arithmetics. The evaluation ov van der Pol’s equation till 28 order took under the rational arithmetic 21 minutes, under approximate rational arithmetic it took 13 minutes with relative precision about 10^{−7} and under float point arithmetic the evaluation took 161 secs.

7. Conclusion

Here we can conclude that the creation of the high order normal forms allows to produce finite formulaes for a quantitative approximation of periodic solutions of nonlinear ODEs.

Acknowledgments:

I am very grateful to A.D.Bruno for important advices and useful discussions, to O.A.Khrustalev for a concentration of my attention on the normal form method at the most beginning, to J.H.Davenport and J.P.Fitch for support at the stage of further investigations. I would like also to express my acknowledgments to H.Melenk who kindly introduced for my needs the ”zeropoint” facility into the Groebner/REDUCE package and to S.Walcher who provided me with his newest papers.

References

Andersen, G.M., Geer, J.F. (1983). Power series expansions for the frequency and period of the limit

cycle of the van der Pol equation. SIAM J. Appl. Math. 42, 678–693. Arnold, V.I., Anosov, D.V. (Eds.) (1988). Dynamical Systems, I. Encyclopaedia of Mathematical Sci-

ences. Springer–Verlag, N.Y.

Boege, W., Gebauer, R., Kredel, H. (1986). Some Examples for Solving Systems of Algebraic Equations

by Calculating Groebner Bases. J. Symbolic Computation 1, 83–98. Bruno(Brjuno), A.D. (1971) Analytical form of differential equations. I. Trans. Mosc. Mat. Soc., 25,

131–288. Bruno(Brjuno), A.D. (1972) Analytical form of differential equations. II. Trans. Mosc. Mat. Soc. 26,

199–239. Bruno, A.D. (1989) Local Method in Nonlinear Differential Equations. Part I - The Local Method of

Nonlinear Analyses of Differential Equations, Part II - The Sets of Analyticity of a Normalizing Transformation. Springer Series in Soviet Mathematics. ISBN 3-540-18926-2, pp. 370.

Bruno, A.D. (1990) The restricted three-body problem: Plane periodic orbits. Moscow, Nauka. ISBN

5-02-000683-1, 295. (in Russian). Bruno, A.D. (1993) Bifurcation of the Periodic Solutions in the Case of a Multiple Pair of Imaginary

Eigenvalues. Selecta Mathematica formerly Sovietica 12, # 1, 1–12. Bruno, A.D., Walcher, S. (to be printed) Symmetries and convergence of normalizing transformations.

J. Math. Anal. Appl. Buchberger, B. (1985). Groebner Bases: An Algorithmic Method in Polynomial Ideal Theory. Progress,

directions and open problems in multidimensional systems theory. Dordrecht: Reidel, edited by N.K.Bose, 184–232.

Deprit, A. (1969) Canonical transformation depending on a small parameter. Celestial Mechanics 1, #

1, 12–30. Godziewski, K., Maciejewski,A.J. (1990) Celest. Mech. Dyn. Astron. 49, 1. Ito, H. (1989) Convergence of Birkhoff normal forms for integrable systems. Comment. Math. Helv. 64,

412–461. Ito, H. (1992) Integrability of Hamiltonian systems and Birkhoff normal forms in the simple resonance

case. Mth. Ann. 292, 411–444. Edneral, V.F., Khrustalev, O.A. (1985). The Normalizing Transformation for Nonlinear Systems of ODE.

The Realization of the Algorithm. Proceedings of International Conference on Computer Algebra and its Application in Theoretical Physics (USSR, Dubna, September 1985). Dubna: JINR publ., 219–224. In Russian.

Edneral, V.F., Khrustalev, O.A. (1992). Program for Recasting ODE Systems in the Normal Form. Sov.

J. Programmirovanie, # 5, 73–80. In Russian. Edneral, V.F. (1993). Computer Generation of Normalizing Transformation for Systems of Nonlinear

ODE. Proceedings of the 1993 International Symposium on Symbolic and Algebraic Computation (Kiev, Ukraine, July 1993). New York: ACM Press, edited by M.Bronstein, 14–19.

Fitch, J.P., Norman, A.C., Moore, P.M.A. (1986). ALKAHEST III: Automatic analysis of periodic weakly

nonlinear ODEs. Proceedings of SYMSAC 86. Waterloo, 34—38. Guckenheimer, J., Holmes, P. (1986) Nonlinear Oscillations, Dynamical Systems and Bifurcations of

Vector Fields. Springer –Verlag, N.Y. Hassard, B.D., Kazarinoff, N.D., Wan, Y.-H. (1981) Theory and Applications of Hopf Bifurcation. Cam-

bridge Univ. Press. 280. Hearn, A.C. (1987). REDUCE. User’s manual. Rand Publication CP78. Henon, M., Heiles,C. (1964). The Applicability of the Third Integral Of Motion: Some Numerical Ex-

periments. Astronomical J. 69, 73–79. Hori, G.I. (1966).Theory of general perturbations with unspecified canonical variables. J. Japan Astron.

Soc. 18, # 4, 287–296. Mersman, W.A. (1970) A new algorithm for Lie transformation. Celestial Mechanics 3, # 1, 81–89. Rand, R., Armbruster D. (1987) Perturbation methods, Bifurcation Theory and Computer Algebra.

Springer–Verlag, N.Y. Rodionov, A.Ya. (in progress). A Program for Rational Arithmetic RATAR. (In Russian). Shevchenko, I.I., Sokolsky, A.G. (1993) Algorithms for normalization of Hamiltonian systems by means

of computer algebra. Comp. Phys. Comm. 77, 11—18. Vallier, L. (1993) An Algorithm for the Computation of Normal Forms and Invariant Manifolds. Proceed-

ings of the 1993 International Symposium on Symbolic and Algebraic Computation (Kiev, Ukraine, July 1993). New York: ACM Press, edited by M.Bronstein, 225—233.

Walcher, S. (1991) On differential equations in normal form. Math. Ann. 291, 293—314. Walcher, S. (1993) On Transformations into Normal Form. J. Math. Analysis and Appl. 180, 617—632.