·

.

.~

~-

-

-

. -

-

~

~

-

~ -~

~.

-

.

' .

.

CD

Hamiltonian Study of the Toroidal
Helical Magnetic Field in the Vacuumwith Symplectic Integrators

A dissertation submitted in partial fulfillmentof the requirement for the degree of Doctor of Science

Giovanni GNUDI

Department of Material Science, Hiroshima University
Higashi- Hiroshima, Kagamiyama, Japan

January 1994

· ',. .

-.

~- .... - --- . .

~.

.. - - .

Abstract
It is well known that the magnetic field bnes system can be tr ated as a IIamilLoniansystem, however a full exploitation of this fact has been hindered by the difficulty offinding the proper canonical variables or by their inconv ni nce for th solution of practicalproblems. There would be many advantages in having a HamiHonian formulaLion of the
problem of the magnetic field lines flow, because for Hamiltonian systems many th Ol'eticai
results, like for example conservation theorems, are known. Moreover, in recent years the socalled symplectic numerical integration schemes have been developed, and th se schemes
are particularly designed for Hamiltonian dynamical systems. These methods are very
important for the magnetic field lines problem, because to obtain a detailed description of
the field lines system it is necessary to follow their winding around the torus many times,
typically for 10 5 rv 106 revolutions, and this results, when traditional integration methods
like the Runge-Kutta are used, in the accumulation of secular errors that make the results
unreliable. On the other hand, symplectic methods are free from secular errors.
In this thesis we investigate possible Hamiltonian formulations for the toroidal helicalmagnetic field in the vacuum and study the obtained models by means of a new linear
symplectic integration scheme developed by us.
We first derive a Hamiltonian for the cylindrical limit approximation of the field. The
Hamiltonian is derived until the second perturbative order and it is obtained by using a
perturbation theory based on the Lie transform. This perturbation theory is a theory inwhich all the relations are expressed in operational form, resulting in superior efficiency
and compadness when compared to the classical perturbation theory. The Lie transform." in the literature are intended for quite special Hamiltonian systems and are therefore notuseful in our case. The ideas underlining lts const.ruction are explain d, and stability andaccuracy analyses are also presented. Finally w present the numerical l' suIts, showing thatthe third order Hamiltoninan is able to fully reproduce the r gular and chaotic regions of
the toroidal helical field. The present study concludes that the IIamiltonian approach isfeasible and that a sui table Hamiltonian integrated by a sympl dic method can b a cheap
way to study with good accuracy the transition region between order and chao.

based perturbation theory forms the core of our theoretical study of all the Hamiltonianmodels treated in this thesis.
We th n study the general toroidal case and derive for it two different Hamiltonianmodels, th fir t on correct until the third perturbative order, the second one correct until
the s condo The second model is an integrable model, while the first one is not. Thisapparent discrepancy is connected with the use of the perturbation expansion and can be
xplained by noticing that the canonical variables used for the two models are different.
rvforeov r, the perturbation th ory yields expressions that are only asymptotic, so some caremust be x rcis d when u ing it. The second, iIltegrable model, of course is not able tor produce the chaotic region of the real field, but in the first model we are able to show thatth charad 1'i tic formation of magnetic surfaces, islands and chaotic regions is reproduced.
\Ne tb r fore proceed to th

d velopment of a new linear symplectic integration method.

Thi has be n motivat d by th ob ervation that most of the symplectic methods existing
11

111

·

.

-..

. - ...

- --

.

..

..

-. ~ - .

Acknow ledgell1en ts

It is my great pleasure and honor to acknowledge Professor Kyoji Nishi kawa and Professor
Tetsuya Sato for giving me the opportunity to study in this beauti I'll I country of Ri singSun. It was thanks to their invitation lelLers that I was abl

to obtain the Monbu ho

Scholarship that financially supported my three years and a half in Japan and without
which it would have been impossible to make a living here. Prof ssor

Ishikawa was thefirst to welcome me, a gaijin, in Hiroshima, and I will never forget hi s gr at kindness andgenerosity. Despite his fully booked schedule he always took time to listen to my probl ms
and to encourage me . I will always remember the six months I spent in Hiroshima as oneof the best times I had.
Thanks are given to all the students and the staff of Nishikawa-ken.
My three years at NIFS have been very fruitful. My deepest gratitude goes to Professor Tadatsugu Hatori, who constantly pushed and encouraged me, and who, throughuncountable discussions, shaped and deepened my understanding of the physics of Hamiltonian systems. He was always keen to listen to my questions and to help me in the manyIV~ _ <O""OO'~ --- - -

- --

-

-

-

-

v

problems that I faced, both scientific and of daily life. My deepest gratitude goes also toProfessor Tsuguhiro Watanabe, who was of invaluable help for deepening my understanding of numerical integration and spent much time in answering to my many naive questions.
\iVithout them the completion of this thesis would not have been possible.
Thanks also to all members of the Theory and Computer Simulation Center of NIFS for
their friendships and supports.

VI

-

-

--~

-

-

Vll

Contents

Abstract

Acknowledgementsv

1 Introduction

3

1.1

The Reasons of this V\Tork3

1.2

Overview . . . .

6

1.2.1

Notation

7

2 The Lie Transform

9

3

Vlll

===----==--- ~-~

-

---

2.1

Introduction . . .

92.2

Hamiltonian Mechanics .

10

2.3

The Lie Transform . . .

17Magnetic Field Lines Hamiltonian

25

3.1

25

Introduction . . . . . . . . . . . .
IX

.

~

~

-.

-

. -

.- - ~~.- -- - --.. --

-.

~

-

~

3.2

Magnetic Potential . . . . . . . .

26

3.3

Taylor Expansion of the Potential

28

6Concluding Remarks

85

87

A Coefficients
3.4
3.5
3.6
3.7

4

5

Cylindrical Limit Approximation..... .... . .... . . .. .. ...

31............... . . .... .. ....

40........ . . .... . ... . . . . . . . .......

45

Conclusion . . . . ... . . . . . . . .. .. . .................

47

Helical Toroidal Pot('ntial
Integrable Model

Symplectic Integration

49

4. 1

Introduc tion . . . . . . .. .................. . ...... ..

49

4.2

Symp lectic Me thods. . . . . . . . . . . . . . . . . . . . . ... . ......

52

4.3

The Di fferen ce Scheme . . .. . . . . . . . . . . . . . . . . . .. . . . . . .

54

4.3.1

Implem ent ation............. . . .. .. .........

57

4.3. 2

Truncation Error. . .............. . . . . . .......

57

4.4

Accuracy and Stabi li ty of the Me thod. .... . ....... ... .. . ..

58

4.5

Runge- K u t ta Form ulation. · ...... . . .... . .... . .... . . .

66

4.6

Linear Symplectic :McL hods

· . . . . . . . . . . . . .. . ... . . ... . .

70

4.7

Concl usion . . . . . .. . . . · .................. . . . ... .

72

Numerical Results
5.1

Introdu ction ..

. . . . . ................. . . . .....

75... . .......... . ......

76

x

~

-

-

87

A.2 T he case n =1

88

. . .. . .

A.3 The cases n = 2, 3, 4, 5, 6

- -

-

89

A.3 .2 n = 3

90

A.3 .3 n = 4

91

A.3 .4 n = 5

92

A.3.5

95

n = 6

R eferences

97

Xl

-

-

-

89

A.3. 1 n = 275

5.2 Toroi dal H amil tonian_~~T

A .I Int rodu ction

1

2

1

~-~-----

~

-

-

-

-

CHAPTER 1. INTRODUCTION

4with the magnetic chaos domain and the magnetic surfaces domain discussed above.

The existing research in thi field is based on numerical calculation about the tructureof the magnetic field, but in these methods numerical error are inherently present andChapter 1

this poses limitations on the accuracy by which the properties of the out.ermost magnetic
surface are investigated. To enhance the accuracy it is n cessary t.o track many times
the magnetic field lines around the torus, typically for 10 5Introduction"-J

10 6 l'('volutions, but there

are indications that this procedure causes a shrink of th surfa e towards th magnetic
axis and that the accumulat.ion of secular errors mak s the interpretation of th r suI ts
problematic. A possible explanation of this fact is that the divergenc -free natur of 1,h -'1.1magnetic field is not fully respected in thes num rical methods. In other word, many

The Reasons of this Worknumerical calculations do not respect the intrinsic Hamiltonian natnr(' or the magnetic field
The toroidal helical magnetic field possesses magnetic surfaces in the proximity of ,thelines system, in the sense that the numerical schemes utilized are not symplectic, while themagnet.ic axi , while in regions far from the magnetic axis such surfaces do not exist. The

flow of the magnetic field lines system is symplectic. With this background in mind, we

domain in which magn tic surfaces do not exist is called magnetic chaos domain, andthink that in order to further advance this kind of research it is necessary to approach the

betw 'n Lh magnetic axis and the magnetic chaos domain lies the outermost magnetic

problem from the point of view of Hamiltonian theory, in which the div rgenc -fr e property

surface, whose position is an important information from the point of view of magnetic

is rigorously respected. Besides, it is important to notice that for Hamil tonian syst ms it isconftn ment.

possible to use symplectic int gration techniques to solve th equation of motion, techniquesin e th di vergenc of the magn tic field is zero, the magnetic field lines system canwhich are free from secular) or dissipative errors. Symplectic integration schemes have b n

be treated a a Hamiltonian system. For Hamiltonian syst ms the phase space structuredeveloped in the past few years and seem to be the best way to numerically investigate theis subdivided into chaotic and regular regions, and this subdivision has a correspondence3----- --"....

~

I

properties of Hamiltonian systems.

5

The research approaching the toroidal heli cal magnetic field problem from the Hamil-

6

CHAPTER 1. fNTRODUCTION

1.2

Overview

tonian point of view ends with works done more than ten years ago (refs. [1], [2] and[3]), partly because of the strong push towards research on tokamaks. However, we have
The remaining part of this thesis is subdivided into fiv chapt rs and an appendix.
two good reasons to revive the Hamiltonian approach to the toroidal helical magnetic field
In chapter 2 we first review those concepts from Hamiltonian mechanics which will be

problem:]). The Japanese Ministry of Education is pushing the research concerning the helicalsystems approach to magnetic confinement and the LHD (Large Helical Device) is in phaserelevant to our work and describe the Hamiltonian nature of the magnetic field lines system.
Then the concept of Lie transform is introduced and its utilization for a p rturbation th ory
that is more efficient than the classical one is described.

of construction.
2). The last ten years have seen an impetuous development of chaos physics, and time

In chapter 3 we present th calculations that allow us 1,0 constru t a Hamiltonian for the'is ripe for research about the chaotic properties of the toroidal helical magnetic field lines

magnetic field lines system. Th Hamiltonian is obtained for Lh cylilldrical approximation

system.

of the field and for the toroidal case.

Recently T. IIatori and T. Watanabe (ref. [4]) have obtained an explicit form for the
In chapter 4 we will discuss in detail symplectic integration schemes and the development
Booz r's magnetic coordinates in the first order toroidal correction to the cylindrical helicalof a new linear symplectic integration scheme which will be used to solve the IIamilton's
magn tic confignration. In this work a more systematic procedure (Lie perturbation exequations derived in chapter 3.
pansion technique) is used to proceed to higher orders and a new linear symplectic method
is onstructed and used to integrate numerically the equations of motion.

In chapter 5 some numerical results and their discussion is presented. Of particularinterest is that our Hamiltonian, yet of simple form, is abl to reproduce fully the regularand the chaotic region of the field.

In chapter 6 the concluding remarks are given.

In the appendix the coefficients that characterize our method are explicitly calculated.

---------

-

....

7

1.2.1

Notation

