Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://acat02.sinp.msu.ru/presentations/edneral/EdneralACAT2002.ps Äàòà èçìåíåíèÿ: Tue Jul 16 17:44:20 2002 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 20:30:36 2012 Êîäèðîâêà:

Application of the Resonant Normal Form to High Order
Nonlinear ODEs Using MATHEMATICA
Victor Edneral
Institute for Nuclear Physics,
Moscow State University, Moscow, Russia
edneral@theory.sinp.msu.ru
Raya Khanin
Department of Applied Mathematics and Theoretical Physics,
University of Cambridge
R.Khanin@damtp.cam.ac.uk
July 16, 2002
Abstract
This paper studies applications of the normal form method to dynamical systems
using Mathematica. It describes the package for power series treatment, and the pro-
cedures for building a normal form and normalizing transformations. Examples of
treatment of a crude one-dimension equation and a double pendulum system are given.
1 Introduction
The concept of using co-ordinate transformation to simplify ordinary dierential equations
has been widely used for a long time. The normal form method which has applications in
many branches of engineering, physics, and applied mathematics is a powerful technique
to simplify the governing ordinary dierential equations. The basic idea of the method of
normal form is the use of \local" coordinate transformation to simplify the equations that
describe the dynamics of the system under consideration. In other words, with the method
of normal forms, one seeks a near-identity coordinate transformation in which the dynamical
system takes the simplest form or so-called normal form. The transformations are called
local because they are generated in a neighbourhood of a known solution such as a xed
point or a periodic orbit of the system. For nonlinear dynamical systems the coordinate
transformation will, in general, be nonlinear. The normalisation is usually carried out with
respect to perturbation parameters. The normal forms of dierential equations are also very
useful in bifurcation analysis [1].
The method of normal forms dates back to the times of Euler, Poincare, Dulac and
Birkho. In this paper, algorithm based on the approach by Bruno [2] is employed. It is an
algorithm for a resonant normal form creation. The important advantage of this approach
is its general algorithmic frame, which allows one to investigate a wide class of autonomous
systems including Hamilton and non- Hamilton cases. The algorithm provides a constructive
way for obtaining the approximations of local families of periodic and conditional periodic
solutions in the form of power (Fourier) series. Another advantage of the approach employed
1

here is the simplicity of creation of the normal form and the corresponding transformations,
which is implemented using recursive formula.
Examples of the lowest orders calculation of the normal forms in computer algebra nec-
essary for bifurcation analysis has been done in [3]. A package for calculating coeôcients
of the normal form has been developed by one of the authors [4]. The package NORT is
written in Standard LISP and contains about 2000 operators.
2 Preliminary treatment of the system
Consider a system of autonomous ordinary dierential equations:
_
x = (x); (1)
where x = (x 1 ; : : : ; xn ) is a vector function of time, (x) is a vector function of x and
other parameters. Such equations arise in numerous scientic, engineering and industrial
applications. It is worth noting that a non-autonomous system of equations
_
x = ^
(x; t): (2)
can be transformed into the autonomous system by introducing an additional variable  =
(t) such that ^
(x; ) becomes (explicitly) independent of t. (t) is dened by adding the
equation _
 = () to system (2). For example cos(!t) can be replaced by the sum ( 1 + 2 )=2
with an additional couple of equations: _
 1 = i
!
 1 and _
 2 = i
!
 2 .
The important condition on a system of ordinary dierential equations (1) is that it
is explicitly solved with respect to the rst derivatives x(t). Sometimes this can be done
approximately.
2.0.1 Preanalysis.
The second important condition on system (1) is that it is studied in the neighbourhood
of stationary point x 0 , (x 0 ) = 0. It is convenient to shift variables by x 0 , so that
(0) = 0. Sometimes the original system has no stationary points or it is important to
investigate it in the domain which is far from such points. In the neighbourhood of regular
points other methods should be used (see for example [5]). Usually, the right-hand side of
the system is approximated by polynomials. For analytical function, it is simply truncated
power series of this function. Finally, the linear part of the right-hand side is transformed
to Jordan's form so that the original system (1) is represented by
_
y i =  i y i +  i y i 1 + ~
 i (y);  1 = 0; i = 1; : : : n (3)
