716

Science in China Series G: Physics, Mechanics & Astronomy 2006 Vol.49 No.6 716—728 DOI: 10.1007/s11433-006-2017-8

Algebraic dynamics solutions and algebraic dynamics algorithm for nonlinear ordinary differential equations WANG Shunjin & ZHANG Hua Center of Theoretical Physics, Sichuan University, Chengdu 610064, China Correspondence should be addressed to Wang Shunjin (email: [email protected]) Received September 19, 2005; accepted July 10, 2006

Abstract The problem of preserving fidelity in numerical computation of nonlinear ordinary differential equations is studied in terms of preserving local differential structure and approximating global integration structure of the dynamical system. The ordinary differential equations are lifted to the corresponding partial differential equations in the framework of algebraic dynamics, and a new algorithm——algebraic dynamics algorithm is proposed based on the exact analytical solutions of the ordinary differential equations by the algebraic dynamics method. In the new algorithm, the time evolution of the ordinary differential system is described locally by the time translation operator and globally by the time evolution operator. The exact analytical piece-like solution of the ordinary differential equations is expressd in terms of Taylor series with a local convergent radius, and its finite order truncation leads to the new numerical algorithm with a controllable precision better than Runge Kutta Algorithm and Symplectic Geometric Algorithm. Keywords: exact algebraic dynamics solutions of ordinary differential equations, algebraic dynamics algorithm, preserving fidelity geometrically and dynamically.

1

Algebraic dynamics solutions to ordinary differential equations

1.1

Multivariable ordinary differential equations and their structures[1,2]

The first order ordinary differential equations with n variables Xi and p parameters αμ read X i = Fi (α μ (t ), X i ), (i = 1, 2,...n; μ = 1, 2... p). (1.1)

Consider the autonomous systems where α μ (t ) is independent of time. The velocity www.scichina.com

www.springerlink.com

Algebraic dynamics solutions & algorithm for nonlinear ordinary differential equations

717

fields Fi define the local differential structure of the dynamical system. The existence of local solutions requires that Fi are continuous, differentiable and compact[1]. 1.2

Algebraic dynamics solution to the first order ordinary differential equations[3

―6]

Based on Fi, we introduce the time translation operator Lˆ ( X , ∂ X ), Fi (α , X ) = Lˆ ( X , ∂ X ) X i , (1.2)

∂ ⎞ n ∂ ⎛ . Lˆ ⎜ X , ⎟ = ∑ Fi (α , X ) ∂X i ⎝ ∂X ⎠ n =1

As Fi, Lˆ ( X , ∂ X ) defines the local differential structure of the dynamical system as well. Then eq. (1.1) can be lifted to the partial differential equations: ∂X ( X , t ) ˆ ⎛ ∂ ⎞ X i = i = L⎜ X , ⎟ Xi. ∂t ∂ X⎠ ⎝

(1.3)

The exact analytical solutions with a local convergent radius are ∞

ˆ t n ˆn L ( X 0 ) X 0i = e L ( X 0 ) t X 0i . ! n n =0

X i ( X 0 , t) = ∑

(1.4)

The analytical solution (1.4) of (1.1) is based on Taylor series with a local convergent radius, in terms of which the algebraic dynamics algorithm is designed. ˆ Like quantum mechanics, we introduce the time evolution operator Uˆ [ X , t ] = e L ( X 0 )t 0

and express the solution in terms of this operator: X (t ) = X ( X 0 , t ) = Uˆ [ X 0 , t ] X 0 , X i (t ) = X i ( X 0 , t ) = Uˆ [ X 0 , t ] X 0i ,

(1.5)

where the linear functions of coordinates X 0i play the roles of wave functions. From the above solution one finds that the structure of ordinary differential equations is defined in two folds: The finite global integration structure is described by the time evolution operator Uˆ [ X , t ] ; and the infinitesimal local differential structure is described by the time 0

translation operator Lˆ ( X 0 ), much like the structure of Schrödinger equation in quantum physics. The relation between the two operators is just that of Lie group and Lie algebra.

1.3

Correctness of the algebraic dynamics solutions

The correctness of the algebraic dynamics solution (1.4) to eq. (1.1) can be proved as follows: (i) Proof 1: Based on differential operators. Taking the time derivative of (1.4) or (1.5), one has ∂X i ( X 0 , t )  (1.6) = X i = Uˆ ( X 0 , t ) Lˆ ( X 0 ) X 0i = Uˆ ( X 0 , t ) Fi (α , X 0 ), ∂t where the equation Lˆ ( X ) X = F ( X ) has been used. Substituting the following identities i

i