Vie have tried to avoid a rigorous, but heavy notation , in the the hope that the context
will help the reader to decipher the notation. We do not use boldface character for vectors,
so a symbol like P could indicate a single variable or the vector (PI,' .. ,Pn)' The Einsteinsummation convention is used extensively. When the indexes are Latin the sum is from 1,
when the indexes are Greek the sum is from O. For the partial derivative both the symbols:x and ax are used.

-----

--

-

8

CHAPTER 1. INTRODUCTION

10CHAPTER 2. THE LfE TRANSFORM

2.2

Hamiltonian Mechanics

In classical mechanics a syst m ofparticles is describ d by the NC'vvt.on equations(2.1 )
where i = 1, ... , N and Pi is the total force acting on th ith point of mass m i . If con traint

Chapter 2that limit the motion of the system are present (for xample the motion may be constrained
to occur on some surface or curve), then the coordinates r}, . .. ,rN are not ind pend nt

The Lie Transformof each other. However, if the constraints are holonomic, that is if they are expr ssibl byequations of the type
j=l, ... ,I<,

then it is possible to change the coordinates from r I l " "

2.1

(2.2)

rN to a new set of independent

Introduction
coordinates qI, ... , qn, with n = 3N - I{, which are called generalized coordinates. For
conservative systems, that is for systems such that the forces are derivable from the gradient

In (,his chapLer we will first reVlew those aspects of Hamiltonian mechanics which willof a potential function V(rl,"" rN), then a Lagrangian L( qIl .,. ,qn) = T - V, where T isb r 1 vant t.o our work. In particular the notion of canonical transformation and thethe total kinetic energy of the system, can be constructed, 1 and the equations of motionsym} lectic nature of the flow of Hamiltonian systems will be described. Then the equation

for the generalized coordinates are gi ven by

for the magnetic field lines will be put into a form which is equivalent to Hamilton 's
equation. In thfinal section the concept of Lie transform will be introduced. The Lietransform can b

us d a a basis for a perturbation theory which is more efficient andd 8
8--L--L=O.

dt 8qj8qj

1 vVe observe that a Lagrangian can be constructed also for the electromagnetic forces on moving charges ,though in that case the forces are derivable from a velocity-dependent potential.

1 gant than the clas ical one.
9

(2.3)

11

Once a Lagrangian has been constructed for a given system, we can construct the Hamil-

12

CHAPTER 2. THE LIE TRANSFORM

in the new variables Q, P be in the same form as (2.5), that is ther must exist a function

K( Q, P, t) such that Qi = 8~J{ , Pi = - 8~i K, which is in general not true for arbitrary

Lonian fun ction II , defined as

H (q,p, t) =L qiPi - L(q, q, t),(2.4)

coordinate transformation. A transformation which preserves the form of Hamilton's quations is called a canonical transformation. V-le now state the onditions under which thewhere P =

PI"'"

Pn ar

called generalized momenta and are defined by Pi =8~i L.

tr ansformation
By this procedure the n coordinates ql , "" qn are substituted by the 2n coordinates(2.7)
q}, ... , qn, Pl, ... ,Pn and the n, second order Lagrange equations by the 2n, first orderis canonical. If M is the Jacobian of th transformation (2.7)
II amilton equations

.

8
Pi

qi = -8 H'.

8 H

Pi = --8
qi

.

M =

(2.5)

( ;qQ ;q p),
8 Q

8p

Th IIamilton equations can be derived also from the variational principle

8 P

8p

and J the skew-symmetrical matrix

5

Jdt(L-,Piqi - H(q,p, t)) O.
=

(2.6)
(2.8)

This variational principle will allow us to compare the equations for the magnetic field
line now with the I1amilton equations and to identify the magnetic field line system as athen a necessary and sufficient condition for the transformation (2.7) to be canonical is

Hamiltonian system.

that the condition Mt J M = J is satisfied (ref. [5], pag. 391-397). This condition is called

Sol ving analyLi ally the Hamilton equations is in general not possible except in very

symplectic condition.

jJccial, t.hough v ry important , cases. Often it is useful, even if we want to proceed to

While when changing coordinates it is often convenient to change coordinates by mans

a numeri al inL gration of eqs. (2.5), to make a transformation of variables in order toof a canonical transformation , this is not strictly necessary, for som systems are bettersimplify th

understood using noncanonical variables. However canonical transformations manifest their

quaLion or in order to reveal some symmetry of the physical system.

importance when we consider the time evolution of the system. In fact , let us suppose thatto Q = (Ql Q2, ... Q ), P = (PI P2 , ... , PN) , it is desirable that the equations of motionq = q( qo , Po , t) and p = p( qo , Po, t) are a solution of Hamilton 's equations. Then it can be13

shown that the coordinate transformation (qo, Po) - t (q: p) is a canonical transformation,that is the symplectic condition Mt J M = J holds, where

14

CHAPTER 2. THE LIE TRANSFORM

We conclude this section by putting the variational principl (2.6) in a form which will
be useful later. Let us introduce in the phase space the new variables zi, i = 1, ... , 2n, with

i = qi and zi+n = Pi, for i = 1, ... , n. We also d fin the quantities Ii = Pi, IHn = 0,i = 1, ... ,n. Then, with these notations, (2.6) becomes

The time evolution of a Hamiltonian system satisfies therefore the sympledic condition, and(2.10)

this fact has important practical consequences, as we will see when we discuss symplectic

If we also set zO = t, 10 = -H, then 2.10 takes the very compact form

integration m thods.
An important property of Hamiltonian systems is the conservation of phase space volume ,(2.11)

a fact known as Liouville th orem. Liouville theorem is easily proved: if we define a veloci ty
vector v = (q, p), then the di vergence of this vector iswhere the dot denotes the derivative with respect to time. Hamilton's equations are foundby carrying out the variation in (2.11). Since the int grand of (2.11) is a scalar, if we\1. vchange coordinates,(2.9)
t.hat is, if we consider the motion of the points of the phase space as composing a fluid, thisnuid i incomprcssi ble. Actually the conservation of phase space volume is a particular caseof a more general th orem, which says that the flow of a Hamiltonian system preserves thediff l'ential form w = dpl 1\ dq l + ... + dpn 1\ dqn and the exterior powers w\ ... ,wn . 2 Thisthcorem i equivalent to the fad that th time evolution of a Hamiltonian system satisfies
thc sympl di . condition. The Liouville theorem corresponds to the conservation of w n .
2H i al 0 ajd that the flmv of a Hamiltonian system is symplectic.

(2.12)
then the new equations of motion are again obtained from (2.11). If the coordinate transformation is canonical, then the equations of motion in the new variables will be in Hamiltonian form. If the coordinate transformation is not canonical, the equations of motion
obtained from (2.11) will not be in Hamiltonian form. Under the coordinate transformation(2.12) the I-form lp,dzP, will become fp,dZP, and in general the compon nts fn+i' i = 1, ... ,nwill not be zero. However, if the transformation is canonical then we will have r n+i = 0,i = 1, ... , n. This can also be taken as a definition of canonical transformation: a canonical15transformation is a transformation (2.12) such that r n+i = 0, i = 1, ... , n .

CHAPTER 2. THE LIE TRANSFORM

16which has to be compared with (2.10). We see that -A3 plays the role of the Hamiltonian

Let us now see how the equations for the magnetic field lines flow can be put into

for the magnetic field lines system and that Zl and z2 are canonical1y conjugate variables.

Hamiltonian form . Writing the magnetic field as B = 'V 1\ A, then [or the magnetic field

The variable Z3 is the analogous of time. The variational principl (2.14) can also be writtenlines we have the equations

asdx('V 1\ A) 1\ dA = 0,(2.13)

where A is an arbi trary parameter. If now we consider the variational principle (ref. [11])(2.14)(2.19)

which has the same form as (2.11).

The field line flow is then obtained by applying

Hamilton 's equations to - A 3 .
The vector 1M will be called I-form (more precisely the 1-form i th quantity Ip.dz M ).

then from the calculus of variations we know that we must have
The addition of the derivative of an arbitrary function S to the] -form does not change th

d 8g
8g _ 0dA 8Xk - 8Xk - ,(2.15 )

where 9 = Ap.(x)dxM/dA and the dot denotes the derivative with respect to A. We thus
obtain(2.16)

which is t he sam as (2.13). We now perform a gauge transformation A~ = Ap. + 8p,S inorder to eliminate one of the component of the vector potential, say the component A~.
lhcn , we obtain (omitting the prime)(2.17)and introducing th new variable p = A 1(z\Z2,z3) to replace Z2, we arrive at(2.18)

equations of motion obtain d by carrying out the variation of (2.11), that is the l-forms
IP, and IP, + 8p,S are equivalent.

17

2.3

The Lie Transform

18

CHAPTER 2. THE LIE TRANSFORM

based on the so called Lie transform , which is a sp cial kind of coordinate transformation.
The perturbation theory based on the Lie transform has a number of advantages. It is

In classical perturbation theory for Hamiltonian systems, in general a generating function
canonically invariant and it is possible to give a direct expansion of any function of the old
F is used to perform a canonical transformation from the old variables q, p to the newvariables in terms of the new variables. We present in this section a brief account of th 'vari abIes Q, P. This generating function depends on both the old and the new variables,
perturbation theory based on the Lie transform. 1 he great majority of the results are notand as a consequence the transformation itself and the relation between the old Hamiltonian
derived, and the reader is referred especially to ref. [11] for a rigorous pres ntation of th

II and the new one K appear in a mixed form. For example, if the generating function istheory.
of Lhe type F( q, P, l), then [or the coordinates we have the relations
Let us consider the coordinate transformationoQ(q, P, t) = opF,op(q, P, t) = oq F,(2.20)

and for the IIamiltonians

from the variables Zl, . .. ,zn, to the variables Z\ ... ,zn, wher c is a small pararn ,t '1' andK(Q, P, t) = H(q, p, t) +\¥s~

(2.22)

ato F.

(2.21 )

that (2.20) is in a mixed form, while what we actually need are the relations

Q(q, p, t) and P(q,p, t), and that (2.21) is a relation between the functions Hand I{ atZr, Zf, ... are functions of n variables. W will use the compact notation(2.23)The transformation (2.23) is called a Lie transformation when the differential equationcorr sponding points in phase space. When (2.20) and (2.21) are used as the basis fora p rturbation analysis the result are very lengthy formulas even for very low orders. It
would be prcf ra.ble to hav a formalism that yields directly the new variables in terms

together with the initial condition ZI-i( z, 0) = zl-i, is satisfied by the functions ZI-i. Theof th old one or vic -vel' a, and such that th relation between Hand I{ is a relation

vector gl-i is called the generator of the Lie transform, and compiet ly charact rizes it.

bet vve n fun tion .

Let us now see how coordinates and scalars transform under a Lie transform. We first

This formali m has been developed by researchers working in celestial mechanics ([6],introduce the operator Lg which maps scalar functions into scalar functions according to[7], [] [9]), and ha b en pres nted in the review articles [10] and [11]. It is a formalism(2.24)

19

Notice that this is a functional relationship, so the symbol for the independent variables is20

CHAPTER 2. THE LIE TRANSFORM

expression of the form

a dummy, and it could be z/1- or Z/1- or any other. It is then found that the old and the newcoordinates are related by(2.28)
where r 1 is obtained by a Lie transform characterized by the generator 91(2.25 )

r 2 l' obtainedby a Lie transform characterized by the generator 92 and so on. It is found that, and r
are related bywhile a scalar s( z) transforms in the opposite way, that isr = T, + dS,(2.29)(2.26)
where
W will be particularly interested in the transformation properties of the expansion, according to a smallness parameter c:, of a I-form

with, =, + c:, + c:, +"',o1(2.30)

2 2(2.31 )
(2.27)and dS is a total differential which represents a gauge transformation or the l-form andwhich could represent, for instance, the magnetic vector potential or some other physicaltherefore does not affect the equations of motion. Expanding in powers of c: we obtain

quantity. Th - uperscript in the ,'s indicate the order in L The I-form (2.27) is better
treat d by a composition of individual Lie transforms, each of them characterized by its(2.32)

gen rator g/1-, In the following the lower index will represent the perturbative order, sothat we will write g~ to indicate the generators of the nth Lie transform, of order c:n , and(1- c:L 1 + ~E2Li + ... ), + dS.
Collecting the terms and using also the expansions for, and r we obtain

Ln will nmv indicat the operator such that, if (7 is an arbitrary I-form, then (Ln(7)v =g~(fjlt(7l1-fjv(7/1-)' 3 Jotice that the convention is the opposite for the I-form: for the I-form
th superscript r presents the pcrturbative order. Our aim is to transform (2.27) into an
3We have horten the notation by writing Ln instead of L g ",2- L 2, +,

dS

0

2

-L

1, + 21 L2l'0 ,1(2.33)

21and so on. Here S represents a gauge function that will not alter the equations of motion.

CHAPTER 2. THE LIE TRANSFORM

22

The general form of r is therefore

The coordinate transformation associated with the transformation of the I-form is found(2.41)

observing that the coordinat.e transform controvariantly, that ist, and the r suIts of the preceding lower orderwhere Gn is a I-form calculated from ,.

(2.34)
calculations. The generator is contained in the term Ln/ D• The expressions of G1 , G2 andwhere1
T - 1 -- T-1r.-1r.1
2
3 "',(2.35)(2.4 2)with(2.36)
It is possible to choose g~ = 0 to all orders, and this corresponds to th fad that usually we
Again, expanding and collecting the terms, we obtaindo not want to transform also the time when we change coordinat s. The 2N compon nts(2.37)
To simplify the notation it is convenient to introduce the tensor w, called Lagrange tensor,

g~ and the scalar Sn can be chosen as to bring the 2N+ 1 compon nts of r n , where N is

the number of degrees of freedom, into some desired form. We want a form in which on ly
the temporal component of r , which will be the <p-component in our case, is not zero, andand defined by (~ is a I-form)(2.38)

this can be done choosing(2.43)
so that w can write(2.39)

where J~j is the inverse of the spatial part of th Lagrange tensor. With this choice of thegenerators the temporal component of the I-form becomes\\Then (J i one of the t rms of (2.27), then we will further simplify the notation by writing

w~v inst ad of w~:, so that, [or exan1ple(2.44)(2.40)

where vt is the Poisson vector, defined as vd = J~jwgj' Voo = l. We stress that with the

23choice (2.43) of the generators, the temporal component of the 1- form is the only one
which survives after the transformation.

24

CHAPTER 2. THE LfE TRANSFORNI

26

3.2

CHAPTER 3. NIAGNETIC FIELD LINES HAMILTONIAN

Magnetic Potential

The general scalar potential satisfying the condi tion \7 2 \]"10, regular at ~0, andallowing for an axial current is (refs. [12], [13])+00

Chapter 3\[J = Jcp

+ (1- ~COS7])1/2 L

Qlm Ulm (e) ei l '1+im 'P,(3.1 )

I,m=- oo

where

(3.2)

Magnetic Field Lines Hamiltonian

The coordinates ~, 7] and cp are toroidal coordinat s, in terms of whiclt Lh)

arLesiancoordinates are z = ~sin1]/(l- ~COS7]), x = (1- e)1/2coscp/(1 - ~COSTl), and y = (L e)1/2 sincp/(1- ~COS7]).

3. 1

Units are selected such that J = 1. Thco fficienLs Qlm arearbitrary provided that Qlm = Q~(lm)' The [unction Q~ is the modified Legendre function

Intro d uctio nof second kind.

From the potential \[J the component for the vector potenti al for thedivergence-free and curl-free magnetic field are obtained as explained in ref. [13] and the
In t.his chapter we discuss Hamiltonian formulations for the magnetic field lines system. We

vector potential can be written as

Ilr t u 'e a gaug transforn1ation to eliminate the component 'Ij; of the vector potential andA=then expa.nd t.be non vanishing components in Taylor series. Then we apply perturbation
theory based on the Lie transform to transform syst matically to canonical variables. In(2: QlmA~m, A~ + L QlmA~;\ 2: QIOA~),
I,m

(3.3)

I

where

section 3..1 we how the procedure rather in detail for the cylindrical limit approximation of

A~=(1/2)[(1-cosTl)-lln(
the field. InI,me

1 - ~ cos 7]_
1 - cos 7]
1 - ~ )+(1-COS7]) Iln( l+~ )],(3.4)

ction 3.5 \'I'e d rive a Hamiltonian correct until the third perturbative orderfor the toroi dal h Ii cal ca

and in the following section we discuss a different Hamiltonian.

25(3 ..5)

27

Aim =
'7

'/,

m(l - ecos 17)e(1- e)8c(1 - eCOS7'l)1/2Ui (t)e i1'7+imcp
<,

'/

m

~

,

A~ = - foedx[x(l- XCOS7])t 1U/O(x)8'7(l- XCOs7])1/2e ii7J .

(3.6)

28

CHAPTER 3. lvIAGNETIC FIELD LINES HAMILTONIAN

The I-form is of the type i = iodO + 'Ycpdcp, wh r

(3.7)

(3.13)We select the ha.rmonics (I,m) and (/,0), that is we will consider the potential
and

(3.8)
1'1'
where E = aim, E' = a/O. Now, in order to apply the Lie transform, we consider the I-form

J

E'A~ - (m/I)(A~ + EA~m) - E d~(8cpA~m - (m/I)87JA~n)J= t'A~ - (m/l)'Yo - t d~acpA~m.

(3.14 )In the next section we will expand the I-form in Taylor s ries in Lh variable ~ .

wh re S' is a ga.uge function. Exploiting the gauge freedom, we choose the gauge functionin order to make the ~ component of the I-form to vanish, that is we choose the gauge

3.3

fundion to be(3.10)and ther for the I-form becomes

Taylor Expansion of the Potential

To apply the Lie transformation we need to have the I-form written as an expansion of thetype (2.27). We will consider E, E' and ~ as smallness parameters of the same order, and(3.11)\1\1 now introduce the new variable 0 = 7] + (m/l)cp, where l is the poloidal multipolarity

we will expand the I-form in the variable ~ about the point e = O. W proce d first to the
expansion of the quantities U1m(O. The functions Ulm(O admit the integral repr s ntation(ref. [13])and m is the number of field periods, so that we obtain[A~ + EA~m - t(3.15)

Jde8'7A~mJdO
Expanding (1 - e)-m/2 and (1 + e cos t)m-l-l/2) we can put (3.15) in the form+[E'A~ - (m/I)(A~ + EA;)m)-E

Jd~(8cpA~m - (m/l)8'7A~m)]dcp.

(3.12)(3.16)

29

where the normalization constants 131m have been incorporated into the constant quantities

CHAPTER 3. l\1AGNETIC FIELD LINES HAMILTONIAN

30(3.17) we obtain, up to the third order, the following expansion for If)

U~m' Uim' .... From now on we will set l = 2. For the other terms we obtain

I~

=

1j;,

IJ = (1 /3)(21j; )3/2 cos( (mrp) /2 - 0),(3.18)I~

A 1]2m= (2i/m)U~mee2if) +and for Ir.p(i/m)(2Uim 2 - U~m)~4e2if) + "',(3.17)

A2m e(2/m)U~mee2if)+,~= -(m1j;)/2,

I~= -(1/3)y2m1j;3/2 cos((m<p)/2 - 0)-(l /2m )Ufm (2cos 7] - sin 7] )ee 2i f) + ... ,= -iUfoee 2i(1/5)y21j;5/2(7 cos((mrp)/2 - 0) + cos((3m<p)/2 - 30)+
(1/(4m))c1j;2 (35 U~m - 48 Uim) sin(20),(i/2m )cos 7]U~m(3e2if) +

