Symplectic integrators and the problem of wave propagation in layered media Текст научной статьи по специальности «Физика»

Математическое моделирование

UDC 517.958

Symplectic Integrators and the Problem of Wave Propagation

in Layered Media

M. N. Gevorkyan, J. V. Gladysheva

Telecommunication Systems Department Peoples' Friendship University of Russia 6, Miklukho-Maklaya str., Moscow, 117198, Russia

In this paper numerical methods that preserve the symplectic structure of the Hamiltonian systems are considered. Hamiltonian is constructed for the propagation of electromagnetic waves in a stratified medium without any sources. Hamilton's equations are solved using symplectic second-order Runge-Kutta method.

Key words and phrases: symplectic integrators, symplectic structure, Hamiltonian formalism, Maxwell's equations without sources.

1. Introduction

In this paper, we use a geometric approach to the problem of dynamical systems. In the framework of geometric approach, the basic physical concepts follow from the various transformations and their invariants. For example, when the Hamiltonian formalism is formulated in the language of the symplectic geometry, the symplectic form oj plays a key role.

A lot of numerical methods have been developed to solve the ODE. However, not all of these methods take into account the geometric (symplectic) structure. It is not important for the local solution (short time period). In the case of global solution (long time period), preservation of the symplectic structure is very important because it can give more information about the system, than the equations themselves.

Fortunately, the numerical schemes, known as variational integrators have been developed. Such numerical methods preserve the global structures and give high accuracy.

In the first section of this paper the Hamiltonian formalism is formulated in the language of the symplectic geometry. Then with the use of the simple model — linear oscillator — we illustrate that classical Euler method does not conserve the energy H of the system. After that, the special case of variational integrator is described. In the last section we use it to solve Maxwell's equations without sources.

2. Hamiltonian Formalism with Respect of Cotangent Bundle

and Symplectic Structure

2.1. Lagrangian and Hamiltonian Functions

For a neighborhood U of point P of the bounded manifold M there exists a coordinate map with coordinates q1,... ,qn. In the following when we mention coordinates on manifold M, we mean local coordinates of some point P E M. M is called configu-rational manifold. Some mechanical system has N coordinates ql = (q1(t),... ,qN(t)). Evolution in time is given by curve ql(t) on the manifold M. Hamiltonian formalism starts with Lagrangian L(ql,ql), and momentum definition [1]:

= 0L(g',g') = df (i)

Pz = d? , q = di . (1)

Received 1st November, 2011.

Now we can calculate the Hamiltonian using the usual definition of H as the Legendre transformation of L

particularly strongly nonsingular, when the equation (1) may be continuously and one-one resolved in the form of q1 = vl(ql,pi), "iql,ql.

If Lagrangian L is strongly nonsingular, the Euler-Lagrange equations and the Hamiltonian equations are equivalent [2]:

Let v(t) = (ql,...,qN)T. v(t) is a vector field given on manifold M (the cross-section of the tangent bundle to a manifold M), and p(t) = (p\,... ,pn) is a 1-form field given on manifold M (the section of cotangent bundle to a manifold M), thus:

— the Lagrangian L(ql,ql) is the function on the tangent bundle TM,

— the Hamiltonian H(ql,Pi) is the function on the cotangent bundle T*M.

Manifold T*M is called phase space. Phase space is manifold and it's dimension is 2N.

Def 2. A symplectic form on a manifold M is a closed non-degenerate differential

def ~ •

2-form Co = dql A dpj.

A symplectic manifold consists of a pair (T*M, ¿D), a manifold T*M and a symplectic form Cj. Assigning a symplectic form Co to a manifold T*M is referred to as giving M a symplectic structure. In the following we will consider only 2 dimension phase space with symplectic form Co = dp A d^. The curve p = p(t), q = q(t) is considered. Tangent vector to this curve is

Now as dw = ddp A dq + dp A ddg = 0, then Lie derivative La of the symplectic form along the vector field u is 0.

Def 3. Vector field u, that satisfied L^Co = 0 is called Hamiltonian vector field.

To find the convolution Cj(u) the definition Co = dq ® <ip — <ip ® dq, thus should be used:

Co(u) = dq(u)dp — <ip(u)<iq = (d^, u)dp — (dp, u)(°q =

_d 0L 0L _ dpi _ d^'

dtW - W _0 ^ _ - dT' Wi _dT'

2.2. Phase Space and Symplectic Structure

d dq d dp d di dt dq di dp '

Co(Ü' •) _ dH(•) ^ Q(u) _dH ^ dw _ d(dH) _ 0.

2.3. Canonical Transformation

Def 4. Canonical transformation is a transformation that preserves the symplectic structure ti.

The new coordinates P = P(q,p) and Q = Q(q,p) will be canonical if dq A dp = dQ A dP,

-¡^ -¡^ dQ~\ idP~ dP~\ dQdP ~ ~ dQ dP ~ ~

dQ A dP = ^ d, + ^ ^ ^dq + ^ = ^ A + ^dq A dp+

dQdP d d dQdP d d f dQ dP dQ dP\ d d

+ "5"" TTdP A dq + -g- -¿—dp A dq = —5--dg A dp.

dp dq op op \oq Op dp dq I

This canonical transformation should satisfy the following condition:

(9QdP_dQdP\ = ^ det d(Q,P) = 1. (2)

I dq dp dp dq I d(q,p) '

2.4. Poisson Brackets and ti

One of the most important points in the geometrical approach to Hamiltonian dynamics is the role of 2-form ti. This form plays the same role in symplectic manifold as a metric tensor in Riemann's manifold. ti imposes one-one mapping for vectors and

one-forms. Let V be a vector field on M. One-form field V can be defined with the following formula:

V = V(■) =f ti(V, ■) = ti(V) — non-coordinate form, (V)i = (ti)ij (V) = tiijV° — coordinate form.

In the same way on the base of the one-form field a can be defined vector field a (one-one definition): a = ti(a, ■) = ti(a).

-t def

Let f and g be functions on the manifold. Vector fields can be introduced: Xf =

^ ^ def ^ ->

df and Xg = dg, where df is the vector gradient. It can be defined from the expression df = Q(df ) = Q(Xg ) = Xg.

Def 5. Following scalar is called Poisson bracket: {f,g} =f ti(Xf ,Xg).

u(Xf, •) = ù(df, •) = df (•) ^ ü(Xf ,Xg ) = u(df, dg) — df (dg) = (df, dg) — (df,Xa ),

—d — dq dp"

Xo = d9 = dQ + dP.

Now, one can find vector Xg :

dg d dg d d ( ^ ^ ) — d f ( Y ) — Y f = dg df dg df

dq dp - d^dq ^ W(A f ,A a ) = d (A a ) = Aa 1 = dqd^p - d^dq '

So, we define Poisson brackets of f and g. The definition clarifies geometric sense of the Poisson brackets. It is easy to see, that they don't depend on the coordinate system but only on two-form w. Canonical transformation does not affect on the symplectic form w, thus Poisson brackets remain the same. A similar statement is true for all constants of motion (zero Poisson bracket also is invariant under canonical transform).

3. Numerical Methods and Symplectic Structure

3.1. Family of Runge-Kutta Methods

Anyone who has dealt with ODEs is aware that for the most of them it is not always possible to find an analytical solution. That's why we have to use various numerical methods to get a solution. One of the most well-known numerical schemes may be the family of Runge-Kutta methods, which can be defined by the following formulas [3]. There is given the array of the coefficients c2,... ,cs; bi,... ,bs; [aij],i,j = 1, For each step s + 1 one should calculate:

s.

ki(h) = hf (x + Ci, y + ^ aijki(h)), y(xi+i) = y(x3) + biks(h).

3 = 1

i=1

It is important to note, that the values of the coefficients c1,... ,cs; b1,... ,bs; [aij] are chosen based on the traditions and practical use of methods [3]. These data (the coefficients) are usually arranged in a mnemonic device, known as a Butcher tableau.

Ci an a-12 0>13 0,14 . . au

C2 0>21 0>22 a23 0,24 . . a,2s

c3 131 132 a33 0,34 . . 03s

C4 a 4i a>42 a 43 O44 . . a 4S

c5 a>51 0*52 a>53 054 . . 05S

Cs a si aS2 aS3 aS4 . . 0Ss

b1 b2 b3 b4 . bs

c A

bT

The above Butcher tableau describes the family of implicit Runge-Kutta methods. But more often the explicit methods are used, where 0 < i < j < s:

ki{h) = hf {x,y),

k2(h) = hf (x + C2,y + a,2iki(h)),

k3(h) = hf (x + c3,y + a31k1(h) + a32k2(h)),

ki(h) = hf (x + C4,y + (i4iki(h) + a42k2(h) + a,43k3(h)),

ks(h) = hf (x + cs,y + asiki(h) + (h) + as3k3(h) + ... as,s-iks-i(h)),

s

y(xj+i) = y(x3) bMh).

i=1

The Butcher tableau for the explicit method is simplified and all elements of matrix A, for which 0 < i < j ^ s are equal 0.

Tables for the methods of Euler, Runge-Kutta 2nd, 3rd and 4th order for specific values of the coefficients have the following form:

0 0

1

0

1/2 1/2

0 1

0

1/2 1/2

1 -1 2

1/6 4/6 1/6

0

1/2 1/2

1/2 0 1/2 .

1 0 01

1/6 1/3 1/3 1/6

3.2. Runge-Kutta Methods and Symplectic Form

Before turning to the symplectic methods, it is worth to study a simple example — harmonic oscillator with unit mass.

H(",q) = T + 2 H m = -M.

To solve the ODE system the Euler's method should be applied, where t0 < t < T, h = tk — tk-\.

( Qk+i = Qk + hpk, \pk+i = Pk - hqk.

So it can be seen, that the total energy H of the system does not change under the time

(q(t) = Ci cost + C2 sin t, in o

\fl Ci t C2 • : ^H(p, q) = 2 (C2 + C22) = const.

yp(t) = C2 cost — C\ sini 2

But from the discrete system we get + 0,2+i) = (pi + 9fc)(h2 + 1), thus with each

new iteration (pi + q\) increases (h2 + 1) times. So, the total energy of the system is not conserved, what can be clearily seen on the Figure 1.

Figure 1. Phase portrait for the exact solution and the solution by Euler's method

<7(0) = 2, v(0) = 1

The above example shows that some numerical methods don't respect global characteristics and structures of the problem. In our case, this structure is a symplectic form. We mentioned its importance above. There is a problem of finding methods, which respect a symplectic form. Such methods exist and are called symplectic integrators [4-7].

3.3. Symplectic Integrators

Symplectic integrators usually belong to the family of implicit numerical methods. Software implementation of such methods is difficult. However, for separable Hamiltonian function H(p, q): T*M — R

H(p, q) = T(p) + U(q), q = (q1 (t),... ,qN(t)), p = (Pl(t),... ,pN(t)), (3)

it is possible to construct explicit numerical scheme. The separable Hamiltonian describes a conservative dynamic system of Newton mechanics. In the framework of relativistic mechanics it couldn't be represented. We divide the segment [to, to + T] into m equal parts with step h = ti+1 — ti, i = 1,... ,m. At each step h some auxiliary quantities are calculated. Let (p(io), q(io)) = (p0, q0). Before a step i + 1, n auxiliary quantities should be calculated.

(p^ qi) — (pi+1 , qi+1) — ... ^ (pi+-L, qi+_L) ^ ...I^(p(ti+l), q(ti+i)).

1 n 1 n 1 n 1 n \_ _/

V

Intermediate calculations

The scheme is very similar to Runge-Kutta method (for RK method one calculates ks(h) at every step). That's why this symplectic integrator is called symplectic

Runge-Kutta method. But the main point is that every transformation (p i, qi) -+1 (piq«+1) is canonical. In the article [5] for the case (3) the second, third and fourth order symplectic integrators were obtained. We write them for the 2D phase space (I — order of method) as following:

qi = qo + amVPmT(qm), i = 1,...,l

m=1

Pi = Po — bqm-iu Gv-i), i = 1,..., I = 2, 3,4. We write out formulas for the second order integrator [5]:

Pk+2 = Pk - Mi %(,k),

Qk+2 7 dT ! \ = Qk + hai dp [Pk+2),

Pk+i = Pk - h \ dU , , , dU, / bi (Qk) + b2 {%+ 2) ,

Qk+1 = Qk + h dT dT °1 dp (pk+ i) + a2 dp (pk+1)

Coefficients (a1, a2,b1,b2) are not uniquely determined. They can be calculated from the indeterminate system of equations. In the paper [5] it is noted that the following two cases are of the most interest:

(a1,a2,b1,b2) = Q, 2,0, l) — leapfrog,

(a1,a2,b1,b2) = ^l, 0,1, ^ —pseudo-leapfrog.

For the leapfrog case the system we have: ' Pk+i = Pk,

To prove the canonicity of the transformation from (pk, Qk) to (pk+i, Qk+i) the formule (2) should be used (for the leapfrog case):

Exactly in the same way we prove the canonicity of the transformation from (Pk+i/2, Qk+i/2) to (pk+i, qk+i) and from (pk, qk) to (pk+i, qk+i) for the 4th order method. If we require the time reversibility of the numerical solution, we can determine ai, a2, a3, a4, bi, b2, b3, b4 uniquely [5,7]:

4. Construction of the Hamiltonian Function for a Layered

Medium without Sources

Let us use Maxwell's equations for the isotropic medium without sources insofar

as:

— in case of currents existence we get Hamiltonian formalism with links (see [8]);

— in case of wave propagation in waveguide there are no currents in waveguide.

The medium which properties are constant on each plane perpendicular to the fixed direction (we get Oz for this direction) is called layered medium. We will consider a plane linearly polarized monochromatic electromagnetic wave propagating in the layered medium.

— The wave is called plane, if the solution of the wave equation has the form E(r-s, t) (the same for H), r — radius vector and s — wave propagation direction. The quantity E(r ■ s, t) for each moment of time is constant on the plane r ■ s = const.

— The wave is called monochromatic if fields vectors are harmonic functions of time.

— When the wave is linearly polarized and it's electric field intensity vector is perpendicular to the incidence plane we will call it Transverse Electric (TE-mode).

— When the wave is linearly polarized and it's magnetic field intensity vector is perpendicular to the incidence plane we will call it Transverse Magnetic (TM-mode).

If the medium is linear, it is possible to decompose any plane polarized wave into two waves, one of which is TE-mode wave and the other one is TM-mode wave. So, we study the plane monochromatic electromagnetic wave. Let us introduce the Cartesian

d qk+1/2 d qk

OP fc+l/2

d qk

dpk+1/2 dp k

d qk+1/2 dp k

coordinate system. The plane of the wave propagation is (xOz) plane and

k = k0n = n—, k0 = — , - = n = \/e(z)a(z), k = ks,

where k is a phase vector, k — wavenumber, e(z) and ¡i(z) - electric and magnetic constants. Assume that e(z) and ¡i(z) change along the z axis. The vector k depends only on z and its component kx is constant. A complex form of the plane monochromatic electromagnetic wave equation is [9,10]

( E = №E0 exp(-i(—t — k ■ r)), E0 = (E0, E0,E0) — complex constant vector,

\ H = №H0 exp(—i(—t — k ■ r)), H0 = ( H0,H0, H0) — complex constant vector ,

( kx, ky, kz) — k £ (xOz) —^ ky — 0 —^ k ■ r — kxX + kyy + k%z — kxX + k%z,

The symbol № is usually omitted during the calculations, because it is much easier to work with complex forms of E and H. After all calculations are done, we can write them out in their real form. The quantities E and H depend only on x and z, thus the following conditions are fulfilled:

oe H Oh h dE dH dH . aE

— = 0, — = 0, — = I kxE, — = I kxH, — = — I—H, — = — I—E. (4)

ay ay ox ox at at

Considering the formula (4) let us write out six equations:

yd H "c dt

' dEE± dEy

dy dz

dEx dEz

d dx

dEy dEx

„ dx dy

e0E c dt

i y k0Hz

( dHz dHy

dy dz

dHx dHz

d dx

dHv dHx

d x dy

= -iek0Ex,

= — i s k0Ey,

= —i e k0Ez.

Considering the equation (4) and the following conditions for the TE-mode: Ex Ez = 0, Hy = 0, and for the TM-mode: Hx = Hz = 0, Ey = 0 we get six equations:

TE-mode: <

9Ey a

dHx a

—i [i koHx,

( dHn

= ikxHz — iek0Ey, TM-mode: <

1 kx

Hz = — -¡—Ey, [ ko

a

dEx a

= i e koEx,

= i kxEz + i[k0Hy

rr — 1 kx n

Ez =---¡—Hy.

s ko y

In each system there is one algebraic equation and it can be used to reduce the number of ODE. Thus we get two systems of two equations.

TE-mode:

0HX a

dE„

= k0

(e — I li)Ey,

TM-mode:

—i [ koHx,

0EX a

dH,

= k0

y = i £ koEx.

dz \ dz

There are two invariant combinations of the electromagnetic field components.

Ii = eE2 — [H2 and I2 = e[E ■ H.

x

y

The Lagrange function L = ^I which is written with the use of the Cartesian components of electromagnetic field has the following form:

L = 1 {e(E2x + E2 + E2z) — „(H2X + Hi + H2Z)) ,

for our case we get L as following

L = 1

2

1 kx\ TTI2 tt2

£-I у-ßH

1

+ 2

- (" -1 % )Hi + 'E

= Li + L2

L\ TE—mode

L2 TM-mode

For the transition to the Hamiltonian function the canonical variables should be chosen. It can be done using the following replacement q1 = aHx and q2 = aHy where a = const G C.

d q\ dHx . —— = a^— = — гакп dz dz

dLi

d

d qi dz

-(g- ßt) Е г ak° (е- Ш)

ЕУ

i akn

dq2 дНу

-= a--

d д

i as k0Ex, p2

д Li

д

dq2 dz

£Ex

i as k0

Now it is possible to write out Hamiltonians Hi and H2 :

Hi

ж-Li = H£-тД)Еу + >0-

H = - L = Uß - 1%\ну + \eEl > 0.

1 k.

1

dz ~2 2 Г ek2) "у ' 2~

Let a =1. For subsequent calculations we should use the real quantities. For the real qi, q2, pi,p2 ,Ex,Ey ,Hx,Hy we get:

TT Ey

qi = Hx, pi = - —, q2

KQ

Ну

P2

ß ko

Hamiltonian equations [11]:

1 - 2 ( _ — ~2 , 1..„2 tj _ 1 Л. 1 kX „2 , 17„2^j2

V ßko)

H1 = 2k2 ( e " -Û\Pi + й^ь H2 = 0 ( ß - 7 tDq^ + т^ЭД.

2

-1 kl )

2

TE-mode:

&h =k2

dz =k0

d i

(£-1 МЛ

V ßkU

pi,

TM-mode:

d

= - ß (Ii-,

d 2 2 dl = £ koP2-dp2

d

= - ß

( 1kl\ - щ)д2

K

0

The Hamiltonians for TE and TM-mode are separable, thus there are different sim-plectic Runge-Kutta methods for Hamiltonian equation. Using the leapfrog Runge-Kutta method we get the following numerical schemes:

Pk+1 = Pk - hj

К 2

Qk+1 = Qk + tA

h 2

Qk +

(е- 1 ^ V J ко /

k

e--то (Pk + Pk+1 )

(е- 1 -)

Pk+1 =Pk - h{J - щ)

h2

Qk + 2е ko Pk

h

Qk+1 = Qk + 2еко (Pk + Pk+i)

- TE-mode,

TM-mode.

The Figure 2 shows the phase portrait for e = ei(1 + m cos(2-^z)) and m = 0,1, e1 = 1,5, ^ = 1, = 0, Hx(0) = 1. The coordinate z changes from 0

to 10, number of points is N = 2000, and the iteration step is h = 0.005

— TE-mode I

л

л л A A

2 4

Figure 2. Phase portrait drawn for simplectic 2d order Runge-Kutta method

Figure 3. Modulations of q

LU

^ 0

6 8

z

L0

L2

Conclusion

We have reviewed the simplectic Runge-Kutta methods for the case of the separable Hamiltonian H(p, q) = T(p) + V(q). Just for these types of Hamiltonians symplectic methods are well developed.Still there remains an open question about the numerical methods for the Hamiltonians of the general form. It is also interesting to study the classical numerical methods in their symplectic form. The applications often need to obtain a solution with sufficient accuracy. Perhaps the existing numerical schemes can provide an acceptable 2 - form uj preservation error. This can help to avoid the development of new methods in areas where accuracy, given by the classical scheme is sufficient. For the classical scheme there exist already well established and optimized software implementations. The foregoing does not cancel important theoretical significance of symplectic methods and their application in areas where the accuracy of any of the classical methods will not suffice.

References

1. Шутц Б. Геометрические методы математической физики. — Платон, 1995. [Shutc B. Geometricheskie metodih matematicheskoyj fiziki. — Platon, 1995. ]

2. Дубровин Б. А., Новиков С. П., Фоменко А. Т. Современная геометрия. — 1 издание. — Москва: Наука, 1979. [Dubrovin B. A., Novikov S. P., Fomenko A. T. Sovremennaya geometriya. — 1 издание. — Moskva: Nauka, 1979. ]

3. Бахвалов Н. С., Жидков Н. П., Кобельков Г. М. Численные методы. — 6 издание. — Москва: Бином, 2008. — ISBN 978-5-94774-815-4. [Bakhvalov N. S., Zhidkov N. P., Kobeljkov G. M. Chislennihe metodih. — 6 издание. — Moskva: Binom, 2008. — ISBN 978-5-94774-815-4. ]

4. Budd C. J., Piggott M. D. Geometric Integration and its Applications. — 2004.

5. Candy J., Rozmus W. A Symplectic Integration Algorithm for Separable Hamiltonian Functions // Comput. Phys. — 1991. — No 92. — Pp. 230-256.

6. Kinoshita H., Yoshida H., Nakai H. Integrators and their Application to Dynamical Astronomy // Celestial Mechanics and Dynamical Astronomy. — 1991. — No 50. — Pp. 59-71.

7. Forest E., Ruth R. D. Forth-Order Symplectic Integration // Submitted to Physica D. — 1989.

8. Дирак П. А. М. Лекции по теоретической физике. — Ижевск: Регулярная и хаотическая динамика, 2001. — ISBN 5-93972-026-9. [Dirak P. A. M. Lekcii po teoreticheskoyj fizike. — Izhevsk: Regulyarnaya i khaoticheskaya dinamika, 2001. — ISBN 5-93972-026-9. ]

9. Sevastianov L. A. The System of Hamilton Equations for the Modes of the Electromagnetic Field in a Stratified Isotropic Medium (in russian) // Bulletin of Peoples' Friendship University of Russia. Series "Mathematics. Information Sciences. Physics". — 2011. — No 2. — Pp. 169-171.

10. Sevastianov L. A., Kulyabov D. S. The System of Hamilton Equations for Normal Waves of the Electromagnetic Field in a Stratified Anisotropic Medium // The 12th small triangle meeting of theoretical physics. — Stakcin: Institute of Experimental Physics. Slovak Academy of Sciences, 2010. — Pp. 82-86.

11. Gevorkyan M. N., Kulyabov D. S., Sevastyanov L. A. A Study of Impedance, Frequency and Parametric Excitation of Oscillators (in russian) // Bulletin of Peoples Friendship University of Russia. Series "Mathematics. Information Sciences. Physics". — 2009. — Vol. 4. — Pp. 56-62.

УДК 517.958

Симплектические интеграторы и задача распространения

волн в слоистой среде

М. Н. Геворкян, Ю. В. Гладышева

Кафедра систем телекоммуникаций Российский университет дружбы народов ул. Миклуха-Маклая, д. 6, Москва, 117198, Россия

Рассмотрены численные методы, сохраняющие симплектическую структуру гамиль-тоновой системы. Построен гамильтониан для случая распространения электромагнитной волны в стратифицированной среде без источников. Решены уравнения Гамильтона с помощью вариационного метода Рунге—Кутта 2-го порядка.

Ключевые слова: симплектические интеграторы, симплектическая структура, формализм Гамильтона, уравнения Максвелла без источников.