Here  = ( 1 ; : : : ; n ) is the vector of eigenvalues of the matrix of the linear part of the
system and ~
 is the vector of polynomials obtained from approximating  by series of
polynomials. ~
 does not contain constant or linear terms. Constant terms have been
removed by shifting into the stationary point while linear terms are in the linear part of the
system, which is written in Jordan form.
2.0.2 Normal form transformation.
For the clarity of exposition let us consider pure diagonal systems. (The general normal
form method can be applied to a full Jordan case as well.) Equations (3) can be written in
the form:
_
y i =  i y i + y i
X
q2N i
f i;q y q ; i = 1; : : : ; n ; (4)
2

where we use the multi-index notation:
y q =
n
Y
j=1
y j q j ;
with the power exponent vector q = (q 1 ; : : : ; q n ) and the sets:
N i = fq 2 Z n
: q i  1 and q j  0 ; if j 6= i ; j = 1; : : : ; ng ;
because the factor y i has been moved out of the sum in (4).
The normalization is done with a near-identity transformation:
y i = z i + z i
X
q2N i
h i;q z q ; i = 1; : : : ; n (5)
and then equation (4) will have the form:
_
z i = i (z) def
=  i z i + z i
X
hq; i = 0
q 2 N i
g i;q z q ; i = 1; : : : ; n : (6)
The important dierence between (4) and (6) is a restriction on the range of the sum-
mation which is dened by the equation:
hq; i def
=
n
X
j=1
q j  j = 0 : (7)
The h and g coeôcients in (5) and (6) are found by using the recurrence formula:
g i;q + hq; i
h i;q =
n
X
j=1
X
p + r = q
q 2 N i
(p j + ô ij )
h i;p
g j;r + ~
 i;q ; (8)
where the second summation in the right-hand side is over all integer vectors satisfying the
constraint p + r = q, and ~
 i;q is a coeôcient of the factor z i z q in the polynomial ~
 i in
(3), arguments of which have been transformed by (5). Here jjpjj and jjrjj < jjqjj, where
jjqjj def
= q 1 + : : : + q n .
The ambiguity in (8) is usually xed by the conventions:
h i;q = 0; if hq; i = 0;
g i;q = 0; if hq; i 6= 0 (9)
and then the normalizing transformation is called a \basic".
3 Package for working with series
As one of the steps in the preliminary analysis involves approximating right-hand side by
polynomials, one needs an eective tool for dealing with power series. Surprisingly, Math-
ematica lacks tools for working with multivariate power series. To ll this gap, a package
PolynomSeries has been developed by the authors. The name of the package arises from an
idea of a transformation between polynomial and specic power representation in Mathemat-
ica. Below is a brief description of the functionalities of the package. Consider a multivariate
polynomial
3