A20
r.p

=

2 c'1jJ U~o sin(m <p - 20) + 2 c 1j; U~m sin(20),

1]-

I;

(1/6)U~o(sin7] + 2icOS7])ee 2i I7 + ....

=-(3/4)rn1jJ2 - (1/4)m1jJ2 cos(m <p - 20)(3y2)-lc1jJ3/2U~m sin((mrp)/2 - 30) - (3y2)-lc'1jJ3/2Ufo sin((3m<p)/2 - 30)-

(1/ V2) c'1jJ3/2U~O sin( (mrp )/2 - 0) + (1/ V2)c1jJ3/2U~m sin( (m<p )/2 + 0),I~-(I/(5y!2))m1jJ5/2 cos((3m<p)/2 - 30) - ((7m)/(5V2))1jJ5/2 cos((m<p)/2 - 0)+(3/8) ell} U~m sin( m <p) - (3/8) c'1jJ2 U~a sin( m <p - 20)The quanLiLi

lm are constant, and Lhe only values that we will need are Ufo = 1/2,

Ufm = (m 2)/ , Uio = 0, Ulm = 0 and Uia = 35/96. Putting all the terms together,
introducing Lb n w variable 1jJ = (1/2)e and taking the real part of the various terms in

2 c'1j;2 Uia sin(m <p - 20) - 4 c 1jJ2 Ufm sin(2 0) + 8 c 1jJ2 Uim sin(20).
(3.19)

31

3.4

Cylindrical Limit Approximation

32

CHAPTER 3. MAGNETIC FIELD LINES HAMILTONIAN

and, for I = 2, we have for the first two terms
1(0) =In this section we apply the Lie transform technique to the cylindrical limit approximation

2

~(m)2
2! 2

'

of the magnetic field, which is the approximation that neglects the toroidicity of the field.
The field becomes therefore a straight helical field. The calculabons will be carried out(3.26)(3.27)For the cylindrical limit case it is convenient to u e the formulain great detail in order to illustrate with a simple example the use of the Lie transform.

[A~ + tA~m - t

l h calculation scheme described below will be used, unchanged, when we will take intoJd~8!]A~m]dO+[t'A~ - (m/l)(A~ + tA~m)

i1.ccount also the toroidicity of the field. In the case of the cylindrical limit approximation-t

we haveJd~(8t,OA~m - (m/1)8!]Akm)]dcp,(3.28)and observe that 8t,OA~m - (m/ l)8'11At = O. Then, expanding th quanLit,jes A~, A~m, A~me

T A '11-2"'(3.20)(3.21 )
(3.22)(3.23)in Taylor series about the point ~ = 0 , taking the real parts and introducing the variables'I/J = (1/2)~2 and () = 7] + (m/l)c.p, we obLain

I~ = 1jJ,

I~ = -(m/2)'I/J,
I~ = c1jJ(m 2 /4) sin(2()),(3.29)'I'll ft( mO are Lhe modi ned Bessel functions, that is solutions of the equation
From the expansion above we see that the lowest order part of the vector potential is(3 .24)(3.30)and they admi L Lhe s ries expansion
This expression is already in canonical form, () and cp play the role of canonical conjugate2
k== ( T )2k( !?!.£.
2 ) L:k=O k!f(k+l+l)ft(m~)t2(j(O)
I

~

-

~

---

+ t2j2(I) + t412(I) + •.• ,~~(3.25)

variables and the unperturbed Hamiltonian is H = (m/2)1jJ. Actually it was the introduction of the variable 'I/J = (l/2)e that put the lowest order potential in canonical form. The

33

lowest order of the vector potential, and the corresponding Hamiltonian H = (m/2)'lj; will

34

CHAPTER 3. MAGVETIC FIELD LINES HAMILTONIAN

The Lagrange tensor is thcrefore the 3-dimensional matrix

bc considered the unperturbed part of the potential.
We now recall, for convenience, the formulas for the generators of the Lie transform andfor thc temporal part of the transformed I-form. We have (see eqs. (2.43) and (2.44))(3.31 )

o m/2

o

o -1

- m/2 10

and its spatial part is(~ ~l 1

(3.32)
whcre J~j is thc inverse of the spatial part of the Lagrange tensor and vt is the Poissonofrom which we obtain its inverse, the matrix J~jv cLor, defmcd as vd = J~jwgj' Va° = 1. In our calculations we have three variables: !.p,o and 'cPo Th variable!.p corresponds to the time , and therefore the temporal part of the
I-form is tbe one corr sponding to !.p. Let us now see how eqs. (3.31) and (3 .32) appear in
The Poisson vedor is now obtained from the formula

vd = J~jwgj ' Accordingly

our case . First we notc that the zeroth component of the transformed I-form is equal to
Yo - 0,th unperturb d one, that is

1 71/.1 -r O = 'lj;dB - (m/2) 'lj; d!.p.

(3.33)

VVe use this expr ssion to calculate the Lagrange tensor wo, the inverse of its spatial part

v,60
- m/2
,

and we remind the reader that by definition Vact' = 1. Using these results, from eq. (3.31)
we obtain

J~j and th Poi sson vector vt. For the Lagrange tensor we obtain(3.38)(3.39)(3 .34)(3.35 )o =w01/.l

-

- --

-~--

--

a0,,,,, - a1/.116 - - 1•00_(3.37)(3.36)and g'!t = 0 for every n. The expression (3.32) for r~ is given by(3.40)

35

36

CHAPTER 3. MAG1VETIC FIELD LINES HAMILTONIAN

We now proceed to the calculation of the first. order contribution to ro. We start byand substituting the Poisson vector, we obtain(3.41 )

taking r~ as the average over () of the quantity VtC1~ = ,~ + (m/2),J, which vanishes,smce

VVe see that, for every order n, the basic formulas (3.38), (3.39) and (3.41) suggest two(3.43)

possible strategies:
We have now to calculate the first order gauge function and the first order Lie generators,
1) first we choose the g nerators. Then (3.38) and (3.39) give us the equations by which
which will be necessary in order to calculate the second order component of the transformed
we determine the gauge function S. Substituting in (3.41) we obtain the transformed
I-form. To calculate the first order gauge function, we have to solve the equation
L-form.
( .44)

2) first we choose the gauge function S, then from (3.38), (3.39) and (3.41) we obtain
A solution of this equation is a gauge function which does not depend on 'P, that isth generators and the transformed l-form.
The two strategies are completely equivalent, and using one or the other is a matter of

51 = (1/4)t'l/;cos20.

(3A5)

choice. We will proceed according to 2), because choosing first the gauge function allows
The generators are obtained from eqs. (3.38) and (3.39) which explicitly readon to specify more directly the form of the transformed I-form. A typical choice would
he to choose 5 n such that the equation

