Integrator for Lagrangian Dynamics

C. Udrişte, M. Postolache and I. Ţevy

Abstract

In §1 we bring the reader gently along from the theory of discrete singletime Lagrangian dynamics (for details see [3]). §2 develops the theory of discrete multi-time Lagrangian dynamics, emphasizing the possibility of computer modelling via Newton method (for details see [8], [9]).

Mathematics Subject Classification: 49KXX, 65P10 Key words: Lagrangian energy, extremals, discrete Euler-Lagrange equations

1 Single-time Lagrangian dynamics

Let x ∈ R^{n} and let γ = {(t, x)|x = x(t), t_{0} ≤ t ≤ t_{1}} be a curve in the space R × R^{n}. Let Γ be the set of all curves γ which join the point (t_{0}, x_{0}), (t_{1}, x_{1}). A differentiable function

L: R × R^{n} × R^{n} → R

is called Lagrangian density energy.

Theorem 1.1 A C^{2} curve γ is an extremal of the functional

E(γ) = ∫_{t_{0}}^{t_{1}} L (t, x, \frac{dx}{dt}) dt, γ ∈ Γ

iff it is a solution of Euler-Lagrange equations

\frac{∂L}{∂x^{k}} − \frac{d}{dt} \frac{∂L}{∂ (\frac{dx^{k}}{dt} )} = 0, k = 1, …, n.

This is a differential system of n equations of second order and therefore the solutions depend on 2n arbitrary constants. For fixing one solution there are used 2n boundary conditions

x(t_{0}) = x_{0}, x(t_{1}) = x_{1}.

Balkan Journal of Geometry and Its Applications, Vol.6, No.2, 2001, pp. 109-115 \begin{matrix}c \\ ◯ \end{matrix}Balkan Society of Geometers, Geometry Balkan Press

Since the property for γ to be an extremal of the functional E does not depend on the choice of the system of coordinates, the previous explanations extend easily to (Riemannian) manifolds.

The Lagrangian L produces the Hamiltonian

H(t, x, ẋ) = ẋ^{i} \frac{∂L}{∂ẋ^{i}} (t, x, ẋ) − L(t, x, ẋ),

which is conserved along the extremals only if it does not depend on the parameter t.

The discretization of L can be made by using the midpoint rule which consists in the substitution of t by h, of the point x with \frac{x_{k+1} + x_{k}}{2} , and of the velocity ẋ by

\frac{x_{k+1} − x_{k}}{h}, where h is the time step. One obtains the discrete Lagrangian

L_{d}: R × R^{n} × R^{n} → R, L_{d}(h, u, v) = L (h, \frac{u + v}{2} , \frac{v − u}{h}) .

This determines the action

S: R × (R^{n})^{N+1} → R, S(h; x_{0}, x_{1}, …, x_{N}) = ∑_{k=0}^{N−1} L(h, x_{k}, x_{k+1}), x_{k} ∈ R^{n}.

The discrete variational principle characterizes the sequence (x_{0}, x_{1}, …, x_{N}) ∈ (R^{n})^{N+1} for which the action S is stationary, for any family x_{k}(ε) ∈ R^{n}, ε ∈ I ⊂ R, 0 ∈ I, with x_{k}(0) = x_{k}, and x_{0}, x_{N} fixed points. Using the variation of the first order we find

Theorem 1.2 The sequence (x_{k}), k = 0, 1, …, N is stationary for the action S iff it is generated by the discrete Euler-Lagrange equations (variational integrator)

\frac{∂L_{d}}{∂x_{k}^{i}} (h, x_{k−1}, x_{k}) + \frac{∂L_{d}}{∂x_{k}^{i}} (h, x_{k}, x_{k+1}) = 0, i = 1, …, n; k = 1, …, N − 1.

Denoting

A_{i}(k) = \frac{∂L_{d}}{∂x_{k}^{i}} (h, x_{k−1}, x_{k}), f_{i}(u) = \frac{∂L_{d}}{∂x_{k}^{i}} (h, x_{k}, u) + A_{i}(k),