ˆ )] = ( Lˆ2 X ) , [ Lˆ , ( Lˆn X )] = ( Lˆn +1 X ) , ˆ ) , [ Lˆ , ( LX [ Lˆ ( X ), X ] = ( LX

(1.7a)

718

Science in China Series G: Physics, Mechanics & Astronomy ∞

n

t Uˆ ( X 0 , t ) X 0Uˆ −1 ( X 0 , t ) = ∑ ( Lˆn ( X 0 ) X 0 ) = Uˆ ( X 0 , t ) X 0 = X ( X 0 , t ), n =0 n!

(1.7b)

Uˆ ( X 0 , t ) Fi (α , X 0 ) = Fi (α , X (t ))

(1.7c)

into eq. (1.6) yields X i = Fi (α , X ) which is exactly the ordinary differential equation (1.1). Thus one has succeeded in proving that the solution (1.4) or (1.5) to the partial differential equations (1.3) is also the solution to the ordinary differential equations (1.1), where the initial coordinates { X 0i } are variables of the partial differential equations (1.3) besides of the time variable. (ii) Proof 2: Based on the generalized Liuoville equation. We introduce the Hermite conjugate operator Lˆ+ of Lˆ : ( Lˆ+ )+ = Lˆ and the probability density function ρ ( X , t ) in X-space, which satisfies the partial differential equation: ∂ρ ( X , t ) ˆ+ ρ ( X , t ) = = L ( X ) ρ ( X , t ). ∂t For the Hamiltonian system Lˆ+ = − Lˆ , (1.8) becomes Liuoville equation:

(1.8)

∂ρ ( X , t ) = − Lˆ ( X ) ρ ( X , t ). (1.9) ∂t Therefore, eq. (1.8) is a generalized Liuoville equation. The solution of (1.8) reads ˆ+ ρ ( X , t ) = Uˆ + ( X , t ) ρ ( X , t = 0), Uˆ + ( X , t ) = e L t . (1.10)

ρ ( X , t ) =

Considering the time evolution of the first order moment of X i , X i (t ) = ∫ X i ρ ( X , t )dX = ∫ X iUˆ + ( X , t ) ρ ( X , t = 0)dX ˆ = ∫ (Uˆ ( X , t ) X i ) ρ ( X , t = 0)dX = ∫ (e Lt X i )ρ ( X , t = 0)dX .

(1.11)

For the point-like initial density distribution ρ ( X , t = 0) = δ ( X − X 0 ) , (1.11) turns to be ˆ

X i (t ) = e L ( X 0 )t X i 0 = X i ( X 0 , t ) = X i (t ).

(1.12)

By virtue of (1.7) and (1.12), the equation of motion for X i (t ) is ˆ ˆ X i (t ) = X i (t ) = e L ( X 0 )t Lˆ ( X 0 ) X i 0 = e L ( X 0 )t Fi (α , X 0 ) = Fi (α , X (t )),

(1.13)

which is exactly the equation of motion (1.1). Thus we have finished our second proof that the partial differential solution (1.12) or (1.4) is exactly the solution of the ordinary differential equation (1.13) or (1.1). (iii) Proof 3: Based on the solution of (1,1) in Taylor series of time. For the ordinary differential equations (1.1), if F i (α , X ) are nonsingular, compact, and differentiable to all orders, its solution is unique and can be expressed in terms of Taylor series of time: ∞

∞ n t n (n) t X i (t = 0) = X i 0 + ∑ X i ( n ) (t = 0), n =0 n ! n =1 n !

X i (t ) = ∑ Xi

(n)

⎛ d n X i (t ) ⎞ (t = 0) = ⎜⎜ ⎟⎟ . n ⎝ dt ⎠t = 0

(1.14)

Algebraic dynamics solutions & algorithm for nonlinear ordinary differential equations

719

On the other hand, from (1.1) and (1.2) one has ˆ , X i = X i (1) = Fi (α , X ) = LX i

X i(2) =

∂Fi  ∂Fi X j = Fj = Lˆ2 X i . ∂X j ∂X j

(1.15)

In the following, we shall prove by induction that if X i ( n −1) = Lˆn −1 X i , then X i ( n ) = Lˆn X i . In fact, from X i ( n −1) = Lˆn −1 X i and (1.1), one has dX i ( n −1) ˆn = L Xi. dt Since (1.15) is valid, (1.16) is also valid in general by induction. From (1.1), eq. (1.16) can be rewritten as X i(n) =

X i ( n ) = Fi ( n −1) = Lˆn X i ,

Fi ( n −1) =

d n −1 Fi dt n −1

.

(1.16)

(1.17)

At t = t0 = 0 , X i ( n ) (t = 0) = Lˆn ( X 0 ) X 0i . The Taylor series solutions (1.14) turn to be ∞

∞ n ˆ t n (n) t X i (t = 0) = ∑ Lˆn ( X 0 ) X 0i = L ( X 0 )t X 0i = Uˆ ( X 0 , t ) X 0i . n =0 n ! n =0 n !

X i (t ) = ∑

(1.18)

Thus we have proved that the time series solutions (1.18) of the ordinary differential equations (1.1) and the solutions (1.4) of the partial differential equations (1.3) are the same and that X ( n ) (t = 0) are equal to Lˆn X to all orders. 0i

i

The above results imply that any nonlinear ordinary differential equation of autonomous systems, if Lˆn X exist to all orders, always has piece-like and locally convergent 0i

analytical solutions in terms of Taylor series of time and that these solutions can be obtained by solving the corresponding partial differential equations. In the above sense, we conclude that any nonlinear ordinary differential equation of autonomous systems, if Lˆn X exist to all orders, is always integrable locally, and that the concept of integrabil0i

ity and non-integrability of a nonlinear system is meaningful only globally, while the distinction between them is meaningless locally. 1.4

Examples for algebraic dynamics solutions to ordinary differential equations

For further check of the correctness of the algebraic dynamics solutions to ordinary differential equations, we give some examples which include Hamiltonian and non-Hamiltonian systems. (i) One-dimensional harmonic oscillator: Its Hamiltonian and equations of motion read 1 H = ( p2 + q2 ) , 2 ∂H ∂H = Fq = p , p = − = Fp = − q. (1.19) q = ∂p ∂q ∂H ∂ ∂H ∂ ∂ ∂ ∂ ∂ The time translation operator is Lˆ = − = Fq + Fp = p −q . ∂p ∂q ∂q ∂p ∂q ∂p ∂q ∂p (1.19) can be lifted to the partial differential equations:

720

Science in China Series G: Physics, Mechanics & Astronomy

X = Lˆ X ,

X = ( q, p ) T ,

(1.20)

which have the solutions: ∞

t n ˆn L (q0 , p0 )q0 , n =0 n !

q(q0 , p0 ; t ) = e L0 ( q0 , p0 )t q0 = ∑ ˆ

∞

t n ˆn L (q0 , p0 ) p0 . n =0 n !

p (q0 , p0 ; t ) = e L0 ( q0 , p0 )t p0 = ∑ ˆ

(1.21a) (1.21b)

Lˆn ( X 0 )q0 and Lˆn ( X 0 ) p0 can be calculated as follows: Lˆ2 n ( X 0 )q0 = (−1)n q0 ,

Lˆ2 n +1 ( X 0 )q0 = (−1)n p0 ,

(1.22a)

Lˆ2 n ( X 0 ) p0 = (−1)n p0 ,

Lˆ2 n +1 ( X 0 ) p0 = −(−1) n q0 .

(1.22b)

Inserting (1.22a,b) into (1.21a,b), one obtains the exact solutions: q(q0 , p0 ; t ) = q0 cos t + p0 sin t , p (q0 , p0 ; t ) = p0 cos t − q0 sin t. (ii) A free particle: Its Hamiltonian and equations of motion are 1 H = p2 , 2 ∂H ∂H = Fq = p , p = − = Fp = 0. q = ∂p ∂q

(1.23)

∂ ∂ The time translation operator Lˆ = p − 0 . From the partial differential equations ∂q ∂p one has the solutions. From Lˆ ( X 0 )q0 = p0 , Lˆn ( X 0 )q0 = 0 (n ≥ 2) ,

Lˆn ( X 0 ) p0 = 0 (n ≥ 1) ,

(1.24)

one obtains q(q0 , p0 ; t ) = q0 + p0 t ,

p (q0 , p0 ; t ) = p0 .

(iii) Free fall of a particle: Its Hamiltonian and equations of motion are 1 H = p 2 + gq , 2 ∂H ∂H = Fq = p , p = − = Fp = − g . q = ∂p ∂q

(1.25)

∂ ∂ . Lˆn ( X 0 )q0 and Lˆn ( X 0 ) p0 read The time translation operator Lˆ = p − g ∂q ∂p Lˆ ( X 0 )q0 = p0 ,

Lˆ2 ( X 0 )q0 = − g ,

Lˆ ( X 0 ) p0 = − g ,

Lˆn ( X 0 )q0 = 0 (n ≥ 3),

Lˆn ( X 0 ) p0 = 0 (n ≥ 2).

(1.26a) (1.26b)

The exact solutions are 1 2 gt , p (q0 , p0 ; t ) = p0 − gt . 2 (iv) 1-D harmonic oscillator with a momentum-dependent potential: Its Hamiltonian and equations of motion are q(q0 , p0 ; t ) = q0 + p0t −

Algebraic dynamics solutions & algorithm for nonlinear ordinary differential equations

1 H = ( p 2 + q 2 ) + kqp, 2 ∂H = Fq = p + kq , q = ∂p

p = −

721

∂H = Fp = −(q + kp ). ∂q

(1.27)

∂ ∂ The time translation operator is Lˆ = ( p + kq ) − (q + kp ) , and Lˆn ( X 0 )q0 and ∂q ∂p Lˆn ( X 0 ) p0 read Lˆ2 n ( X 0 )q0 = (k 2 − 1)n q0 ,

Lˆ2 n +1 ( X 0 )q0 = (k 2 − 1)n ( p0 + kq0 ),

(1.28a)

Lˆ2 n ( X 0 ) p0 = (k 2 − 1) n p0 ,

Lˆ2 n +1 ( X 0 ) p0 = −(k 2 − 1) n (q0 + kp0 ).

(1.28b)

The solutions are q(q0 , p0 ; t ) = q0 cos ω t +

p0 + kq0

ω

sin ω t ,

p (q0 , p0 ; t ) = p0 cos ω t −

q0 + kp0

ω

sin ω t ,

where ω = 1 − k 2 . (v) Particle in a periodic potential: Its Hamiltonian and equations of motion are 1 H = p 2 − k cos q , 2 ∂H ∂H = Fq = p , p = − = Fp = − k sin q, (1.29) q = ∂p ∂q ∂ ∂ Lˆ = p − k sin q . Lˆn ( X 0 )q0 and Lˆn ( X 0 ) p0 up to the 5-th order are ∂q ∂p ˆ = − k sin q , Lˆ ( X 0 )q0 = p0 = Lˆ 0 p0 , Lˆ2 ( X 0 )q0 = Lp 0 0 2

k Lˆ3 ( X 0 )q0 = Lˆ2 p0 = − kp0 sin q0 , Lˆ4 ( X 0 )q0 = Lˆ3 p0 = kp0 2 sin q0 + sin 2q0 , 2 (1.30) 5 4 3 2 2 ˆ ˆ L ( X 0 )q0 = L p0 = kp0 cos q0 + k p0 (cos 2q0 − 2sin q0 ), 11 Lˆ5 ( X 0 ) p0 = − kp0 4 sin q0 − k 2 p0 2 sin 2q0 − k 2 sin q0 (cos 2q0 − 2sin 2 q0 ). 2 The analytical solutions can be obtained by solving the Hamilton-Jacobi equation 1) as follows: q(q0 , p0 , t ) = 2sin −1 (mSN [ kt + δ , m]), p (q0 , p0 , t ) = 2m gCN [ kt + δ , m], m=

⎡ ⎤ E0 + k p2 p0 , E0 = 0 − k cos[q0 ], δ = CN −1 ⎢ , m⎥ , 2k 2 ⎢⎣ 2( E0 + g ) ⎥⎦

(1.31)

where SN and CN are the first and second kind of Jacobi elliptic functions, m and δ are functions of ( q0 , p0 ). The following relations are checked by calculations:

1) Zhang H, Wang S J. Analytical solutions to a particle motion in a periodic potential. J Math Phys (in Chinese), 2006 (to be published)

722

Science in China Series G: Physics, Mechanics & Astronomy

⎡ ∂ n q (q0 , p0 , t ) ⎤ ⎡ n ⎤ ˆn p = ⎢ ∂ p (q0 , p0 , t ) ⎥ . , Lˆn q0 = ⎢ L ⎥ 0 n n ∂t ∂t ⎣ ⎦ t =0 ⎣ ⎦ t =0 (vi) Huygens oscillator: Its Hamiltonian and equations of motions are 1 H = p 2 − 2q 2 + q 4 , 2 ∂H ∂H = Fq = p , p = − = Fp = 4q − 4q3 , q = ∂p ∂q

(1.32)

(1.33)

∂ ∂ Lˆ = p − (4q3 − 4q ) . Lˆn ( X 0 )q0 and Lˆn ( X 0 ) p0 up to the 3-rd order are ∂q ∂p Lˆ ( X 0 )q0 = p0 ,

ˆ = 2q − 2q 3 , Lˆ2 ( X 0 )q0 = Lp 0 0 0

2p Lˆ3 ( X 0 )q0 = Lˆ2 p0 = 0 − 2 p0 q02 , 3

(1.34)

2 ⎞ ⎛2 ⎞⎛ 2 Lˆ3 ( X 0 ) p0 = −4 p0 2 q − ⎜ − 2q02 ⎟⎜ − q02 + q3 ⎟ . 3 ⎠ ⎝3 ⎠⎝ 3 The analytical solutions can be obtained by solving the Hamilton-Jacobi equation 1): q(q0 , p0 , t ) =