p = Plus@@Flatten[Table[y1^i*y2^j*y3^k,
{i, 0, 1}, {j, 1, 2}, {k, 1, 3}]]
y2 y3 + y1 y2 y3 + y2 2 y3 + y1 y2 2 y3 + y2 y3 2 + y1 y2 y3 2 + y2 2 y3 2 + y1 y2 2 y3 2
Transforming the polynomial into a series with respect to variables y 1 ; y 2 ; y 3 yields:
R = Polynom2Series[p, {y1, y2, y3}]
ff2; y2 y3g; f3; y1 y2 y3+y2 2 y3+y2 y3 2 g; f4; y1 y2 2 y3+y1 y2 y3 2 +y2 2 y3 2 g; f5; y1 y2 2 y3 2 gg
The above list, R, is the main object of the package. Each element of R consists of two
entries: the rst one is dened as \index". It is the common power  =  1 +  2 : : : n
of the monomial y 1
1 y 2
2 : : : y n
n . The monomial itself is the second entry. Here y j are the
elements of the second input parameter of the function. One can consider transformation of
the polynomial into the series with respect to some (not all) variables (y 2 ; y 3 below):
Polynom2Series[p, {y1, y3}]
ff1; y2 y3+y2 2 y3g; f2; y1 y2 y3+y1 y2 2 y3+y2 y3 2 +y2 2 y3 2 g; f3; y1 y2 y3 2 +y1 y2 2 y3 2 gg
Several arithmetic operations on series of type R are dened in the package, including
addition of two series (AddR[R1, R2]), multiplying two series and retaining terms up to or-
der (index) (MultRR[R1, R2, n]) or returning only terms of index n (MultRRLayer[R,
R1, n]), ordering series by their indices (OrderR[R])
It is often in the applications that one needs to substitute a new series for some (or
all) variables in the given series. This can be done using SubstR[R, RV, xvar, yvar,
n] command. Here RV is a vector of series with respect to new variables yvar which are
substituted into the series R with respect to variables xvar. To substitute a simple series
(z1+z2; z2+z3; z1+z3) to every (y1; y2; y3) in the series R, a vector of series substitutions
is to generated rst:
subst = {z1 + z2, z2+z3, z1+z3};
RVy = MapThread[Polynom2Series[#1, #2]&,
{subst, {{z1, z2, z3}, {z1,z2,z3}, {z1,z2,z3}}}]
fff1; z1 + z2gg; ff1; z2 + z3gg; f1; z1 + z3ggg
Rxy = SubstR[R, RVy, {y1, y2, y3}, {z1, z2, z3}, 3]
ff2; z1z2 + z1z3 + z2z3 + z3 2 g; f3; 2z1 2 z2 + 2z1z2 2 + 2z1 2 z3 + 6z1z2z3 (10)
+ 2z2 2 z3 = 4z1z3 2 + 4z2z3 2 + 2z3 3 gg
Note that we only retain terms up to index 3. If the user is interested only in terms of index,
say, 4, SubstRLayer[ ] is to be used:
Rxy = SubstRLayer[R1, RVy, {y1, y2, y3}, {z1, z2, z3}, 4]}
ff4; z1 3 z2 + 3z1 2 z2 2 + z1z2 3 + z1 3 z3 + 7z1 2 z2z3 + 7z1z2 2 z3 + z2 3 z3 + 4z1 2 z3 2
(11)
+ 10z1z2z3 2 + 4z2 2 z3 2 + 4z1z3 3 + 4z2z3 3 + z3 4 gg
There are separate commands for two commonly used functions, which are singular in
zero:square root of series and inverse function. The correct usage of these commands can
only be guaranteed if the series contains a constant term (i.e. term with index 0). To nd
square root with terms up to index 4 use
4

SqrtR[R, 4]
ff0;
p
7; f1; x12
p
7g; f2; 27 x1 2 56
p
7g; f3; 365 x1 3 784
p
7g; f4; 19763 x1 4 43904
p
7gg
Inverse series (1=R) is found by
InverseR[R, 4]
ff0; 17g; f1; x149g; f2; 6 x1 2 343g; f3; 36 x1 3 2401g; f4; 216 x1 4 16807gg
4 Mathematica NormalForm package
The procedures for building a normal form and a normalizing transformation till the xed
order are implemented in the package NormalForm. It includes procedures for a linear
normalization, for a nonlinear normalization and for evaluating a periodicity condition. The
input for investigation is the system of ODEs.
The package was tested by comparing the results for approximated solutions of Duông's,
van-der-Pol's, Henon-Heiles' systems with solutions obtained before by V.Edneral with the
LISP based package NORT [4, 7]. The results are fully identical.
4.1 Example of one-dimension system.
To demonstrate the main functions of the package, let us consider a trivial example of a
one-dimension crude system:
_
y = y + y 2 : (12)
This equation has an exact solution
y(t) = c
exp( t)
1 + c
exp( t)
; (13)
where c here is an integration constant.
Now let us calculate an approximation of this solution by the normal form method using
the package NormalForm. The equation (12) is entered in standard form:
Eqs = f _
y[t] == y[t] + y[t] 2 g;
The rst step of a linear normalization is performed by the procedure LinearNormaliza-
tion[Eqs List, OVars List, NVars List]. Here Eqs is a list of equations resolved with
respect to the rst derivatives of variables, OVars (original variables) and NVars is a list
of the new variables:
{Atr, Btr, Jrd, NonLinV} = LinearNormalization[Eqs, {y[t]}, {u[t]}]
ff1g; f1g; f1g; fu[t] 2 gg
The result is four matrices: Atr and Btr are diagonalizing (right and inverse) matrices
and Jrd is a matrix of the linear part in Jordan's form, i.e. Jrd = Atr . LinearPart .
Btr, where LinearPart is a matrix of coeôcients of the linear part of the original system.
NonLinV is a vector of nonlinear parts in new variables
u[t] = Atr:fy[t]g:
The system is in a diagonal form already, so we have an identical transformation. Now we
look at the linear part matrix. It is a 1 by 1 matrix and it has a single element y, i.e. the
5