i = 1, …, n; u = (u^{1}, …, u^{n}); F = (f_{1}, …, f_{n}),

the preceding equations tansfer into a nonlinear equation system

F (u) = 0, (1)

at each step k. The solution of the system (1) can be approximated by using u(e) = x_{k} in the Newton method

J_{F} (u(e)) ( \begin{matrix}u^{1}(e + 1) \\ ⋯ \\ u^{n}(e + 1) \end{matrix} ) = J_{F} (u(e)) ( \begin{matrix}u^{1}(e) \\ ⋯ \\ u^{n}(e) \end{matrix} ) − ( \begin{matrix}f_{1}(u(e)) \\ ⋯ \\ f_{n}(u(e)) \end{matrix} ) ,

e = 1, 2, …, ē. The matrix J_{F} is the Jacobi matrix of the function F, and the integer number ē is determined by the condition

\frac{|∑_{i=1}^{n} f_{i}^{2}(u(e)) − ∑_{i=1}^{n} f_{i}^{2}(u(e + 1))|}{1 + ∑_{i=1}^{n} f_{i}^{2}(u(e + 1))} < ε.

Next, we put x_{k+1} = u(ē).

The discrete Lagrangian L_{d} produces the discrete Hamiltonian

H_{d}(k) = \frac{x_{k}^{i} − x_{k−1}^{i}}{h} \frac{∂L_{d}}{∂ẋ^{i}} (h, x_{k−1}, x_{k}) − L_{d}(h, x_{k−1}, x_{k}).

2 Multi-time Lagrangian dynamics

In the preceding section the space coordinate and time played quite distinct roles: the space coordinate was merely an index numbering freedom degrees, and the time coordinate was the usual physical time in which the system evolves. According Dickey [1], such a theory is satisfactory unless we turn our attention to relativistic invariant equations, e.g. chiral fields, sine-Gordon, and others. Also considering the KP-hierarchy for arbitrary m and n, the variables x_{m} and x_{n} are quite equal in rights and there is no reason to prefer one to the other by choosing it as time. In such cases a new field theory is useful which involves many time variables. The multi-time formalism is of interest even for all examples where the variables are involved in a distinctly asymmetric way.

The multi-time calculus of variations is of course not new. Let T_{0} ⊂ R^{p} be a relatively compact domain. Let φ = {(t, x)|x = x(t), t ∈ T_{0}, x ∈ R^{n}} be a parametrized sheet in the space R^{p} × R^{n}. Let Φ be the set of all parametrized sheets φ satisfying the boundary condition φ|_{∂T_{0}} = f, where f is a given continuous function. A differentiable function

L: R^{p} × R^{n} × R^{np} → R

is called Lagrangian density energy.

Theorem 2.1 A C^{2} parametrized sheet is an extremal of the functional

E(φ) = ∫_{T_{0}} L (t^{α}, x^{i}, \frac{∂x^{i}}{∂t^{α}} ) dt^{1} ∧ ⋯ ∧ dt^{p}

iff it is a solution of Euler-Lagrange equations

\frac{∂L}{∂x^{k}} − \frac{∂}{∂t^{α}} \frac{∂L}{∂x_{α}^{k}} = 0, x_{α}^{k} = \frac{∂x^{k}}{∂t^{α}} , k = 1, …, n.

This is a PDEs system with n equations of second order. Therefore the solutions depend on 2n arbitrary functions. For fixing one solution we use the boundary condition.

Since the property of φ to be an extremal of the functional E does not depend on the choice of the system of coordinates, the previous explanations extend easily to Riemannian manifolds.

The Lagrangian L produces the Hamiltonian

H(t^{α}, x^{i}, x_{α}^{i}) = x_{β}^{j} \frac{∂L}{∂x_{β}^{j}} (t^{α}, x^{i}, x_{α}^{i}) − L(t^{α}, x^{i}, x_{α}^{i}),