f DN [ 2 f t + δ , k ],

p (q0 , p0 , t ) = − 2(1 + E0 ) SN [ 2 f t + δ , k ]CN [ 2 f t + δ , k ], f = 1 + 1 + E0 , E0 =

(1.35)

⎡q ⎤ p0 2 1 DN −1 ⎢ 0 , m ⎥ , − 2q0 2 + q0 4 , δ = 2 2f ⎢⎣ f ⎥⎦

where SN , CN , DN are Jacobi elliptic functions; f , δ , E0 are functions of ( q0 , p0 ). The following relations are checked by calculations: ⎡ ∂ n q (q0 , p0 , t ) ⎤ Lˆn q0 = ⎢ ⎥ , ∂t n ⎣ ⎦ t =0

⎡ ∂ n p (q0 , p0 , t ) ⎤ Lˆn p0 = ⎢ ⎥ . ∂t n ⎣ ⎦ t =0

(1.36)

(vii) Arnold’s anti-example I[1]. It is a non-Hamiltonian system: x = v( x) = x 2 / 3 , where v( x) is not differentiable and the solution is not unique at the point x0 = 0. The t3 . Since Lˆn x |x =0 (n = 1, 2,3...∞) exist, the 27 ˆ , Lˆ = x 2 / 3 ∂ , and its solution solutions can be obtained by algebraic dynamics: x = Lx ∂x