eigenvalues vector is -1. Because scalar production (7) can not be zero in this case, the right
hand side of the normalized equation (normal form) will be linear:
_
z[t] = z[t]: (14)
It is remarkable that we can write down the normal form without any calculations. It can
be done for any eigenvalues for which equation (7) has nite numbers of solutions in natural
integers. Such cases are called crude ones. Let us use the NormalForm procedure for the
corresponding calculation:
{GRv, uN} = NormalForm[Jrd, NonLinV, u, z, 3]
ff1g; fz1[t] z1[t] 2 + z1[t] 3 z1[t] 4 gg
The result above (4:1) consists of a vector of the normalized right hand side GRv and
a normalizing transformation uN (till the third order over identical transformation, i.e.
3 + 1 = 4) such that the normalized equation has the form:
_
z i [t] = GRv i
z i [t]; i = 1; : : : ; n (15)
y i [t] =
X
j=1;n
Atr i;j u j [t] =
X
j=1;n
Atr i;j uN j ; i = 1; : : : ; n (16)
Here n is a dimension of the system, which is 1 in our case.
Solution of (14,4:1,15) is z(t) = c
exp( t). For an approximation of a solution of (12)
we need now to substitute this solution into the corresponding normalizing transformation
of the type (5,4:1,16). After that we will have a truncated expansion of (13) in the series of
constant c:
y[t]  c
exp( t) [c
exp( t)] 2 + [c
exp( t)] 3 [c
exp( t)] 4
4.2 Example of four-dimension system. Double pendulum.
Consider now a system of the 4-th order which describe equations of motion for the double
pendulum:
( m 1
3 +m 2 )l 2
1 + m 2 l 1 l 2 cos  2 + m 2 l 2
2
3 ) 
 1 + ( m 2 l 1 l 2
2 cos  2 + m 2 l 2
2
3 ) 
 2 (17)
= m 2 l 1 l 2
2 ( _
 1 + _
 2 ) 2 sin  2
 m 1
2 +m 2

gl 1 sin  1
m 2 l 1 l 2
2
_
 2
1 sin  2
m 2 gl 2
2 sin( 1 +  2 ):
( m 2 l 1 l 2
2 cos  2 + m 2 l 2
2
3 ) 
 1 + m 2 l 2
2
3
  2 (18)
= m 2 l 1 l 2
2
_
 2
1 sin  2
m 2 gl 2
2 sin( 1 +  2 ):
Equations in this form have been created by the automatic equation generator by D.Forehand
et al [8]. This system is not initially solved with respect to derivatives. Assuming that 
and  are small and approximating Sin up to the third order, a system in polynomials with
respect to unknown variables is obtained (not shown). Solving the system with respect to
higher derivatives and expanding the denominators in the right hand side in the power series
6

of  and  till the third order yields the following right-hand sides for the system of the
fourth order with respect to new variables 1[t] = [t], 2[t] =  0 [t], 1[t] = [t], 2[t] =  0 [t]:
NewEqs =

(19)
 0
2 [t] == 18 1 [t]
7 + 183 1 [t] 3
49
9
7  1 [t] 2 [t] 2 + 9 1 [t]
7
873
98  1 [t] 2
 1 [t]
+ 9
7  2 [t] 2  1 [t] + 387
49  1 [t] 1 [t] 2 123 1 [t] 3
49
6
7  1 [t] 2 [t] 2 + 6
7  1 [t] 2 [t] 2
 0
2 [t] == 27 1 [t]
7
369 1 [t] 3
49 + 24
7  1 [t] 2 [t] 2 24 1 [t]
7 + 891
49  1 [t] 2  1 [t]
24
7  2 [t] 2  1 [t] 1539
98  1 [t] 1 [t] 2 + 244 1 [t] 3
49 + 9
7  1 [t] 2 [t] 2 9
7  1 [t]2[t] 2
 1