but this function is not conserved along the extremals.

Suppose p = 2. The discretization of L can be made by using the centroid rule which consists in the substitution of the point t with the step (h_{1}, h_{2}), of the point

x with \frac{x_{kℓ} + x_{k+1ℓ} + x_{kℓ+1}}{3} , and of partial velocities x_{α}, α = 1, 2, by \frac{x_{k+1ℓ} − x_{kℓ}}{h_{1}} ,

\frac{x_{kℓ+1} − x_{kℓ}}{h_{2}} . One obtains the discrete Lagrangian

L_{d}: R^{2} × R^{n} × R^{n} × R^{n} → R, L_{d}(u, v, w) = L (h_{1}, h_{2}, \frac{u + v + w}{3} , \frac{v − u}{h_{1}}, \frac{w − u}{h_{2}}) .

This determines the 2-dimensional discrete action

S: R^{2} × (R^{n})^{(M+1)(N+1)} → R,

S(h_{1}, h_{2}; A) = ∑_{k=0}^{M−1} ∑_{ℓ=0}^{N−1} L(h_{1}, h_{2}; x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1}), x_{kℓ} ∈ (R^{n})^{(M+1)(N+1)},

A = ( \begin{matrix}x_{00} & x_{01} & ⋯ & x_{0N} \\ x_{10} & x_{11} & ⋯ & x_{1N} \\ ⋯ \\ x_{M0} & x_{M2} & ⋯ & x_{MN} \end{matrix} ) .

We fix the 2-step (h_{1}, h_{2}). The discrete variational principle consists in the characterization of the matrix A for which the action S is stationary, for any family x_{kℓ}(ε) ∈ R^{n}, k ∈ {0, 1, …, M − 1}, ℓ ∈ {0, 1, …, N − 1}, ε ∈ I ⊂ R, 0 ∈ I with x_{kℓ}(0) = x_{kℓ}, and fixed lines (x_{00}, x_{01}, …, x_{0N}), (x_{M0}, x_{M2}, …, x_{MN}), fixed columns t(x_{00}, x_{10}, …, x_{M0}), ^{t}(x_{0N}, x_{1N}, …, x_{MN}). The discrete variational principle is obtained using the variation of S of the first order.

Theorem 2.2 The first variation of the discrete action S is

δS(A)(η) = \frac{∂L}{∂x_{00}^{i}} (x_{00}, x_{10}, x_{01})η_{00}^{i} + \frac{∂L}{∂x_{10}^{i}} (x_{10}, x_{20}, x_{11})η_{10}^{i} + \frac{∂L}{∂x_{01}^{i}} (x_{01}, x_{11}, x_{02})η_{01}^{i}

+ ∑_{k=1}^{M_{1}} ∑_{ℓ=1}^{N−1} [\frac{∂L}{∂x_{kℓ}^{i}} (x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1}) + \frac{∂L}{∂x_{kℓ}^{i}} (x_{k−1ℓ}, x_{kℓ}, x_{k−1ℓ+1})

+ \frac{∂L}{∂x_{kℓ^{i}}} (x_{kℓ−1}, x_{k+1ℓ−1}, x_{kℓ})] η_{kℓ}^{i} + ∑_{k=1}^{M−1} \frac{∂L}{∂x_{k0}^{i}} (x_{k−10}, x_{k0}, x_{k−11})η_{k0}^{i}

+ ∑_{ℓ=0}^{N−1} \frac{∂L}{∂x_{Mℓ}^{i}} (x_{M−1ℓ}, x_{Mℓ}, x_{M−1ℓ+1})η_{Mℓ}^{i} + ∑_{ℓ=1}^{N−1} \frac{∂L}{∂x_{0ℓ}^{i}} (x_{0ℓ−1}, x_{1ℓ−1}, x_{0ℓ})η_{0ℓ}^{i}

+ ∑_{k=0}^{M−1} \frac{∂L}{∂x_{kN}^{i}} (x_{kN−1}, x_{k+1N−1}, x_{kN})η_{kN}^{i},