two solutions are x1 (t ) = 0 and x2 (t ) =

∞

t n ˆn L x0 . As x0 = 0, the first solution obviously reads x1 (t ) = 0. n =0 n !

is x(t ) = e Lt x0 = ∑ ˆ

ˆ = x 2 / 3 , Lˆ2 x = 2 x 1/ 3 , As to the second solution, a direct calculation yields: Lx 0 0 0 0 3

1) Zhang H, Wang S J. Analytical solutions to Huygens oscillator. J Math Phys (in Chinese), 2006 (to be published)

Algebraic dynamics solutions & algorithm for nonlinear ordinary differential equations

2 Lˆ3 x0 = , 9

723

(c + t )3 Lˆn x0 = 0 (n ≥ 4) and x(t ) = . As c = 0, the second solution thus 27

t3 . Algebraic dynamics method can apply for more general equations: 27 q n −1 x = x q / p , = . p n

reads x2 (t ) =

(viii) Arnold’s anti-example II[1].

x = v( x) = x 2 . Since v( x) is non-compact, the

mapping x(t ) = ϕ (t ) = g t ( x0 ) : R → R is not diff-isomorphism, and there are two pieces of disconnected solutions. Algebraic dynamics solutions are as follows: At the neighborhood of (t = 0, x(t = 0) = x0 ) , the first piece of solutions can be obtained as ˆ , Lˆ = x 2 ∂ , Lˆn x = n ! x n +1 , x = Lx 0 ∂x ∞