gf = -81/151 = -(1/4)tmcos 20,(3.46)

gt = 8eS1 + ,J = -(1/2)tm~sin 2().

(3.4 7)(3.42)At the first order we see therefore that there is no contribution to the unperturbed Hamilwh re the brackets denote average, is satisfied. This corresponds to taking rn as the average

tonian. Let us proceed to the second order calculations. From (3.41) w obtainof entp + (m/2)CnB over the sam variabl s as (3.42).
(3.48)
Finally we tr ss that (3.3 ), (3.39) and (3.41) are derived under the requirement that
only th t 111poral part of r that is rtp is not zero after the Lie transformation of coordinates.

and therefore we now have to evaluate the 0 and 'P components of the quanti ty(3.49)

37

38

CHAPTER 3. MAGNETIC FIELD LINES HAMILTONIAN

We have(3.55)(3.50)(3.56)
Substituting these relations into (3.52) and (3.53) we obtain(3.57)=

g1/11 a 1/1 "VI
I c.p+ gl(Ja(Jlc.p1(3.51 )

(3.58)

Beside, we have
Using these results we obtain(3.59)(3.52)(3 .60)(3.53)

Averaging over e the quantity(3.61 )and ther [ore w hav to calculate first the components of the operator L1/°. We obtainwe obtain the second order contribution to the transformed I-form(3.54)-------(3.62)

39

Accordingly, the expression of the transformed I-form up to the second order is40

3.5(3.63)

Tn this expression the magnetic potential is expressed in canonical form, with 'ljJ and ()playing the role of action angle variables. The Hamiltonian up to the second order is

CHAPTER 3. MAGNETIC FIELD LINES HAMILTONIAN

Helical Toroidal Potential

We now proceed to the calculation of the transformed I-form in the toroidal ca ,along thesame line followed for the cylindrical limit approximatlon. The starting pOlnt is again the
expansion of the magnetic potential, which we r port for conv niencc until the first order(3.64)(3.65)it(l/3)(2'ljJ?/2 cos((mcp)/2 - ()) ,The equations of the magnetic field lines flow are obtained directly from the above expression. This re ult is consistent with the result obtained in ref. [2].

-(m'ljJ)/2,i~= -(1/3)V2m'ljJ3/2 cos ((mcp)/2 - 0) -(3.66)

2 t ' 'ljJ U~o sin(m cp - 20) + 2 t 'ljJ U~m sin(20).
We note that the zeroth order term is equal to the zeroth order term of the cy lindrical limit
approximation. The zeroth order in the expansion is the unperturbed part of the potential

,

o

m

= 'ljJ d8 - - 'ljJdcp ,2(3.67)and it is in canonical form , with 8 and 'ljJ playing th rol of action angle variables and -i~
playing the role of the Hamilt onian for the unperturbed syst m. We now proceed to t hecalculations of the transformed I-form to the various perturbative orders. For the zerothorder we have that the I-form is left unchanged, that isr O = ,0 = 'ljJ d8 - (m/2)'ljJdcp.

-(3.68)

41\Ve also note that this is the same expression that we obtained in the cylindrical limit
approximation , therefore the Lagrange tensor w O, the Poisson vector VIJo and the tensor J~j42

CHAPTER 3. MAGNETIC FIELD LINES HAMILTONIAN

Since Sl = 0 and C1 1/J = 0 \ive obtain

gf = 0,(3.76)

are the same as for the cylindrical limit approximation of the field, that is we have

w~o = 0, w~1/J = m/2, w~1/J = -1,(3.69)(3.77)

andThe equation for the second order contribution is(3.70)(3.78)
This means that eqs. (3.41) , (3.38) and (3.39) are left unchanged, and therefore for the
first order transformed I-form we have to calculate the expression

which becomes, after substituting the values of the Poisson vector(3.71 )

(3.79)In the above equation we can choose the gauge function to be zero, so that the expression
For the second order we need therefore to evaluate the quantity C 2 = ,2 - L 1[or r ~ becomes,1 +(l/2)Li,o.

The calculations of the explicit expression of the operators L's are straightforward but ther~ = ,~ + (m/2),~,(3.72)
algebra is quite lengthy, so we give only the results

and, using th T'aylor expansion for, we obtainr~ = -21f;c'Ugo sin(mcp - 2B) + EU~m sin(2B) .

_(2 1 / 2/3)7jJ3/2 cos(m/2cp - B) X
(3.73)(2 1 / 2m1f;1/2 cos( m/2cp - B) + 4tU~o sin( mcp - 20) l·or the fir t order generators we have
4t'Ugm sin(20)),(3.80)(3.74)(3 .75)-(3.81 )

43(1/2)( Li'/) <p= - (2/3)m'lj} cos( (m/2)r.p - fJ)2,

CHAPTER 3. MAGNETIC FIELD LINES HAMILTONIAN

44(3.89)(3.82)

To proceed to the calculations of the third order we have to calculate now the quantity(3.90)(4/3)m1j} cos( (m/2)r.p - fJ)2.

(3 .83)Calculating the above expression and choosing again the gauge [unction to vanish, we find
The equation for r~ isfor r 3 the following expression(3.84)(3.91 )
and, choosing again the second order gauge function to vanish, we obtain

2E'lj;2Uim sin(2fJ).
Up to the third order the Hamiltonian is given by

JI = - r <p = - r~ - r~ - r~ - r!.

(3.85)For convenience we collect below the various contributions.
JIo =For th second order g nerators we obtain from eq. (2.43)(3 .86)(3.92)

m'lj;(3.93)

2 '2 E' 'Ij; ugo sin( m r.p - 2 fJ) - 2 E'Ij; ugm sin(2 fJ),E

'Ij;~ U~m sin(¥ - 30)

J2

E' 'Ij;~(3.94)

ugo sin(~ - 3 fJ)
Vi(3.9,5)(3.87)

Since"2 = 0 and C ,p = 0 we obtain2

g~ = 0,(3.88)

45

3.6

Integrable Model

46

CHAPTER 3. NIAGNETIC FIELD LINES HAMILTONIAN

is satisfied. We then obtain(3.98)In this section we discuss a dilTerent Hamiltonian for the toroidal helical magnetic field. We
and, consequently, introducing the variable 7] = () - (m/2)<p, we obtainobserve that the perturbation theory based on the Lie transform allows a certain freedom ofr~ = 2c'U~o1j;sin 27].

choice in selecting the transformed variables. In fact it can be noticed from eq. (3.41) thatchoosing a different gauge function would lead to different Hamiltonians. However there(3.99)For the first order generators we obtain
()

0

2t

0

91 = -U2m -cos 2 ,mis no contradiction in this fact, since the different Hamiltonians are expressed in terms of(3.100)differenL variables. In the same sense, when we perform a canonical transformation weare free to s I ct arbitrarily the generating function, 1 leading to different Hamiltonians
expressed in tern1S of different canonical variables. For the same physical system, dependingon th choice of the canonical variables, we can have completely different Hamiltonians,

cP _

0

91 - -U2m

4c1j; .
3/2 cos 7]--;:;:;-Slll
20 + (21j;)3'

For the second order we need therefore Lo evaluate Lhe quantity C 2 = ,2 -(3.10 I)
fJl, 1+ (1 /2)Li,o.

After some lengthy algebra we arrive at the result 3

one being tim -dependent and one being time-independent, or one being the energy of the(3.102)syst m and one being not, and so on. 2. All these different Hamiltonians describe the same
Therefore, up to the second order, we have the following expression for the Hami Itonian

pbysical syst m.

- H =Let 11S consider now eq. (3.41) for n = 1 and let us choose the gauge function such thatthe quationr = r~ + r~ + r~,
cP(3.103)that is(3.97)If

ITh choice of the gauge function is, however, subject LO the restriction that it must be bounded in the(3.101)variabl <po
2 For a nice di cu

Ion ee [5], sec. 8.2.

3We develop these calculations in ref. [14].

47

It can be seen that this Hamiltonian is time-independent, and therefore is integrable andcannot exhibit chaotic behavior. For another discussion on this problem, see ref. [1 5], pag.

130

3.7

Conclusion'vVe have derived in this chapter an expression for t he Hamiltonian for the toroidal helicaln1agnctic field li nes system up to the second perturbative order for the cylindrical limit
approximation of t he field and up t o t he third perturbative order in the toroidal case.
Going to higher or ders will be a straightforward application of the same procedure used inthi s work, and t he only difficulty to be expected is algebraic complication. The parti cular lI arn ilton ian we have deri ved is not the only possible choice, since we can manipulate
the Lie transformation in order to get different Hamiltonians, and we have discussed thispoint in section 3.6 , where a different Hamiltonian has been given up to the second order.
Di fferent expressions for the Hamiltonian could be of interest when analyzing particular
problems, or [or Lh appli cation of symplectic integration schemes to the solution of HamilLon' equations.

48

CHAPTER 3. lvIAGNETIC FIELD LINES IIAlvIILT ONIAN

50

CHAPTER 4. SYMPLECTIC INTEGRATION

and it is natural to look for numerical methods that guarantee that apart from roundoff errors, the conserved quantities of a giv n Hamiltonian syst mare conservcd also in
the numerical integration process. Sincc the differential form w is an invariant for all
Hamiltonian systems, in the recent years much attention has b n given to numerical

Chapter 4schemes that conserve w (see refs. [16], [17], [18], [19], [20], [21], [22]). These num rical
schemes are usually called symplectic schemes, b caus the cons rvation of w is quivalclIt
to the request that the transformation (qn, Pn) -4 (qn+l, Pn+l) b tw en two uccessive tirne-

Symplectic Integration

steps is canonical (canonical transformations are also called symplectic transformations).

It should be noted here that for non integrabl Hamiltonian systems it i. not possible
to have symplectic algorithms that conserve also th energy. In fact, if II is conserved by

4.1

Introduction

the numerical flow and if the original system has no other indep nclent int.egrals then for a
symplectic algorithm the approximate and real flows are the same up to a r parametrization

Ilamiltonian systems possess in general many invariant, or conserved, quantities, which,of time, and we would thus have solved the equations of motion (ref. [23]). 1
depending on the particular system under consideration, may be the energy, t he angular
However, even if the energy is not conserved by the sympl ctic method, there is no
momentum and so on.

Another invariant which is common to all Hamiltoni an systemsecular error in the energy, that is the error on the energy do s not increase steadily wi th

i the phase space volume, which is a consequence of the conservation of differential formw = dpll\. dql+... +dpn I\. dqn and the exterior powers w 2 , ••• ,w n by the flow of Hamiltonian

time, but is bounded, as a consequence of the fact that symplectic numerical schemes areby themselves Hamiltonian systems.

systems. Notice that th conservation of w prevents Hamiltonian systems from exhibiting
The symplectic integration of Hamiltonian systems is by now an established technique,di ipativc dynamics.

IThis remark does not apply to linear Hamiltonian systems and to one-degree of freedom, Lime-

Since only very seldom it i po sible to find an analytical solution for the Hami lton

independent Hamiltonians, which are integrable. For these systems it is Lherefore possible to have sym-

quations, (h - uplectic algorithms that conserve also the energy.

of numerical methods to integrate the equations of motion is common,
49

--

-

......-......o...-..oo.....~~ .... _ _ _ _ _ _

~

_~_

~

... ~~L....-

, ..

51

52

CHAPTER 4. SYMPLECTIC INTEGRATION

however many symplectic schemes are applicable only when the Hamiltonian is separable

in a Runge-Kutta form. Another characteristic of our cherne is that it iymplectic with

+ V(q)), or are not of very high order, and often arespect to linear Hamiltonian systems. This is discussed in section 1.6.

h(that is of the form H(p~ q) = T(p)discussion of the numerical stability properLies of the method lacks (a comparison of thecoefficientsthat characterize our methods are given in the Appendix.

accuracies of some symplectic integrators can be found in ref. [30]). While some meritsof symplectic integration schemes are clearly demonstrated (refs.

[16], [17], [20]) in a

4.2

Symplectic Methods

constant tep-size environment, it is still not clear if it will be possible to keep the benefits
Consider the Hamiltonian systemof symplccticity in a variable step-size environment (ref. [22]). With this background in
mind, we think that it is useful to look for new numerical methods applicable to Hamiltoniansystems, as well as to ordinary differential equations.

dqi

oH

dt - OPi'

i = 1,2, ... , N.

(4.1 )

We know that when we transform the coordinates from q = (q], q2, ... ,qn), P = (PJ, P2, "·,Pn)]n tllis chapt r we will describe a new numerical scheme (ref. [24]) for the solution of

ini tial value problems for ordinary differential equations. This work represents an attempt

ical if and only if the condition Mt J M = J is atisfi d. Here M is the' Jacobian of theto modify, and possibly to improve, an algorithm developed by T. Watanabe et al. (refs.

transformation and J is the matrix that has been defined in eq. (2.8). Also, we discuss d in[25], [26]), call'd HID MAS, with an eye towards applications in the field of Hamiltonian

section 2.2 that, if q = q(qo, Po, t), P = p(qo,Po, t) is a solution of Hamilton's equations, thensyst m. The result of our efforts is a new numerical scheme of very high order and inpossess of good stability properties. The principal difference with respect to HIDMAS is

the coordinate transformation (qo,Po) ~ (q,p) is a canonical transformation, and therefore
the symplectic condition J\1 t J M = J is satisfied. In this case the Jacobian is gi ven bythat while in IllDMAS all the grid points are chosen in a way similar to a collocation RungeKuLta, in our schem we have a number of fixed grid points and other grid points chosen
in analogy to HIDMAS. This results in a higher order of accuracy with less computational

In numerical integration methods, the coordinates after an integration step are given by a

eiTort. \~1 will de crib practically our algorithm until the order 14, and we will show thatmapping Wh in phase space that effects the transition (Pn+b qn+l) = Wh(Pn, qn)' Thus theit i A-stable until th ord r 12. In section 4.5 we show how it is possible to put the scheme

method is canonlcal if and only if the Jacobian of the above transformation is such that

-

53

54

CHAPTER 4. SYA1PLECTIC INTEGRATIONWe see that at each time step the value of the energy i multiplied by the factor (1 + h 2 ),the symplectic condition holds, that is it must bethat is the system has been artificially excited by th

num rical n1cthod. This kind of(4.2)
error, increasing steadily with time, is called a secular lTor. The use of a higher order,
where we have introduced the notation z = (q, p). When we are dealing with the numericalnon symplectic method, like the classical fourth order Rung -Kutta, docs not make things

integration of Hamilton's equations, the use of non symplectic numerical methods may

much better. It is in fact found in that case that (ref. [31))have a number of undesired consequences. Let us consider a simple example. For the(4.7)
one-dirnensional harmonic oscillator with Hamiltonian II = (1/2)(p2 + q2) we know theanalytic solution and we also know that the energy is conserved. Suppose now that weso that we have an artificial damping of the system. Th artificial damping is of much morewant to integrate numerically the equation of motion for q and p and let us use the Euler

small amplitude than the artificial excitation provok d by themethod (non sympledic), which for a system of differential equations of the typemany integration steps its effects will become evident anyway. For symp lccLi

a
8x Y = F(y),uler method is, but after
chemes ingeneral the energy is not conserved, but there is no s cular incr as in th error. This is(4.3)

because a symplectic scheme is a Hamiltonian system itself, and therefore the muneri al
trajectory must lie on a surface of constant H N, where HN is the numerical Hamiltonian.

makes usc of the mappingyl = Y + hF(y).

(4.4)

Then, w find that an application of the Euler method to the harmonic oscillator problem
gi yes, aft r on time step of length h

4.3

The Difference Scheme

To illustrate the difference scheme of our method, let us suppose that the equation to
be integrated is of the form yl(X) = F(x, y(x)).

ql = q + hp,

pI = P - hq,

Then, first the integration time step

(4.5)

is subdivided into n subintervals of equal length by placing in the interval the points

0, h, 2h, ... , nh . The number of unknowns is therefore 2(n + 1), since we consider unknowns

and th I' for we have

(4.6)

---

---

the quantities Yo, ylo, y}, yll, ... , yIn. We have one initial condition (thvalue of yo) and

55

56

CHAPTER 4. SYMPLECTIC INTEGRATIONn+ 1 equations ylo = F(yo), ... , yIn = F(Yn), and therefore we need n more relations to close
our scheme. To obtain the n relations we need to choose n points slh, ... , snh, each for every(4.13)
sub-interval, and we write the equations yl(slh) = F(y(sth)), ... ,yl(Snh) = F(y(snh)). Vve
now write the quantities y(slh), yl(Slh), ... , yl(snh) asny( sh) =L [Cj(s)y(j h) + hdj(s)yl(jh)] ,(4 .8)

1 = 0 .. , 2n + 1

j=O

(4.9)

where 8 is the Kronecker symbol, since we need 2(2n + 2) equations in order to determine
all the coefficients in eqs. (4.10) and (4.11). This yields y( h) as an approximation of order

The coeffici nts obey the relations

O(h2n+2) to the true solution and yl(sh) as an approximation of order o(h2n+l). To makethe approximation of y( sh) to be of the same order of tbe approximation of yl( sh) , we
We now determine the coefficients by expansion of the right hand side of the preceding

add one more equation for yl( sh), precisely the equation obtained letting I = 2n + 2 in eq.

equations about the point sh, and then by equating the coefficients of the derivatives of y.

(4.13)The Taylor expansion yi lds(4.10)
Since the coefficients dj ( s) and h( s) have already been determined (as functions of s ) frorn
eqs. (4.12) and (4.13), eq. (4.14) is an equation of degree 2n + 1 in th variable s. The(4.11 )

2n + 1 solutions are all real, and of these n + 1 correspond to the points s = 0,1, ... , n , sincein these points eq. (4.14) is automatically satisfied. The remaining n solutions correspondh coefficients hould be therefore det rmined by the relations

to points which lie in the n sub-intervals, as can be seen from Fig. 4.1. These n solutions(4.12)
give us n sets of coefficients, which characterize the method for the given order n.

57

4.3.1

Implementation

58

CHAPTER 4. SYMPLECTIC INTEGRATION

where

[1

1

1

n
8y(sh) = ~Cj(s) (2n + 2)! (j - s)2n+2 + dj ( ) (2n + l)!(j - s)2n+l ,Once the coefficients have been obtained, we have to solve simultaneously the following
system of equati.onsnIi = 2:: [CikYk + hdikF(Xk) Yk)] ,(4.15)k=O

F( ih, Y:) =are problem independent coefficients. Plots for n=1,2,3,4,5,6 of 8y( h) and 8Yf(sh) aret[*9ikYk + fikF( Xb Yk)] ,shown in Fig. 4.2 and Fig. 4.3.
(4.16)

k=O

with i = 1, .. . ,n, and the notation Y(Sih) = li is used. That is, after substitution ofeq. (4 .1 5) into eq. (4. 16) , we have to solve simultaneously a system of n equations from\vhi ch we get the n values Yb ... , Yn' As discussed in section 4.3 the order of the solution isO(h2n+2), and this is true for every value Yl , ... , Yn, in contrast with Runge-Kutta methods

4.4

Accuracy and Stability of the Method

To test the accuracy and the stability of the method we have appli d it to tb solution of
the test equation

like the Gauss-Legendre method or the Lobatto methods, for which the intermediate values

yf(X) = -AY(X),

y(O) = 1,

A E C,(4.17)

ar of lower order wi th respect to the final value Yn' A more detailed comparison betweent. bese methods and our method is presented in section 4.4.

for the cases n = 1,2,3,4,5,6. Though eq. (4.17) is very simple, its us as a mod 1
to predict the stability behavior of numerical methods for general nonlinear systems is

4.3.2

Truncation Error

Th expre sions for the truncation errors of y , yf arevery widespread. The exact solution of (4.17) is y(x) = exp(-Re(Ax)), and its module isgiven by 1 y( x) 1= exp(-Re(Ax)) . The module of the exact solution is not increasing if

R e( AX) 2: 0, therefore it is natural to require tbat a numerical method, when applied to(4.17), gives
1 Yn+l 1:::;1 Yn I,(4.18)

59

CHAPTER 4. SYNIPLECTIC I TEGRATION60

if Re(),h) 2: o. By definition the region of absolute stability of the method is the region

ourselves to the Gauss-Legendre and the Lobatto methods, it is known (ref. [33]) that< 1. If this region comprises the half-plane Re( h),) > 0, then the method isthe n-stages Gauss-Legendre and the (n + l )-stages LobatLo methods both require thefor which1y1said to be A-stable (see ref. [32], pag. 374-378).

simultaneous solution of n equations and have order of accuracy 2n. ThC'refore, for theVole give below the expression for yl( x) in the case n = 5. The leng th of the step in thiscase is 5h. The solution yl( x) can be written as

same computational effort , our method provides a higher ord r than the aforernentioned
methods. Besides, for our method the order is the same also for all th inLcrmeJiate values,A(),)

in contrast to the Gauss-Legendre and LobaLto methods for which th order drops in they(5h) = A(-)')'

intermediate stages.
where
We note that when), is pure imaginary, then I A~~l)A(),)