where

η = ( \begin{matrix}η_{00} & η_{01} & ⋯ & η_{0N} \\ η_{10} & η_{11} & ⋯ & η_{1N} \\ ⋯ \\ η_{M0} & η_{M1} & ⋯ & η_{MN} \end{matrix} ) , η_{kℓ} = \frac{∂x_{kℓ}}{∂ε} |_{ε=0} ,

x_{kℓ}(ε) ∈ Q, ε ∈ I, 0 ∈ I, x_{kℓ}(0) = x_{kℓ}.

Proof. We use the family of curves x_{kℓ}(ε) ∈ R^{n}, ε ∈ I. Then

S(A(ε)) = ∑_{k=0}^{M−1} ∑_{ℓ=0}^{N−1} L(x_{kℓ}(ε), x_{k+1ℓ}(ε), x_{kℓ+1}(ε)).

We find

δS(A)(η) = \frac{∂}{∂ε}S(A(ε))|_{ε=0} = ∑_{k=0}^{M−1} ∑_{ℓ=0}^{N−1} \frac{∂L}{∂x_{kℓ}^{i}} (x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1})η_{kℓ}^{i}

+ ∑_{k=0}^{M−1} ∑_{ℓ=0}^{N−1} \frac{∂L}{∂x_{k+1ℓ}^{i}} (x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1})η_{k+1ℓ}^{i}

+ ∑_{k=0}^{M−1} ∑_{ℓ=0}^{N−1} \frac{∂L}{∂x_{kℓ}^{i}} (x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1})η_{kℓ+1}^{i}.

Corollary 2.1 (discrete Euler-Lagrange equations). The matrix A = (x_{kℓ}) is stationary for the action S iff

\begin{matrix}\frac{∂L}{∂x_{kℓ}^{i}} (x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1}) + \frac{∂L}{x_{kℓ}^{i}} (x_{k−1ℓ}, x_{kℓ}, x_{k−1ℓ+1}) + \frac{∂L}{∂x_{kℓ}^{i}} (x_{kℓ−1}, x_{k+1ℓ−1}, x_{kℓ}) = 0, \\ i = 1, …, n; k ∈ \{1, …, M − 1\}, ℓ ∈ \{1, …, N − 1\}. \end{matrix}

Proof. Since the boundary of the grid is fixed, the lines

(η_{00}, η_{01}, …, η_{0N}), (η_{M0}, η_{M1}, …, η_{MN}),

and the columns

t(η_{00}, η_{10}, …, η_{M0}), ^{t}(η_{0N}, η_{1N}, …, η_{MN})

must be zero, and η_{kℓ}, k ∈ {1, …, M − 1}, ℓ ∈ {1, …, N − 1} are arbitrary.

The variational integrator described by the discrete Euler-Lagrange equations works as follows:

- we give the lines (x_{00}, x_{01}, …, x_{0N}), (x_{10}, x_{11}, …, x_{1N});

- we denote

u = x_{kℓ+1},

A_{i}(kℓ) = \frac{∂L}{∂x_{kℓ}^{i}} (x_{k−1ℓ}, x_{kℓ}, x_{k−1ℓ+1}), B_{i}(kℓ) = \frac{∂L}{∂x_{kℓ}^{i}} (x_{kℓ−1}, x_{k+1ℓ−1}, x_{kℓ}),

f_{i}(u) = \frac{∂L}{∂x_{kℓ}^{i}} (x_{kℓ}, x_{k+1ℓ}, u) + A_{i}(kℓ) + B_{i}(kℓ), F = (f_{1}, …, f_{n});

- we solve the nonlinear system (4) F (u) = 0, at each step (k, ℓ) using six points of starting as shown a part of the grid,

\begin{matrix}• & •x_{k−1ℓ} & •x_{k−1ℓ+1} \\ •x_{kℓ−1} & •x_{kℓ} & ∗u = x_{kℓ+1} \\ •x_{k+1ℓ−1} & •x_{k+1ℓ} & ∘ \end{matrix}