∞ x0 t n ˆn L x0 = x0 ∑ ( x0 t )n = 1 − x0 t n =0 n ! n =0

x(t ) = e Lt x0 = ∑ ˆ

(| x0 | t < 1 ).

The above solutions can also be obtained by solving a nonautonomous system. Let 1 dτ dx 1 t = , dt = − 2 , the equation of motion becomes nonautonomous: = − 2 x2 = τ τ d τ τ 1 ∂ Lˆ (τ ) x , Lˆ (τ ) = − 2 x 2 . The solution is ∂x τ τ

τ n−1 τ1 ∞ τ ∫ dτ 'Lˆ (τ ') x(t ) = Tˆe∞ x0 = ∑ ∫ dτ1 ∫ dτ 2 ....... ∫ dτ n Lˆ (τ1 ) Lˆ (τ 2 )....Lˆ (τ n −1 ) Lˆ (τ n ) x0 . n=0 ∞

∞

∞

Tˆ is a time-ordering operator. A direct calculation yields τ

1 2 ∫ dτ1Lˆ (τ1 ) x0 = x0 ,

τ

∞

τ

τ1

τ n−1

∞

∞

∞

τ

τ1

∫ dτ1 ∫ dτ 2 Lˆ (τ1 ) Lˆ (τ 2 )x0 =

∞

∞

1

1

τ2

∫ dτ1 ∫ dτ 2 ..... ∫ dτ n Lˆ (τ1 ) Lˆ (τ 2 ).....Lˆ (τ n )x0 = τ n x0

x03 ,

n +1

,

n

∞ x x0 ⎛x ⎞ x(t ) = x0 ∑ ⎜ 0 ⎟ = 0 = x 1 − x0t n =0 ⎝ τ ⎠ 1− 0

(| x0 | / τ < 1 or | x0 | t < 1 ).

τ

ˆ , For the second piece of solutions ( | x0 | t > 1 ), one sets ( t0 , x(t0 ) = x0 ), x = Lx ∞ (t − t0 )n ˆn ∂ ˆ Lˆ = x 2 , and the solution is x(t ) = e L (t −t0 ) x0 = ∑ L x0 , Lˆn x0 = n ! x0 n +1 , ∂x n ! n =0 ∞

x(t ) = x0 ∑ [ x0 (t − t0 )]n = n =0

x0 , (| x0 | t > 1, | x0 | (t − t0 ) < 1 ). 1 − x0 (t − t0 )

724

1.5

Science in China Series G: Physics, Mechanics & Astronomy

Convergence of Taylor series

The ratio between the (n+1)-th term and the n-th term is ∂ ln( Lˆn X 0i ) Δt Lˆn +1 X 0i Δt N ( ) . rin ( X 0 ) = = F X ∑ j 0 ∂X n + 1 Lˆn X 0i n + 1 j =1 0j

(1.37)

n →∞ → ri ( X 0 ) < 1 yields the convergent radius Ti : The convergent condition rin ( X 0 ) ⎯⎯⎯