0 [t] ==  2 [t];
 1
0 [t] ==  2 [t]

An assumption that the lengths and masses of two pendulums are equal l21 = 1; m21 = 1
has been made to simplify computations: The system is then recasted in the Jordan form.
Float point arithmetic was used for calculation of radicals in eigenvalues. The calculation
of the normal form up to the 5-th order yields a vector of the normalized right-hand sides:
 i + 0:13 iz 1 [t]z 2 [t] 0:08 iz 1 [t] 2 z 2 [t] 2 0:08 iz 3 [t]z 4 [t] (20)
+ 0:21 iz 1 [t]z 2 [t]z 3 [t]z 4 [t] 3:65 iz 3 [t] 2 z 4 [t] 2 ;
0:85 i 0:13 iz 1 [t]z 2 [t] + 0:08 iz 1 [t] 2 z 2 [t] 2 + 0:08 iz 3 [t]z 4 [t]
0:21 iz 1 [t]z 2 [t]z 3 [t]z 4 [t] + 3:65 iz 3 [t] 2 z 4 [t] 2 ;
2:29 i 0:30 iz 1 [t]z 2 [t] + 1:52 iz 1 [t] 2 z 2 [t] 2 + 7:87 iz 3 [t]z 4 [t]
+ 6:22 iz 1 [t]z 2 [t]z 3 [t]z 4 [t] + 71:86 iz 3 [t] 2 z 4 [t] 2 ;
2:29 i + 0:3 iz 1 [t]z 2 [t] 1:52 iz 1 [t] 2 z 2 [t] 2 7:87 iz 3 [t]z 4 [t]
6:22 iz 1 [t]z 2 [t]z 3 [t]z 4 [t] 71:86 iz 3 [t] 2 z 4 [t] 2

In new variables z1[t]; : : : ; z4[t] the system has the form (6). After the standard treatment
of the normalized system [4] we get approximations for conditionally periodic families of
solutions, where a; b; t1; t2 are integration constants below (the result is truncated at order
3 in small a; b):
[t] = (21)
1:4a cos[(t + t 1 )! 1 ] + 0:09a 3 cos[(t + t 1 )! 1 ] 1:80ab 2 cos[(t + t 1 )! 1 ]
+ 0:65a 3 cos[3(t + t 1 )! 1 ]
0:95b cos[(t + t 2 )! 2 ] + 0:07a 2 b cos[(t + t 2 )! 2 ] 1:65b 3 cos[(t + t 2 )! 2 ]
+ 0:8b 3 cos[3(t + t 2 )! 2 ]
+ 0:8ab 2 cos[t! 1 + t 1 ! 1 2t! 2 2t 2 ! 2 ] 0:24a 2 b cos[2t! 1 + 2t 1 ! 1 t! 2 t 2 ! 2 ]
+ 0:61a 2 b cos[2t! 1 + 2t 1 ! 1 + t! 2 + t 2 ! 2 ] + 0:65ab 2 cos[t! 1 + t 1 ! 1 + 2t! 2 + 2t 2 ! 2 ]
[t] = (22)
2:a cos[(t + t 1 )! 1 ] + 0:18a 3 cos[(t + t 1 )! 1 ] + 3:57ab 2 cos[(t + t 1 )! 1 ]
1:39a 3 cos[3(t + t 1 )! 1 ]
7

+ 2:b cos[(t + t 2 )! 2 ] 0:12a 2 b cos[(t + t 2 )! 2 ] + 3:39b 3 cos[(t + t 2 )! 2 ]
1:59b 3 cos[3(t + t 2 )! 2 ]
2:16ab 2 cos[t! 1 + t 1 ! 1 2t! 2 2t 2 ! 2 ] + 3:28a 2 b cos[2t! 1 + 2t 1 ! 1 t! 2 t 2 ! 2 ]
1:12a 2 b cos[2t! 1 + 2t 1 ! 1 + t! 2 + t 2 ! 2 ] 0:97ab 2 cos[t! 1 + t 1 ! 1 + 2t! 2 + 2t 2 ! 2 ]
where frequencies are given by the following formulas
!1 = 0:85 + 0:08a 4 + 0:08b 2 + 3:65b 4 + a 2 0:13 0:21b 2