1= 1. This is reOect d also in

9979200 - 24948000 h ), + 29862000 h2 ), 2 - 22680000 h3 ), 3
Fig. 4.4, where we plotted the contours of y as a function of the complex argument )"'h+12199320 h4),4 -4904550 h 5 ),5 +

1512075 h6 ),6corresponding to the values 1 y 1 = 0.5,1 , 1.5.

In the figure the line corresponding to

1y 1= 1 is the straight line which corresponds to a pure imaginary )'h. It can ben that

The Taylor expansion of y(5h) - exp( -5h),) yields
the method is A-stable for n = 1,2,3 ,4,5, while the A-stabili ty property is lost for n = 6.
This is a consequence of the fact that the zeros of the denomin ator of the solution cross thev"hich di ffers from th expansion of exp( - 5h),) to the order o( h 13). The accuracy is of thesam ord r also for the intermediate values of y. The discretization for n = 5 is of orderO(h12) , but in this case when Vile solve for y the accuracy is raised by one. So, with respect toth tt equation, for a given n our method is accurate up to o(h 2n+3 ), and we say that the

met.hod ha order 2n + 2. Since for a given n we have to solve n simultaneous equations,as explain d in section 4.3.1, we can compare our method with Runge-Kutta schemes
which r quire the simultan 01.1 solution of the same number of equations. Restricting----