n +1

Δt < Ti ,

∂ ln( Lˆn X 0i ) F ( X ) ∑ j 0 ∂X 0j j =1 N

n →∞ ⎯⎯⎯ → Ti .

(1.38)

For the power function Fi ( X 0 ) with its maximum power Fi max ( X 0 ) ~ X 0i m , N

∑ Fj ( X 0 ) j =1

∂ ln( Lˆn X 0i ) ~ X 0i m −1 = const., ∂X 0 j

Ti → ∞.

(1.39)

For the harmonic oscillator, ∂ ln( Lˆn X 0i ) p ⎛ q ⎞ ~ 0 ⎜ 0 ⎟ = const., Ti → ∞. (1.40) ∂X 0 j q0 ⎜⎝ p0 ⎟⎠ j =1 For the numerical calculation of the n-th order truncation, the convergent radius is n +1 Ti ( X 0 , n) = . N ∂ ln( Lˆn X 0i ) ∑ Fj ( X 0 ) ∂X 0j j =1 N

∑ Fj ( X 0 )

1.6

Eigen equation of Lˆ and constants of motion

The equation of motion for any observable is ∂O ˆ = LO , O ( X ) = ∑ X i ∂X i i ˆ

with the solution O[ X (t )] = e L ( X 0 )t O[ X 0 ] . For the constant of motion one has ˆ =0. C ( X ) = LC The eigen equation of Lˆ is ˆ (X ) = λ O (X ) . LO n n n ˆ The set of all the eigen vectors O0c with zero eigen value ( LO 0 c ( X ) = 0 ) span the sub-space of all conserved quantitities of the system. 2

2.1

Truncation of the time evolution operator and algebraic dynamics algorithm

N-th order truncation and N-th order algebraic dynamics algorithm

It has been pointed out that the structure of the dynamical evolution of a nonlinear system is two-fold: the local differential structure described by the time translation op-

Algebraic dynamics solutions & algorithm for nonlinear ordinary differential equations

725

erator Lˆ and the global integration structure described by the time evolution operator Uˆ ( Lˆ Δt ) . To preserve scientific fidelity of a dynamical system in numerical computation, one should preserve the above two structures as accurate as possible. From this view of point, the strategy of algebraic dynamics algothrim is as follows: (1) preserving the local differential structure completely, namely, keeping the integration of Lˆ ; (2) approximating the global integration structure at the N-th order of Δt . This strategy leads to the n

(Δt ) ˆn naive algebraic dynamical algorithm Uˆ N ( Lˆ Δt ) as follows: Uˆ N ( Lˆ Δt ) = ∑ L . n =0 n! N

2.2

Fidelity of Uˆ N (Δ t )

(i) Precision of orbital is (Δt ) N : X i (t + Δt ) = Uˆ N (Δt ) X i (t ) + O(Δt N +1 ) = X i (Δt ) N + O(Δt N +1 ).

(2.1)

(ii) Precision of energy and constants of motion is (Δt ) N : E = H [Uˆ N (Δt ) X (t ) + O (Δt N +1 )] = EN (Δt ) + EOE (Δt N +! ) , EN (Δt ) = E[1 − OE (Δt N +! )] .

(2.2)

(iii) Precision of algebraic relations and Casimirs is (Δt ) N : From X i (q, p ) = X i (q0 , p0 , Δt ) N + Oi (Δt N +1 ) , one has { X i (q0 , p0 , Δt ), X j (q0 , p0 , Δt )}N + Oij (Δt N +1 ) = Cij k [ X k (q0 , p0 , Δt ) N + Ok (Δt N +! )] → { X i (Δt ), X j (Δt )}N = Cij k X k (Δt ) N + [Cij k Ok (Δt N +! ) − Oij (Δt N +! )] ,

(2.3a)

Cn [ X i (Δt )]N = Cn [ X i (0)] + On (Δt N +1 ) .

(2.3b)

(iv) For Hamiltonian systems: (1) Both Uˆ 2 p (Δt ) and Uˆ 2 p +1 (Δt ) preserve symplectic structure at (2p+1) order: Uˆ T2 p ( z , Δt ) JUˆ 2 p ( z , Δt ) = J + O(Δt 2 p + 2 ) , Uˆ T2 p +1 ( z , Δt ) JUˆ 2 p +1 ( z , Δt ) = J + O(Δt 2 p + 2 ). 6

Δt ˆ6 ˆ T (2) Uˆ 4 −6 (Δt ) = Uˆ 4 (Δt ) − LH , U 4−6 ( X , Δt ) JUˆ 4 −6 ( X , Δt ) = J + O(Δt 8 ) . 144 Δt 6 ˆ6 3Δt 8 ˆ8 Uˆ 4 −8 (Δt ) = Uˆ 4 (Δt ) − LH + LH , 144 8 × 144 Uˆ T ( X , Δt ) JUˆ ( X , Δt ) = J + O(Δt10 ). 4 −8