(23)
!2 = 2:3 1:52a 4 7:87b 2 71:86b 4 + a 2 0:30 6:22b 2


(24)
To generate the periodicity condition ([2], [4]) the procedure ConditionForPeriodici-
ty[GRv, Jrd, Var List] is employed. The meaning of GRv and Jrd arguments are
explained in the previous section and Var is a list of all small variables. In this example it is
fz1[t]; z2[t]; z3[t]; z4[t]g. Approximations (21) are conditionally periodic and they formally
satisfy equation (19), but the procedure ConditionForPeriodicity results the condition
for an orbital stability. A separate paper discussing this subject is in preparation. In the
case of the double pendulum case, this procedure results in a condition:
0:564279z1[t]z2[t] + 1:48528z1[t] 2 z2[t] 2 + 0:579792z1[t] 3 z2[t] 3 + (25)
6:91055z3[t]z4[t] + 4:84095z1[t]z2[t]z3[t]z4[t] + 58:0101z1[t] 2 z2[t] 2 z3[t]z4[t] +
69:8654z3[t] 2 z4[t] 2 + 268:305z1[t]z2[t]z3[t] 2 z4[t] 2 + 746:679579792z3[t] 3 z4[t] 3 = 0
This condition can be solved in a; b constants giving:
a = b(3:49953 + 89:1055b 2 + 8868:38b 4 ):
5 Summary
The authors revisit the normal form method and discuss its implementation in Mathematica.
The high order normal form for systems of ordinary dierential equations can be calculated.
A separate package enabling one to work with truncated multivariate power series is
discussed. The demo version of PolynomialSeries package can be accessed on
www.mathsource.com. Authors are going to develop new features of the package. The power
series package should use the O(x n ) mechanism. Truncation of power series allows one to
control the length of truncated series by a correct way as it is done in Maple. Therefore an
existing version is enough for a support of a normal form method. An important procedure
for solving equations in power series is planned to be added to the package [6].
The method produces approximate formulas for the approximation of crude and peri-
odic solutions of autonomous nonlinear ODEs. The comparison of Mathematica package
with an earlier version of normal form package (NORT) written in LISP demonstrates that
the calculations within the Mathematica system are more exible and convenient but are
considerably slower then under the LISP.
6 Acknowledgement
This research was supported by London Mathematical Society grant 02-5717 and the RF-
BR program for supporting leading scientic schools, grant 00-15-96-577. Raya Khanin is
Dorothy Hodgkin Royal Society Research Fellow.
8

References
[1] J.Guckenheimer, P.Holmes Nonlinear Oscillations, Dynamical Systems and Bifur-
cations of Vector Fields. Springer{Verlag, 1986, N.Y.
[2] A.D.Bruno Local Method in Nonlinear Dierential Equations, Springer-Verlag, Berlin
(1989)
[3] R.Rand, D.Armbruster Perturbation methods, Bifurcation Theory and Computer
Algebra. Springer{Verlag, 1987, N.Y.
[4] V.Edneral About Normal Form Method. Proceed. of the Second Workshop on Com-
puter Algebra in Scientic Computing (CASC'99, Munich, Germany, 1999), ed.by
Ganzha et al.,Springer, pp. 51{66.
[5] B.J.Dup ee, J.H.Davenport An Automatic Symbolic-Numeric Taylor Series ODE
Solver. Proceed. of the Second Workshop on Computer Algebra in Scientic Computing
(CASC'99, Munich, Germany, 1999), ed.by Ganzha et al.,Springer, pp. 37{50.
[6] V.F.Edneral, N.N.Vassiliev Approach to solving equations in the ring of formal
power series. Proceedings of XIIth Workshop on High Energy Physics and Quantum
Field Theory (Samara, Russia, 1997), ed. by B.B.Levtchenko, Moscow, 1999, p.462{465.
[7] V.F.Edneral V.F. A symbolic approximation of periodic solutions of the Henon{
Heiles system by the normal form method. J.Mathematics and Computers in Simulation,
Elsevier, v. 45, pp.445{463. Edited by A.Bruno, V.Edneral, S.Steinberg.
[8] D.I.M.Forehand, M.P.Cartmell, and R.Khanin Initial development towards an
integrated symbolic-analytical multibody code (submitted)
9