Newton Method

The solution of the nonlinear system (4) can be approximated by using u(ℓ) = x_{kℓ+1} in the Newton method

J_{F} (u(e)) ( u^{1}(e + 1) ; ⋮ ; u^{n}(e + 1) ) = J_{F} (u(e)) ( u^{1}(e) ; ⋮ ; u^{n}(e) ) − ( f_{1}(u(e)) ; ⋮ ; f_{n}(u(e)) ) ,

e = 1, 2, …, ē. The matrix J_{F} is the Jacobi matrix of the function F, and the integer ē is determined by the condition

\frac{|∑_{i=1}^{n} f_{i}^{2}(u(e)) − ∑_{i=1}^{2} f_{i}^{2}(u(e + 1))|}{1 + ∑_{i=1}^{n} f_{i}^{2}(u(e + 1))} < ε.

Next, we put x_{kℓ+2} = u(ē).

The discrete Lagrangian L_{d} produces the discrete Hamiltonian

H_{d}(k, ℓ) = \frac{x_{k+1ℓ}^{i} − x_{kℓ}^{i}}{h_{1}} \frac{∂L}{∂x_{1}^{i}} (h_{1}, h_{2}, x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1})
+ \frac{x_{kℓ+1}^{i} − x_{kℓ}^{i}}{h_{2}} \frac{∂L}{∂x_{2}^{i}} (h_{1}, h_{2}, x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1}) − L_{d} (h_{1}, h_{2}, x_{kℓ}, x_{k+1ℓ}, x_{kℓ+1}) .

Acknowledgements. A version of this paper was presented by the first two authors at the third Conference of Balkan Society of Geometers, Workshop on Electromagnetic Flows and Dynamics, July 31-August 3, 2000, University POLITEHNICA of Bucharest, Romania.

References

[1] L.A. Dickey: Soliton equations and Hamiltonian systems, Advanced Series in Mathematical Physics, Vol. 12, World Scientific Publishing, Singapore, 1991.

[2] H.Y. Guo, Y.Q. Li, K. Wu and S.K. Wang: Difference discrete variational principle, Euler-Lagrange cohomology and symplectic, multisymplectic structures, arXiv:math-ph/0106001 v1 1 Jun 2001.

[3] D. Opriş and C. Giulvezan: Geometric Integrators and Applications, Mirton Publishing House, Timişoara, 1999 (in Romanian).

[4] C. Udrişte: Electromagnetic dynamical systems and new variants of Lorentz law, Research Grant, University Politehnica of Bucharest-Ministry of National Education, 1998-2001.

[5] C. Udrişte and M. Postolache: Least squares problem and geometric dynamics, Ital. J. Pure Appl. Math., 15 (2000) (in printing).

[6] C. Udrişte and M. Postolache: Least squares problem and magnetic dynamics, in Global Analysis, Differential Geometry and Lie Algebras (Gr. Tsagas (Ed.)); Workshop on Global Analysis, Differential Geometry and Lie Algebras, Thessaloniki, June 24-28 (2000); pp. 169-238, 2000.

[7] C. Udrişte: Geometric Dynamics, Southeast Asian Bulletin of Mathematics, Springer-Verlag 24 (2000), 313-322; Kluwer Academic Publishers, Dordrecht/Boston/London, Mathematics and Its Applications, 513, 2000.

[8] C. Udrişte: Multi-time dynamics dnduced by 1-forms and metrics, in Global Analysis, Differential Geometry and Lie Algebras, BSG Proceedings 5, Ed. Grigorios Tsagas, Geometry Balkan Press (2001), 169-178.

[9] C. Udrişte and M. Postolache: Atlas of Magnetic Dynamics, Geometry Balkan Press, 2001.

Politehnica University of Bucharest Department of Mathematics I Splaiul Independenţei 313, 77206 Bucharest, Romania E-mail:udriste@mathem.pub.ro