4 −8

(v) Energy correction to raise the precision of energy.

Let E ( X i ) = E0 = const. ,

X i (t + Δt ) = Uˆ N ( X (t ), Δt ) X i (t ) , N

EN (t + Δt ) = E ( X i N (t + Δt )) ,

ΔEN (t ) = EN (t + Δt ) − E0 ,

726

Science in China Series G: Physics, Mechanics & Astronomy

⎛ ∂E ( X i ) ⎞ there is ε i (Δt N +1 ) = ΔEN (t ) ⎜ . ⎟ ⎝ ∂X i ⎠ X N (t +Δt )

The modified orbital is X i N (t + Δt ) = X i N (t + Δt ) − ε i (Δt N +1 ) . 2.3 One parameter algebraic dynamics algorithm and composite algebraic dynamics algorithm The above algorithm is the one parameter algebraic dynamics algorithm. If the system has a composite algebra, one can design a composite algebraic dynamics algorithm as follows: (i) Finding the composite algebraic structure of the system: H (q, p) = H [ X i (q, p )], { X i , X j } = Cij k X k , where { A, B} is Poisson bracket.

(ii) Establishing equations of motion for { X i } : X i = { X i , H [ X }} = Fi [ X ] , where the nonlinear degree of Fi ( X ) with respect to { X i } is lower than that with respect to (q, p). (iii) Lifting equations of motion for { X i } to the partial differential equations of motion: ∂X [ X , t ] ˆ = L[ X ] X i , X i = i ∂t

∂ Lˆ = ∑ Fi [ X ] . X ∂ i i

ˆ (iv) Obtaining the solution: X i [ X 0 , t ] = Uˆ [ X 0 , t ] X i 0 = e L[ X 0 ]t X i 0 .

(v) Designing the composite algebraic dynamics algorithm according to truncation apˆ proximation of Uˆ [ X , Δt ] = e L[ X 0 ]Δt : the N-th order algorithm Uˆ is Uˆ ( Lˆ Δt ) = 0

N

N

n

(Δt ) ˆn L . The precision of Uˆ N is (Δt ) N : ! n n =0 N

∑

X i (t + Δt ) = Uˆ N (Δt ) X i (t ) + O(Δt N +1 ) , { X i N , X Nj } = Cij k X kN + Oij N +! (Δt ) , EN = H (qN , pN ) = E[1 − OE N +1 (Δt )] , Cn [ X i (Δt )] = Cn [ X i (0)][1 − On N +1 (Δt )] , E = H [Uˆ N (Δt ) X (t ) + O (Δt N +1 )] = H [Uˆ N (Δt ) X (t )] + EOE (Δt N +! ) = EN (Δt ) + EOE (Δt N +! ). 3

Summary: Features of algebraic dynamical algorithm

The algebraic dynamics algorithm is one kind of Lie method to the solution of differ－ ential equations. Compared with symplectic geometry method[7 14] and other Lie meth－ ods such as Lie algebra method[15 18], Lie group method[19,20], Lie series and Lie transformation method[21], and formal series and formal dynamics theory[22,23], the algebraic dynamics algorithm has its own features. Algebraic dynamics[3,4] was proposed in 1993 for solving the problem of nonautono-

Algebraic dynamics solutions & algorithm for nonlinear ordinary differential equations

727

mous quantum systems by virtue of Lie algebras and Lie groups, addressing the symmetry and conserved quantities of the system. Then it was applied to solving the evolution problems of classical dynamical systems, and the algebraic dynamics algorithm was established as an approximate numerical method which has the following features. (i) The ordinary differential evolution equations are lifted to the partial differential evolution equations, and the time translation partial differential operator Lˆ is introduced, which describes the local differential structure of the dynamical system and conˆ ) describing the stitutes the infinitesimal generator of the one-parameter Lie group Uˆ ( Lt global integration structure of the dynamical system. (ii) The exact analytical solutions of the dynamical system can be obtained by applyˆ ) on the initial coordinates X (t ) = Uˆ ( Lt ˆ ) X (0) which are piece-like funcing Uˆ ( Lt i

i

tions expressed in terms of Taylor series of t , with local convergent radius and coefficients calculated by partial differential operations. ˆ ) lead to some kinds of algebraic dynamics algorithms (iii) Approximations of Uˆ ( Lt which all preserve the integration of the time translation operator Lˆ ―the local differential structure of the dynamical system completely and approximate the time evolution ˆ ) ―the global integration structure with high precision. There are three operator Uˆ ( Lt kinds of algebraic dynamics algorithms: (1) The naive algebraic dynamics algorithm n

(Δt ) ˆn Uˆ N ( Lˆ Δt ) = ∑ L with precision of (Δt ) n . (2) the symplectic algebraic dynamics n =0 n ! N

algorithm sUˆ N ( Lˆ Δt ) which preserves a local symplectic geometric structure of a Hamiltonian system [7 −14] , with precision of (Δt ) n [24]; (3) the exact algebraic dynamics algorithm eUˆ N ( Lˆ Δt ) with precision of computer floating errors 1). (iv) In the algebraic dynamics solution and algebraic dynamics algorithm, the time evolution problem of classical systems has been cast into the same framework of partial differential evolution equations, and then the useful concepts and methods in quantum physics can be used straightforward. Acknowledgements Valuable discussion with Profs. Guo Hanying, Qin Mengzhao, Tang Yifa, Sun Geng, Wu Ke, Du Mengli, and Ju Changsheng is grateful. This work was supported in part by the National Natural Science Foundation of China (Grant Nos. 10375039 and 90503008), the Doctoral Program Foundation from the Ministry of Education of China, and the Center of Nuclear Physics of HIRFL of China.