axis Re(h),) = 0 for n = 6 (see Fig. 4.5). Investigation about the A-stability property forn > 6 has not been done, but Fjg. 4.5 suggests that the A-stability property may be lostfor n > 6. This contrasts with Gauss-Legendre and Lobatto methods, which are A-stablefor every order (ref. [33]).

61

5

5

2.5
r-i
O~------+-------t1
r-i-2 . 5

10"'I" a2.5

r-i
x-7.5li-_~~_________~

a-0 . 5(V)I-7.5u-~

-1

0.5

I

a

r-i
x - 0.5
0.5

I

-2

-1

I

a

a

a

r-i

r-i
x - 1.S

r-i

-2

~-1.5

r-!

x

X

f

(\
a

N

V
1

2

-6_ 3

-2

a

0.5

-2

-1

a

1

2

3

a

-1

n=3

n=4

0

x -3
0
n=5

r-!

-2 . 5- 1. 5 -1 - 0.5

3

r--

-1

'"'I a

a

6

-2

a

n=4

1. 5

-3

a

-0.5

n=2

r-i

-1

3

r-i

-1

a

n=3

I

0.5

to

1.5

1

0

-1

n=l

-1.5

- 1

U)

x

a

x

a

X

- 2.5
- 0 .5r-!

-1.5 - 1 - 0 . 5

r-!

a

r-i

a

r-!

.-l

1.5
1

If)

a

. - 1. 5________~____~~__~
-0. 5
a
0.5
1
n=2

1

r-i

I

r-i

n=l

I

<:1'

a

-5

0.5

a

a

a

r-i- 2 . 5

-5

o:r

ClIAPTER 4. SYMPLECTIC INTEGRATION

7.5

7.5
(V)

62

0

r--

r--

I

I

0
r-!

0
r-!

-3

'-0 . 6

r-!

r-!

X

x

n=6

-1
-2

-1

0

2

-3

-2

-1

n=5

Figure 4.1: P lots of the r.h .s. of eq. 4. 14 for n = 1, . . . ,6.
Figure 4.2: Plots of the truncation error for y.

0

n= 6

1

2

3

64

CHAPTER 4. SYMPLECTIC INTEGRATION

63

o

1
0

~

M

10

10

5

5

0

0

-5

-5

I

I

0
r-i

-1e}

-3

r-i

-6

x

-4

-2

-10

4

2

0

6

-3 -2 -1

n=l , y [l )

0

1

4

2

n-2, y(2)

-5

-2

o

- 0. 5

-1

- 0.5

0

0.5

1

10

10

5

5

0

0

-5

-5

n=2

n=l
1
'<:l'

LI)

0

I

0

I

0

0

r-i

r-i

r-i

r-i

.-1

x

x
-1

-2
- 10

- 1. 5

-1

- 0.5

0

1.5

-1

-2

1

0

n=3

-3 -2 - 1

2

(\

0.5
U)

I

o

o
r-i

rl

x

0

/\

o
n=S

3

4

-3 -2 -1

1

2

0

1

2

3

4

n = 4 , y [ 4J

10

10

5

5

o

o

-5

-5

\

-0.5

\

-1.5 v

-6
-2

2

n=4

\.D

-1

1

n=3, y (3)3 ~----------,-------~--,

I

- 10
0

-3

-2

-1

0

n= 6

1

2

3
-10L-____________- 3 -2 - 1 0
1
2~

n=5 , y [l)

~

__
3
4
~

-10L-____________-3 -2 - 1 0
1
2
~

~

__
3
4
~

n -6 , y ( 6)

Figure 4.4: Plots of the contours of 1y I. On the left of each figure the contour corresponding
Figure 4.3: Plots of the truncation error for yl.

to 1y 1= 1.5, in the center 1y 1= 1, on the right 1y 1= 0.5. Horizontal axis: Re()"h). Verticalaxis : I m()"h).

CHAPTER 4. SYMPLECTIC INTEGRATION

66
65

4.5

5

Runge-Kutta Formulation

It is possible to put our method in the form of an implicit Runge-Kutta scheme, and thisis particularly useful because in doing so we will be abl to u e many of the t.heoretical4

results concerning the class of Runge-Kutta type integration methods. Let us consider aCD~system of D differential equations

CD
<D

dy
dt = F(y, L),( 4.19)3CDcD