References 1 Arnold V I. Ordinary Differntial Equations. Beijing: Science Press, 2001 2 Arnold V I. Geometrical Methods in the Theory of Ordinary Differential Equations. New York: Springer-Verlag, 1988

1) Wang S J, Zhang H. Exact algebraic dynamics algorithm. Internal report of Sichuan University, 2005

728

Science in China Series G: Physics, Mechanics & Astronomy

3 Wang S J, Li F L, Weiguny A. Algebraic dynamics and time-dependent dynamical symmetry of nonautonomous systems. Phys Lett A, 1993, 180: 189 4 Wang S J. Theoretical studies of man-made quantum systems and algebraic dynamics. Progr Phys (in Chinese), 1999, 19(4): 331―370 5 Wang S J. Advanced Quantum Theory and Quantum Many-Body Theory, Chaper 10 (in Chinese). Chengdu: Sichuan University Press, 2005 6 Wang P, Wang S J. Exact analytical solutions to one-dimensional nonautonomous harmonic system. High Energy Phys Nucl Phys (in Chinese), 2005, 29(7): 651―656 7 Feng K. Difference scheme for Hamiltonian formalism and symplectic geometry. J Comput Math, 1986, 4(3): 279―289 8 Feng K, Wu H M, Qin M Z, et al. Construction of canonical difference scheme for Hamiltonian formalism via generating functions. J Comput Math, 1989, 7(1): 71―96 9 Feng K, Qing M Z. Symplectic Geometric Algorithm (in Chinese). Hangzhou: Zhejiang Science and Technology Press, 2003 10 Sun G. A simple way constructing symplectic Runge-Kutta methods. J Comput Math, 2000, 18: 61 11 Sun G. Construction of high order symplectic Runge-Kutta methods. J Comput Math, 1993, b(11): 365 12 Guo H Y, et al. Difference discrete variational principles, Euler―Lagrange cohomology and symplectic, multisymplectic structures (I): Difference discrete variational principle. Commun Theor Phys, 2002, 37: 1 13 Guo H Y, et al. Difference discrete variational principle, Euler―Lagrange cohomology and symplectic, multisymplectic structures (II): Euler―Lagrange cohomology commun. Theor Phys, 2002, 37: 129 14 Guo H Y, et al. Difference discrete variational principles, Euler―Lagrange cohomology and symplectic, multisymplectic structures (III): Application to symplectic and multisymplectic algorithms. Commun Theor Phys, 2002, 37: 257 15 Dragt A J, Finn J M. Lie series and invariant functions for analytic symplectic maps. J Math Phys, 1976, 17: 2215 16 Dragt A J, Forest E. Computation of nonlinear behavior of Hamiltonian systems using Lie algebraic methods. J Math Phys, 1983, 24: 12 17 Healy L M, Dragt A J, Gjaja I M. Computation of error effects in nonlinear Hamiltonian systems using Lie algebraic methods. J Math Phys, 1992, 33: 6 18 Cary J R. Lie transform perturbation theory for Hamiltonian systems. Phys Rep, 1981, 79: 129 19 Iserles A, Munthe-Kaas H, Norsett S P, et al. Acta Numerica. Cambridge: Cambridge University Press, 2000, 215―365 20 Feng K. Formal Dynamical Systems and Numerical Algorithm (Feng Kang’s Selected Papers II)(in Chinese). Beijing: National Defense Industry Press, 1995. 226―234 21 Feng K, Qing M Z. Symplectic Geometric Algorithm, Chapter 9: Formal Power Series. Hangzhou: Zhejiang Science and Technology Press, 2003 22 Olver P J. Applications of Lie Groups to Differential Equations. New York: Springer-Verlag, 1993 23 Steinberg S. Lie Transformations and Their Applications. In: Mondrangon J S, Wolf K B, eds. Lecture Notes in Physics. Berlin: Springer-Verlag, 1985. 45―103 24 Wang S J, Zhang H. Algebraic dynamics algorithm: Numerical comparison with Runge-Kutta algorithm and symplectic geometric algorithm. Sci China Ser G-Phys Mech Astron(in Chinese), 2006, 36(1): 14―37