where y = (Yl, ... ,YD) are the unknown functions and F = (F" ... ,FD). ~or an -stage
Runge-Kutta method the time-stepping from tn t.o t n +1 is given by

(!)

CD

2

(4.20)<D(i)

CD

with i = 1, .. . , s. The vectors Zk are called internal st.ages, and t.he (s + 1) x s r al constants(!)

bk , aik are the coefficients that completely characterize the method. Once the coefficients

CD

1

<r>

<D

of a Runge-Kutta method are given, many theoretical results concerning the stahility,<D

accuracy and other properties of the method are available (see for examp1e ref. [33]). Of

(!)

<D

CD

particular interest for us are results concerning the symplectic properties of the method.
4

-3

-2

-1

o

1

It has been shown in refs. [27], [29] and [28], that if the coefficients satisfy th relation(4.21 )

Figure 4.5: Distribution of the zeros of the denominator of A('\)/ A( -.\), solution of thetest equation yf = -.\y, for n=1, ... ,6. Inside the circles the corresponding n is indicated.
The vertical axis is Re('\h) the horizontal axis is Im('\h). The zeros are symmetrical withrespect to the horizontal axis.

then the method is symplectic.

67

68

CHAPTER 4. SYMPLECTIC INTEGRATIONLet us now show how our method can be put in a Runge-Kutta form. V,/e start by writing
eqs. (4.15) and (4.16) in the form

(4.22)

Now, remembering that Yz = F(xo + slh,}/) y~ = F(xo + {Lh, Y~) and defining the quanti ties Va = Yo, Vi = Yi, i = 1, ... , n, we obtain(4.23)
with the convention that we have to sum from 1 to n if two indexes are equal.

Vj = Yo + h(g;/(F(xa + slh, Yi) - fl~F(xo + {Lh, \/~))),

The

summation is frorn 1 to n if the indexes are Latin and from 0 to n if the indexes are Greeks.
ow, if the matrix gik is invertible, multiplying eq. (4.23) by g;/ we obtain

We obtain therefore the following Runge-Kutta schemethat isVa = Yo·
Sub tituting in eq. (4.22), we obtain
Comparing with (4.20) we see that our scheme can also be consid red as a 2n + 1 stages
Runge-Kutta process, with a step-length nh, and coefficients given by~m;v, jr in eqs. (4.15) and (4.J 6) we consider the function Y = 1, we can derive the relations-1

n bi+n = gni ,-1nai,j+n = gij ,gjilgiO = - 1, CiQ - Cikg;;/glO = 1, j = 1, ... ,n, so that we obtainwith i,j = 1, ... ,n, {L = O, ... ,n and aOi = 0 for i = O, ... ,2n. For n=1 and n=2 thecoefficients are given below. In the last row of the matrices are the coefficients bi. It is69interesting to note that for n = 1, a permutation of the rows and columns shows that our

70

4.6

CHAPTER 4. SYMPLECTIC INTEGRATION

Linear Symplectic Methods

method is equivalent to a Runge-KuLLa of the Lobatto III-A type with 3 stages.
In this section we show that our method is indeed a linear sympl ctic method. By linear
symplectic methods we intend integration methods which are symplectic only for the class

o

o

0

1/6

1/6

2/3

of Hamiltonians of the form H

2n X 2n square matrix. In the case of two degrees of freedom thi expression r duc s to

5/24 -1/24 1/3
1/6

1/6

= zTRz, where z = (q1)'" qn,Pll ... ,Pn) and R is a

H = ap2 + f3q2 + ,qp. The one-dimensional harmonic oscillator is a member of this class.

2/3

Let us now proceed to study the linear symplectic properties of our method. The systemof equations that we have to solve is

0

0

0

0

0

31
480

2
15

1
480

3 (8+5V31
160

3 (8-5V31
160

.l.

4
15

1
15

3
10

10

2~9-7v32-9±2v3l§.:tfl

135

1080

120

3 ~ 2-v31
40

81-2 v31080

2 ~ 9+7 V31
135-{ 9+ 2V3
1080

.l.

.i..

.1..

15

81±2 v3
1080

15

n

Yi = ~ [CikYk + hdikF(Xk, Yk)] ,(4.24)

k=O

3

1 3 ~2+V31
40

18-0
120

.l..

.l..

(4.25)
and, when applied to the Hamilton equations, we obtain

10

15

15

10

n

iii =
ow, jf VI'e apply the condiLion 4.21, we see that our method cannot be proven to be

L [Cikqk + hdikG( qk, Pk)] ,(4.26)

k=O

s mplectic. In fact, for n = 1 we obtain(4.27)

-1/36

o

1/36

o

1/36

-1/36
n

1/36

-1/36

o

Pi = L [CikPk + hdikF( qk, Pk)] ,k=O(4.28)

71

CIIAPTER 4. SYMPLECTIC INTEGRATION

72

obtaining(4.36)(4.29)
where we have used the notation p = F(q,p), q = G(q,p) and qi,

qi, Pi, Pi correspond to(4.37)Yi, Yi. Differentiating we obtainfor i = 1, ... , n. Now, for a symplectic method we must have that the relation dqn 1\ dPn =(4.30)

dqo 1\ dpo is satisfied, that is it must be(4.3 )(4.31)

The coefficients a~l' a~2' a~2' a~l depend on the coeffici nts Cik, dib 9ik, fik' and on the deri va(4.32).
of , ~of , ~of , ~,~,~,~
of oe Be Be , ~Be an d t h e reatlOn] . (38)'. C d f or our mel h0 d.
tlves
~
4.
IS not satlsllCvqk
VPk vqk
VPk
vqk
VPk
vqk
VPk

However, for a linear Hamiltonian all these derivatives are constants, and using a computer
algebra program we have verified that (4.38) is satisfied for n = 1,2,3 and 4. We can(4.33)

TOW,

therefore conclude that the integration scheme is linearly symplectic for these values of n,
and we conjecture that it is a linear symplectic schem also for higher values of n.

we havc(4.34)(4.35)

4.7

Conclusion

and therefore, substi tuting these relations into (4.30) rv (4.33), and then substituting (4.30)

We have presented in this chapter a new numerical scheme of very high order, suitable [orinto (4..31) and (4 .32) into (4..33) we obtain a linear system of 2n equations in the 2n + 2

the solution of ini tial values problems for ordinary differ ntial equations. 1 he method is Astable until the order twelve, and this, together with considerations about the distribution of

_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _--..1_ -........ __73

the truncation error, suggests that maybe the orders between eight and twelve are the best
choice for many problems. The method has also proven to be symplectic with respect tolinear Hamiltonians, a characteristic that could be useful even when nonl1near Hamiltoniansystems are integrated. Numerical experiments on the symplectic properties of the method
and its application to a magnetic field flow problem will be presented in the next chapter.
In the Appendix the coefficients until the order fourteen are explicitly and exactly given .

74

CHAPTER 4. SYMPLECTIC INTEGRATION

76

CHAPTER 5. NUNIERICAL RESULTS

body Kepler motion problem on the vibrating beam Hamiltonian (ref. [16]) and on the
two-resonances Hamiltonian H = p2/2 + 27rc( cos(27rq) + cos(27r( q - t))). In the case of
the Kepler problem no secular error on the energy and the angular momentum has been
observed. In the case of the vibrating beams Hamiltonian no secular error on the nergy

Chapter 5

has appeared.

Numerical Results

5.2

Toroidal HamiltonianIn this section we present the numerical results concerning th

Lhird order Bami Itonian

derived in section 3.5. At the first and second order th Hamiltonian do

5.1

Introduction

not reproduc

in a satisfactory manner the chaotic behavior of the toroidal helical magncLi.c field, so we
describe here the third order results. The Hamiltonian contains two smallness parameters, c

Tn this chapter we apply our method to the solution of the Hamilton equations obtainedand E', connected respectively to the inverse aspect ratio and to the intensity of the vertical[or the magnebc field lines system. The numerical results described in this chapter are
field, which we can vary in each calculation. For a given choice of E and E', we present the
obt.ained using th coefficients corresponding to n = 4, which give a method of order 10.
Poincare section taken at <p = 0, that is we plot the intersection points of the magnetic IIeld
These coefficients are given xplicitly in the Appendix. The method is then implementedlines with the plane X Z. The equations for the magnetic field lines are obtained applying
in a program ,·v ritten in RATFOR (ref.

[34]). The step length h is held fixed, h =Hamilton's equations

2.4'" xI 0- 2and for every time-step the program proceeds using the Newton-Raphson

method Lo solve Lh system of equations given by (4.15) and (4.16). For a Hamiltonian
sy ten1 of n degrees of fr.

8f) = 8<p H,.

87j; = --Ii,

8cp

dom this implies that the Newton-Raphson routine must solve

to the third order Hamiltonian derived in section 3.5. The points are plotted at intervals

a sy tern of 4 x 2n equations at every tim -step. VYe tested the algorithm on the two-in <p of length (27r )/5 because of the periodicity of the Hamiltonian. This results in 5 times

75

77

78

'Ij; . The quite scattered points after '!jJ

CHAPTER 5. NUMERICAL RESULTS

f'V

0.1 are due to the fact that the magnetic field

less computer tlme for the same amounts of plot points.
line escapes after a few revolutions. The power law decrease of t agrees quite nicely with
For small values of E and E' the Poincare section does not show any sign of chaotic

the t heoretical prediction of ref. [35]. W can conclude that the third order Hamiltonian

behavior. Typically there is formatlon of magnetic islands as the radial dist ance 'ljJ increases,

derived in section 3.5 reproduces reasonably well the behavior of Lh field lines system andbut when 'Ij; reaches a certain critical value the orbits tend to escape abruptly instead ofthat it can be used as a basis for more detailed numerical investigations.

giving thc characteri stic random distri bution of t he chaotic domain. Thi s can be seen in
Fig. 5.1. The straight line on the right represents the chosen initial conditions. At a
ccrtain criti cal ini tial condit ion the magnetic line does not wind around the torus forrrilnga magnetic surface, but escap es abruptly and is not even seen in the plot , except for a few
scattered points . Tn Fig. 5.2 we can observe the same behavior. We plot here starting fromgroups of scvcn close initial conditions, that is for each initial value of 'Ij; we choose 7 close
valu cs of 0, so that each line in Fig. 5.2 corresponds in reality to 7 different lines. As t or(' or both are in creased, this b havior gradually changes and the formation of a stochastic
layer begins to take place (Figs. 5.3 , 5.4 and 5.5). For E = 0.023 and E' = 0.1 the chaotic
region appears clearly. In Fig. 5.6 a detail of the Poincare section is shown. The fivet raighL lines on t he top right are the initial conditions. Also in this case we have orbits
confmcd on magnetic surfaces and orbits that escape out of the plot when 'Ij; is big enough,but in b twe n ther are orbi ts that intersect the X Z plane in many random points before
eventu ally e cap e. In Fig. 5. 7 we plot only three orbits. While one orbit originates anisland stru cturc, the oth r two wander around quite randomly. In Figs. 5.8 and 5.9 two
examples of th rot at ional transform t are presented. The horizontal axis corresponds to80

CHAPTER 5. NUMERICAL RESULTS

79
0.3
O. )

O. ,

0.'

O. ,

o· ,

o0
O· 0

-0.1
-0. I

-0.2

- O. J- 0.3"L---~-~-------:t.--,--Q':o"------;;'-;--------;;'-:::-----------':-'-:----~- o.-a.lI.

-0.3

-0.2

-0.\

Figure 5.3: Poincare' section for t
Figure 5.]: Poincare' section for t

0.023, t'

0.0

0·1

0.0245, ('

0.2

0 .3

0.4.

0.09. The straight line on the right0.02. The straight line on the right______~________. -______- .________~______- ,________~__- .________

~

-

-

-

-

o .•

~

represents the initial conditions.
represents the initial conditions.

O. J

,.

0.6

0.2
O. "

,

0.'

.'-:.

~,j~~'~~;:~~~~; '.
,>~.

. ~ .. w-."o. ,:::-"'

~\~.-

\~KV"

0.0

~.,:.~:;'!

0.0

-0.1

-0.'

J

~rf~:r':. . ·~,;~l~ ~.A7#~P··~y~~0~";~~"

-0.2
-0.11

-0.3

-0. ,

-0

. , for groups 0 f 7 c1ose initial conditions .
Figure 5.2: Poincal'e' sectl'on for' " = 0.02, 1
c' --. .0.02
I...

:~~.'~-----=<-0~.3)-------:---:0.22-------:-:O0~.,------00':<.0)----,0;':-,"----~0-:-.-:;-2--------:'0.-:-)------~O.•

Figure 5.4: Poincare' section for t = 0.0245, t' = 0,09, for groups of 7 close initial conditlons.

81
O

"

~

-.r-________________' -________- r________~________- r________- .__________r -______- - ,

82
O. J

CHAPTER 5. NUMER.ICAL RESULTS
o .•

0.1

aj

O. ,

- "

a .1

. .~ "> ~~·· · t·,·· • .
I

0. 0

o· \~

/

I,

,

I

:.c:'"\-0.2~ ';: 0 •'- 0. I-. ...

-0 . J-a.2-O;6L.•----------~O.~J--------~-O~.~2---------70~.,r-------~0~.70--------~OL.~,--------~0.~2--------~Ot.~3--------~0. •-a.3

Figure 5.5: Poincare' section for E = 0.023, E' = 0.09, for groups of 7 close initial conditions.

-::-. ~

-IJ.l..........

'

..'

0, a

-0. I

O. I

a .1

O. J

a .•

Figure 5.7: Plot of three initial conditions, E = 0.023, (' = 0.1.

.....
....
..

-0,2

-a. J

.-

,. '

6. a1 r---------.--------.---------.---------.--------'I---------.________-.________-._________

0.25

....
...:........ ..
_.... . . . ":••~~I.

.....

• •. j -: ,:..•...

,:"

••

..~

"

.....

..

-

,.'

-

......
'....

..:...... .

~. a '" •••••••••••
.

....-

..........

-. .• · · 0 '11....... "... ..... ".

-

. ..

' .0-

-

-

~ .~ ~ ~'0t'.,fl1!:1:~ !;~I' \1:ir-,:, ~(, ,·J~~b'{'~: ,.:;,~..' ,.; '),\\:.,,--...

..

I·"~!~};f~~;~t~~I~~~~~~f~i~i:~,;,\i:~
1,1,:;:>:>,:\;:

0_''l5;f{~,(~,·!:\b. ,t!o-.

:~"(,-,,\-!.I,~."!.:?:. ,-~-. - ..

'-~t~~~!~~~~~;~ft~~~ -';:~.;~::/;~~R;~;~·· :;' ,: . ;': :, .

r

0.\

-

.-" -:. .;,.?@;

0.12

O.I~

0.16

0.19

Figure 5.6: Poincare' secLion [or (; = 0.023, (;'

0.2

0.22

0 · 26

0.28

-

O.a-

0.1 (enlargement). The four straight linesI

correspond to the initial conditions.

- I . gL.o::--------~a~.0:::1------~a.~o:. ------~o-;-.~o',.-------;;0':1.o"s------r:o.:1,---------;0-:..,"2------'at.-,;";-----oa:-:-.1<"'------,:0. 18

Figure 5.8: Rotational transform [or E = 0.03, (' = 0.04.

83

84

- : - - - - - . - - - - - , r - - - - . - - - - - , - - - - - . - - - - r - - - - . - - -i

6 .0r--_ _:~-.... . ...............-......................-.....-........... .. .........
"

.'

3.0

",

2·0

, .0

0.0

.--;-;'2:----0r::-.;;-,.--OO~.,"- - - r / O. I a·1. g,L.
0 - - - - 1 0.-02-----:1. ::-:-.---::01-::.
0-:-6- - - ! - O.-;:;;os---;o;L;.' - - - - - ; ; 0
0 0Figure 5.9: Rotational transform for (; = 0.03, (;' = 0.07.

CHA PTER 5. NUAtfERICAL RESULTS86

CIIAPTER 6. CONCLUDING REMARKS

done, and we plan to study the convergence propertie of th

xpansion scheme in the near

future. Also, a renormalization approach which makes use of a more r alistic form for theunperturbed Hamiltonian and gives a better convergence is being currently investigated.
We have then presented the construction of a llnear symplectic integration scheme. The

Chapter 6

scheme is of very high order, and has also good stability properties. Vie have hown that.
it is A-stable until the order twelve. Contrarily to many existing symplectic method ,our method is applicable to every Hamiltonian system 1 and not only to special classes o[Concluding Remarks

Hamiltonian systems. The choice to use a linear symplectic sch m' for th integration of anonlinear system is ultimately judged by numerical experiments. Numerical t.ests on known

The two main results of this thesis are the derivation of a Hamiltonian for the helicalmodels have been carried out and have confirmed its good performances. It application to
the magnetic field lines problem has shown that the Hamiltonian we have derived is ahle

magnetic field in the vacuum and the construction of a new linear symplectic scheme
to reproduce the regular and chaotic behavior of the field.
suita.ble for its numerical integration.
We have deriv d the Hamiltonian applying a perturbation theory based on the Lie transform, a syst matic procedure to calculate the expansion to high orders. After obtaining
a. lIamiltonian for the cylindrical limit approximation of the field, we have derived the
Hamiltonian for the general toroidal case up to the third order. We have also derived an
intcgrabl model for the same toroidal case. The existence of these two different Hamiltonians [or the same physical system is related with the choice of the variables when weapply th Lie perturbation t chnique, and probably also with the asymptotic character ofthe xpansion itself. Research on the relation between these two models is being currently
85

IThe scheme is also applicable to differential equations in general.

APPENDIX A. COEFFICIENTS

88

A.2

The case n=I

For n = 1 we obtain an integration scheme of th fourth order. Eqs. (4.12) and (4.13) now
read
(A.I)

Appendix A(A.2)
where I = 0,1,2,3 . Solving these equations, for th new variable p = s + l/2, we obtain

Coefficientsthe coefficientsco(p) = 1/2 - (3/2)p + 2p3,

Cl(P) = 1/2 + (3/2)p - 2p3 = eo( - p) ,do(p) = 1/8 - p/4 - p2/2 + p3, d] (p) = - 1/8 - p/4 + p2/2 + p3 = - do( -p),A.I

Introduction

9o(P) = -3/2 + 6p2 ,91(P) = 3/2 - 6p2 = - 90( - p),fo(p) = -1/4 - p + 3p2,

f1(p) = -1/4 + p + 3p2 = fo( -p).

(A.3)

Substituting the coefficients into eq. (4 .14) we obtain the equation
In this appendix we giv the explicit form of the coefficients of the integration scheme
p3 _ p/4 = 0,

for n = 1 ... ,6. While in th case n = 1 it is still possible to calculate the coefficients
1 y hands,

a n increas the algebra is still elementary but quite lengthy and it becomes

useful , if not neces ary, to usc a computer algebra program to perform automatically thewhich has solutions p = -1/2,0,1/2. The solution p =87gives us therefore the set of

d1 =_18'

90 = -~, 91 =~ ,fo =-i, fl = - i,which will be used to solve the differential equation.

cas's we giv only the r suIts.

°

coefficientsanalytical calculations. We have used the program MATHEMATICA. In what follows , thecase n = 1 is treat d in detail , wi th the purpose of illustration, while for the remaining

(A.4)

(A.5)

· --

-.

.

~

-~ . -

-

-

~

--

-. -

89

A.3

The cases n

90

APPENDIX A. COEFFICIENTS

2, 3, 4, 5, 6
eo(p) = 5~~p,

The results are obtained following the same procedure described for the case n = 1. Wegive below D(p), which is defined to be tbe right hand side of eq. (4.14), and the solutions

d(o P )-2.=£- 18'

of D(p) = 0 which lie in the sub-intervals, and the coefficients for n = 2, 3, 4, 5. The
variable p is defined as p = s - nj2.

A.3.1

n

(A.7)

2

gl ()
P = ~
3 '

D ( ) = (1 - p) p (1 + p) (-1 + 3 p2)

P

360

'

(A.6)

fo(p) = -~ +~ , h(p) = fo( -p),

JI(p) = -~ .

A.3.2
The coefficients are given below. Attention should be paid to the fact that the expressionsn == 3
D( ) = P (-3 + 2p) (-1 + 2p) (1 + 2p) (3 + 2p) (5 - 4p2)
322560P

b low arc valid only for the particular values PI, P2 , and not for every p in the interval -1,1.

PI P2·

This is not a loss of generality, because the coefficients which characterize themethod ar the ones corresponding to the solutions PI, P2. SO, the first set of coefficients-v'5

For the solution P2 we haveeo(P2) = ~,1

Cl (P2)

dO(P2) = 2 4'

= ~,C2(P2) = ~,d1 (P2) = -~, d2(P2) = ~,C3(P2) = ~,d3(P2) = - d4'(A.9)

is obtained substituting in the xpressions below PI for p, and the second substi tuting P2[or p.

be same kind of remark applies also to the cases n = 3, 4, 5, 6, with the obvious

cliff renee that the numb r of roots incr as s.

(A.8)

PI = -2-' P2 = 0, P3 = - Pl.

This is a consequence of the fact that the full expressions of the coefficients are in generalquite long and we have chosen to use a simplified expression valid only for the solutions'

185
gl(P2) = ~,go (P2) = -216'

fo(P2) = -185
g2(P2) = -~, g3( P2 ) = 216'

:2' Jl(P2) -~, h(P2) -~, h(P2)=-:2==

91

APPENDIX A. COEFFICIENTS

92

The coefficients are

For the solutions PI) P3 we haveCD ( P)-1_£.E- 4
27'

C3(P) = CD( -p),C1(P) = t,

C2(P) = C1(-P),

eo(p) = -416+220p+2060p2_1105 p3
21600

c (p) = 404-560p-95 p2+260p3
1f6,d3(p) = - do ( - p),d1(p) = -~ +~,d2(p) = -d1( -p),do(p) = i4 -

9 (p) = -185+162p

93 (p) = -90 ( -91 (p) = 5 -86 p ,92(P) = -91 (-p),

A.3.3

216'

1350'

_ 44-Sp2C2 (P) -100'

d (p) = -8+4 p+50 p2 - 25 p3o(A.10)
o.

3600'

d (p) = 52 -52 p-55 p2 +55 p3

P),

1

450'

d ( ) = P(44-5P2)
2 P

100

'

Io(p) = -7i;6 P ,h(p) = Io( -p),

= 444-256p-1587 p2+9S2 p
9 (p)I] (p) = -3t 2p ,h(p) = I1( -p) .

91 (p)

3on == 4f

JO

D( ) = (-2 + p)( -1 + p)p(l + p)(2 + p)( -4 + 15p2 - 5p4)
P
1814400
'

4320'

270'

(p) = 36-20p-13Sp2+80p3

I 1 (p)(A.11)= -276+448 p+159 p2 -208 p3

3600'= -48+100 p-60p2+5 p3
450'

f 2 (P) = -44+Sp2
100'
Pl = -

15+vT4s
10,P2 -yhs-vT4s\jlo,( A.12)A.3.4

n == 5P3 = -P2,
The numerical values are

D(p) =P (-5+2p) (-3+2p) (-1 +2p) (1+2 p) (3+2p) (5+2p) x
PI = - 1.64443· . "P2 = -0.543912· ...

( - 259 + 280 p2 - 48 p4) / (122621409600) ,(A.l:3)

93

94

APPENDIX A. COEFFICIENTS

eo(p) = -306125+12887Op+282500p2_119096p3
4374000(A.l4)
1

34992

C (p)
2

The numerical values are

C5(P)= eo( - p),(p) = 12825-11354p-3348p2+3848 p3

C

= 0,

P3

'

'

= 445-392p+68 p2+176p3

2187'

d (p) = -2675+1070p+2540p2 -1016p3
o291600'd5 (p) = - do( - p),

d (p) = 2103-1402p-1212p2+808 p3
1

PI = -2.16345·· "P2= -1.0737···.

11664'

d (p) = -83+166p-52 p2+104 p3
2

1458(A.L6)(p) =9o

The coefficients are, for the solution P3

9 1 (p),1<17
eo (P3) = 262144'() _Cl P3 -

10625
262144'

()
C2 P3

1875
= 4096'9 2 (p)() _C3 P3 -(
) _C4 P3 -

10625
262]44'

()
C5 P3

447
= 262144'f
JO

dO(P3) = 13i~72'

d 1 (P3 ) --1875
131072'd ()5625= 32768'

I 1 (p)d3 (P3 ) --d4 (P3) = -131072'
1875
d5 (P3) = - 13i~72'

1875
4096'

5625
- 32768'

2 P3

4375
= - 131072'() _g2 P3 -

13125
8192 '

483= - 655360'

()
g1 P3

()
93 P3

= 13125
8192'

(
) _
g4 P3 -

Io(P3)= - 65~36'

11 (P3)= - 6~;;6'f 2 (P3 ) -- - 16384
5625'

f 3 ()
5625
P3 = - 16384'

14 (P3)= - 6~;~6'

I5(P3)

4375
13]072'

()
Cq P3

-

483
= 655360'= - 65;36'

g5(P)= - go( - p) ,= -88235+7695Op+25244 p2 -20088p3
69984'= 8015-8430p-I724p2+1464 p3
8748'(p) = 16385-7410p-14468p2+6600p3

583200'!s(p) = Jo(-p),= 563+1494p-302Op2+936p

323328'

h(p) = fz( -p).

(A .l5)()go P32602985-118875Op-2280308p2+1047000p3
8748000
'

95

A.3.5

96

APPENDIX A. COEFFICIENTSn == 6
p3 -134661 p4+50023 p5
go(p) = -96588+36036p+377661 p2-139831
g6(P) = - go( - p),4536000'

D(p) = (-3+p) (-2+p) (-l+p)p (l+p) (2+p) (3+p) x_

547-316p-2184p2+1336p3+309t-188p5
3500
,

_

g3 ( p ) -2 p (-981+506 p2_53 p4)567,.( ( ) =J0 P-2028+756p+7861 p2-2891 p3-2821 p4+1043p!>
1058400'

g1 ( P) -g5(P) = -gl(-P),(36 -147p2 + 70p4 -7 p6)/(43589145600),(A.18)PI = -2.67804· ", P2 = -1.59578· ·· , P3 = - 0.530655·· "P6 = -Pl'f 3 (p)

The coefficienLs are_
Co ( P) -27780-9660p-98987 p2+34349 p3+36743p4.- 12761 p5
5292000'_
Cl ( P ) --415+245p+2366 p2_1393 p 3_399p4.+252 p5
12250'~(p) = C()( -p),

(A.17)d ( ) = 3(8-8p-7 p2+7 p 3 )2 p
196= 4 (-159+77 p2_14 p4)1323'

[10] J.R. Cary, Phys. Rep. 79 No.2 (1981) 129[11] J.R. Cary and R.G. Littlejohn , Ann. Phys. ] 51 (1983) 1[12] P.M. Morse and H. Feshbach, Methods of Theoretical Physics, McGraw-Hill, New

References

York 1953[13] J. R. Cary, Phys. Fluids 27(1984) 119[1] M.N. Rosenbluth, R.Z. Sagdeev, J.B. Taylor and G.M. Zaslavsky, Nucl. Fusion 6,[14] G. Gnudi and T. Hatori J. Phys. Soc. Jpn. 62 No.6, (1993) 2030

(1966) 297[15] A.J. Lichtenberg and M.A. Lieberman, Regular and Stochastic 1otion, Spring r[2]J.N. Filom nko, R.Z. Sagdeev and G.M. Zaslavsky, Nucl. Fusion 7, (1967) 253

Verlag (1983)[3] Y. Tomita, Y. Nomura, II. Momota and R. Itatani, J. Phys. Soc. Jpn. 44, (1978) 637[16] D. Okunbor, Physica D 60, (1992) 314[4] T. Hatori and T. \,yatanabe, "Nonlinear Dynamics and Particle Acceleration" AlP[17] J . Candy and W. Rozmus, J. Compo Phys. 92, (1991) 230

1

0n fcrence Proceedings No. 230, Particles and Fields Series 45 (1990) 79[18} M. Suzuki, Phys. Lett. A 146, (1990) 319[5] H. Goldstein, Classical 1echanics, Addison-Wesley (1950)[19] H. Yoshida, Phys. Lett. A 150, (1990) 262[6] G. Hori, Publ. Astr. Soc. Jpn. 18, (1966) 287[20] P.J. Channell and C. Scovel, Nonlinearity 3, (1990) 231[7] A. D prit, Cel. Mcch. 1, (1969) 12[21} E. Forest and R . D. Ruth, Physica D 43, (1990) 105[8] A. Kamel, Cel. Mech. 1 (1969) 190[22] J.I\1. Sanz-Serna, Physica D 60, (1992) 293[9} J. Henrard, Cel. M ch. 3, (1970) 107[23] J.E. Marsden and Z. Ge, Phys. Lett. A 133 (198 ) 134
97

98

[24] G. Gnudi and T . \i\Tatanabe J . Phys. Soc. Jpn. 62 No. 10, (1993) 3492
[25] T. WaLanabe and M. Takagi, Trans. JP

Soc. Ind. and Appl. Mat . 1, (1991) 135 (inJapanese)[26] K. Abe, A. Ishi da, T. \i\Tatanabe, Y. Kanada, K. Nishikawa, Kakuyugo Kenkyu 57(198 7) 8.5 (in Japanese)[27] J.M. Sanz-Sern a, BIT 28, (1988) 877[28] F .M. Las agni , J . Appl. Math. and Phys . (ZAMP) 39 , (1988) 952[29] Yu B. Suri s, U.S.S.R. Comput. Maths. Math. Phys. 29 No.1 (1989) 138[30] R.I McLachlan and Paul ALela, Nonlinearity 5, (1992) 541[31 ] H. Yoshida, Cel. Mech. and Dyn. Astr. 56 (1993) 27[32] G. Dahlquist and A. Bjorck , Numerical Methods , Prentice-Hall,ew Jersey (1974)[33] IC Dekker and J.G. Verwer , Stability of Runge-Kutta Methods for Stiff Nonlinear
Di:lIerential Equations,orth- Holland (1984)[34] B. "'.,T. I ernighan and P.J. Plauger, Software Tools, Bell Tel. Lab. Inc. and Yourdon
In c. (1976) (in Japanese)[35] Y. Ai zawa, Y. Kiku chi T. Harayama, K. Yamamoto, M. Ota and K. Tanaka Prog.
T heor. Ph. Suppl ement No. 9, (19 9) 36
99

• Hamiltonian
Vacuum

f
o
rt
h
eT
o
r
o
i
d
a
lH
e
l
i
c
a
lMagneticF
i
e
l
dL
i
n
e
si
nt
h
e

(真空中のヘリカル系トーラス磁力線のハミルトニアン)
共著者羽鳥予承

J
.P
h
y
s
.S
o
c
. Jpn・ぅ Vo.
l6
2，
N
o
. 6(
J
u
n
e1
9
9
3
)

(日本物理学会)• A NewNumericalI
n
t
e
g
r
a
t
i
o
nSchemeofVeryHighOrderandA
‑
S
t
a
b
l
e

(A一安定性を有する超高精度の新しい数値積分法)

共著者渡辺二太
J
.P
h
y
s
.S
o
c
. Jpn
，・ Vo.
l6
2，
N
o
.1
0(
O
c
t
o
b
e
r1
9
9
3
)

(日本物理学会)