- Email: [email protected]
Explicit symplectic partitioned Runge–Kutta–Nyström methods for non-autonomous dynamics
Applied Numerical Mathematics 61 (2011) 832–843
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Explicit symplectic partitioned Runge–Kutta–Nyström methods for non-autonomous dynamics ✩ Fasma Diele, Carmela Marangi ∗ Istituto per le Applicazioni del Calcolo M. Picone, CNR, Via Amendola 122, 70126 Bari, Italy
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 23 June 2010 Received in revised form 26 January 2011 Accepted 5 February 2011 Available online 18 February 2011 Keywords: Symplectic partitioned Runge–Kutta methods Order analysis Nyström methods
We consider explicit symplectic partitioned Runge–Kutta (ESPRK) methods for the numerical integration of non-autonomous dynamical systems. It is known that, in general, the accuracy of a numerical method can diminish considerably whenever an explicit time dependence enters the differential equations and the order reduction can depend on the way the time is treated. In the present paper, we demonstrate that explicit symplectic partitioned Runge–Kutta–Nyström (ESPRKN) methods specifically designed for d ( M −1 q˙ ) = f (q, t ), undergo an order reduction when second order differential equations dt M = M (t ), independently of the way the time is approximated. Furthermore, by means of symmetric quadrature formulae of appropriate order, we propose a different but still equivalent formulation of the original non-autonomous problem that treats the time as two added coordinates of an enlarged differential system. In so doing, the order reduction is avoided as confirmed by the presented numerical tests. © 2011 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction Second order differential equations and their numerical solution play a crucial role in a variety of applications in science and engineering. In this paper we focus on the numerical solution of initial value problems for the special second order equations
d dt
M −1 (t )˙q = f (q, t ),
t ∈ [0, T ]
(1)
where the evolution of the variable q = q(t ) ∈ Rd depends on time by means of the time dependent, non-singular, continuously differentiable matrix M : [0, T ] → Rd×d and through the function f : Rd × [0, T ] → Rd sufficiently smooth in order to assure the well-posedness of the problem. By setting p = M −1 (t )˙q the problem can be recast into the following first order differential system
q˙ = M (t ) p ,
p˙ = f (q, t ).
(2)
Both Eqs. (1) and (2) allow to analyze dynamics relevant to nonlinear science whenever time dependent damped or driven harmonic oscillators are considered. Several examples can be found in electronics, in optics and in classical mechanics literature. Single and coupled oscillators with friction (see [10,2,12]) as well as variable mass dynamics (see [16]), can be described by means of the equations above. Moreover, in quantum mechanics they can describe, after a spatial discretization, the evolution of a time dependent Schrödinger equation related to quantum damped oscillators (see [15,21]). ✩
*
This work has been supported in part by National Research Council grant CNR-RSTL id.332/2008. Corresponding author. E-mail addresses: [email protected] (F. Diele), [email protected] (C. Marangi).
0168-9274/$30.00 © 2011 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2011.02.003
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
833
In case of a time independent matrix M (t ) = M 0 , it is easy to see that the original system (1) can be written as an autonomous second order equation Q¨ = F ( Q ) by treating the time in f (q, t ) as an added coordinate q0 = t, setting Q = (q0 , q), and denoting F = (0, M 0 f ( Q )). The corresponding first order differential system, obtained by setting p = q˙ and P = (1, p ) can be written as
Q˙ = P ,
P˙ = F ( Q ).
(3)
A large number of numerical methods for the integration of such equations can be found in literature. They exploit the two main features of such a kind of problems, that are the separability of the variables Q and P , as well as the linear dependence on P . In particular, powerful methods were provided in [7], where an optimal efficiency was obtained through the minimization of the effective error. Unfortunately, in case of a time dependent matrix M (t ), the original non-autonomous problem (1) cannot be reformulated so to recover the useful algebraic structure of (3). In this case, the classical way to face the problem is to further enlarge the system by introducing another time variable p 0 to account for the explicit time dependence in M (t ). By this trick we can gain the separability but we lose the linear dependence on one of the two sets of variables describing the enlarged system. As a consequence, in the Taylor series expansion of the exact solution, some extra elementary differentials appear which determine further order conditions that are not satisfied by an ESPRKN method. Different numerical approaches have been developed in literature in order to design effective numerical methods for non-autonomous dynamics (2) (see [3,4]). Notice that all the methods referred above are symplectic methods. It is known that, for a given problem of a given algebraic structure, symplectic methods turns out to be the most accurate ones with respect to any other non-symplectic PRK method (see e.g. [14]). As a matter of fact, the problems we are dealing with in the present paper, can all be described in an hamiltonian framework by introducing the function
H ( p , q, t ) =
1 2
p T M (t ) p + V (q, t )
where M (t ) is taken as a positive definite symmetric matrix and −∇q V (q, t ) = f (q, t ). In the special case of autonomous systems (3), the Hamiltonian function is conserved and symplectic methods outperform any other method in preserving the qualitative features of the theoretical flow. Whenever non-autonomous systems (2) are considered the corresponding Hamiltonian system loses some of its features as, for example, the conservation of the energy. In this context it is reasonable to wonder about the real advantage of using symplectic methods. The answer lies in the P -series order analysis of PRK methods where the graphical representation of elementary differentials as bicoloured trees is adopted. Indeed, by means of symplecticness one is able to prove useful relations among trees’ elementary weights useful to reduce the number of conditions at a given order [17]. So, even in the context of non-autonomous systems, symplecticness plays a major role in the reduction of the order conditions. In the recent work [6], developed in the framework of splitting and composition techniques, the authors proposed an optimal strategy to deal with the explicit time dependence. Such a strategy, based on the introduction of a single time variable q0 , loses the separability, but regains the linear dependence. Then, starting from those results, one may wonder if there is any special advantage in the linearity on P with respect to separability, that can be exploited to build optimal PRK methods for this class of problems. In order to provide an answer, in this paper we compare the classical approach, where two time coordinates, q0 and p 0 are added to the system, to an alternative one that makes use of a single time variable q0 as in [6]. As a matter of fact, thanks to the symplecticness, in our paper we are able to provide new results that allow us to establish the number of order conditions for the non-separable problem considered. Unfortunately, we find out that the advantages obtained form the linear dependence on P are the same one gets from the separability of the two sets of variables, i.e. the same number of order conditions is obtained in both cases. As a main result of the paper we propose a new strategy to approximate the time, based on the introduction of two time coordinates, as in the classical approach, and the use of symmetric quadrature formulae. We prove that an ESPRKN method applied on a new (still equivalent) formulation of the problem (3), preserves its accuracy avoiding order reduction. Notice that the proposed strategy for the integral approximation is based on symmetric quadrature formulae in order to preserve the symmetry of the resulting ESPRK method. We recall that, for symmetric methods, a 2nth order method can be obtained by satisfying the order conditions for the (2n − 1)th order. In details, in Section 2, the integration of the explicit time dependence in the framework of explicit symplectic partitioned Runge–Kutta methods (ESPRK) is addressed; in Section 3 we perform the P -series analysis of an ESPRK method applied to the non-separable problem obtained by introducing the only variable q0 and we establish the order conditions of a method up to the fifth order. In Section 4 we suggest how to modify the original problem in order to gain the desired order whenever an ESPRKN method is used to approximate the solutions of the non-autonomous system (2). Theoretical results are validated in Section 5 on the solution of perturbed time dependent oscillators as well as on time dependent linear Schrödinger equations chosen as test problems. Finally, in Section 6 we draw our conclusions.
834
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
2. Classical implementation of ESPRK methods for general non-autonomous dynamics ( A)
We recall that a partitioned symplectic k-stages Runge–Kutta method is defined by two RK methods ( A i j , b i
( B i j , b(i B ) , c (i B ) ), which satisfy the symplecticity conditions (see e.g. [14]): ( A)
bi
(B)
(B)
( A) (B )
B i j + b j A i j = bi
bj ,
( A)
, c i ) and
i , j = 1, . . . , k ,
(4)
( A)
i = 1, . . . , k , (5) ( A) (B) where c i = = j B i j are the nodes of the two Runge–Kutta methods. An important feature of PRK methods j A i j and c i bi
= bi ,
is the possibility of having explicit methods for non-autonomous dynamics
q˙ = g ( p , t ),
p˙ = f (q, t ).
(6)
I.e. dynamics which reduce to separable ones whenever there is no explicit time dependence in one or both the equations describing the q’s and p’s evolution. ( A) (B) 0. The These methods are based on two diagonally implicit methods satisfying A ii · B ii = 0 with b i = 0 and b i = classical form of a k = 2s stage method is given in literature by the following Butcher arrays
0
β1 β1 β1 β1 .. . β1 β1
α1 α1 0 α1 β2 α2 α1 β2 α2 0 .. .
.. . α1 β2 α1 β2
.. .
.. . α2 · · · α2 · · ·
..
. βs αs βs αs
β1 β1 β1 β1 .. . β1 β1 β1
0
α1 β2 α1 β2 0 .. .
.. . α1 β2 α1 β2 α1 β2
.. . ··· ··· ···
..
(7)
.
α s −1 β s α s −1 β s 0 α s −1 β s α s
which define the A (left) and B (right) RK methods with αi = 0 and βi = 0 for i = 1, . . . , s. It is easy to verify that the partitioned method (7) satisfies both symplecticity conditions (4)–(5). For non-autonomous dynamics (6), the application of an explicit symplectic Runge–Kutta method, defined by (7), to the partitioned system
q˙ 0 = 1, q˙ = g ( p 0 , p ),
p˙ 0 = 1, p˙ = f (q0 , q)
where (q0 , q T ) ∈ Rd+1 are solved by the A method and ( p 0 , p T ) ∈ Rd+1 variables are solved by the B method, leads to the following steps
P 0 = pn ,
Q 0 = qn ,
Do i = 1, s + 1
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i
(b)
Q i = Q i −1 + hai g P i , tn + hc i enddo p n +1 = P s +1 , where, by setting (b)
and c i
=
i
qn+1 = Q s+1 ,
t n +1 = t n + h , (a)
α0 = βs+1 = 0, we define bi = αi−1 + βi , for i = 1, . . . , s + 1, ai = αi + βi , for i = 1, . . . , s, c i =
k=1 b i .
Moreover, we set as+1 = 0. (a)
Notice that the nodes of the Runge–Kutta formula (7) are different from the coefficients c i (a)
( A)
(a)
( A) = c 2i −2
( A) = c 2i −1 for i
(b)
and c i
(8) i −1 k=1 ai
for i = 1, . . . , s + 1 (b)
(B)
(B)
which appear in (9). They are related as follows: c 1 = c 1 = 0, c i = 2, . . . , s + 1 and c i = c 2i −1 = c 2i for i = 1, . . . , s + 1. This is the classical implementation of non-autonomous problems: the time acts as two additional coordinates treated as p and q respectively. This way the original non-separable problem can be treated as a separable one in an enlarged space. However, the classical partitioning is simply one of the possible forms in which an ESPRK method can be rewritten in case of non-autonomous dynamics. One may as well treat the time as a single additional coordinated evolving with q in both equations of (6) as in [6], where this particular coupling of the time variable is what provides the optimal performance for the splitting and composition methods when applied to non-autonomous problems (1). Since the introduction of a single time variable leads to a non-separable problem, it is not so trivial to infer the performance of an ESPRK method for the case at hand. We will see that, despite what can be proved in the case of splitting and composition methods, this approach does
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
835
not provide advantages with respect to the classical implementation. On the other hand, we are able to demonstrate that the number of order conditions does not increase with respect to the classical partitioning, even if we are dealing with a non-separable problem. 3. A different implementation and relative order analysis Eqs. (2) belongs to the class of non-autonomous dynamics of type (6) where g ( p , t ) = M (t ) p. We search for approximate solution obtained by applying an ESPRK method defined by (7). As we noticed before, different implementations are possible according to how the time variable is approximated. The classical implementation, described in the previous section is provided by the following inner stages:
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i
(b)
Q i = Q i −1 + hai M tn + hc i
,
Pi.
(9)
The implementation we focus on, instead, corresponds to
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i
(a)
Q i = Q i −1 + h βi M tn + hc i
a) + αi M tn + hc i(+ Pi. 1
(10)
The above method describes the application of an ESPRK to the partitioned system
q˙ 0 = 1, q˙ = M (q0 ) p ,
p˙ = f (q0 , q)
(11)
where the time variable is integrated together with the q coordinate in both equations of (6). To perform the order analysis of the partitioned Runge–Kutta (PRK) method (10) applied to the differential system (11), we follow the classical procedure as given in [14]. Notice that, the application of the theory in this context is not straightforward. In order to get the minimal set of order conditions, we will extend the results given in [9], concerning separable systems, to the non-separable problem at hand and provide new relations that allow to properly reduce the number of necessary order conditions. The order of a PRK method can be characterized in terms of its coefficients via bi-coloured rooted trees. We represent with a ‘white’ vertex the elementary tree representing a f = f (q0 , q), and with a ‘black’ vertex the elementary tree representing a g = M (q0 ) p. The rules for constructing trees, for the case at hand, are the following: - ◦ vertices can have only • children. - • vertices can have both ◦ and • children; ◦ children cannot be more than one, • children are elementary trees. As can be seen from the second rule, the introduction of time affects the separability of the problem and is responsible for the peculiar constraint on black children of black vertices. So, unless we are able to prove that such conditions are enough to drastically reduce the number of trees, we should expect a number of order conditions similar to the one that features an ESPRK method when applied to a generic non-separable problem. In order to apply the order theory within the new set of rules, we denote by T P b and T P w all the bi-coloured trees with black and white roots, respectively. For symplectic PRK methods, from condition (4) and for all u ∈ T P b and v ∈ T P w , it follows that [1]
Φ(u ◦ v ) + Φ( v ◦ u ) = Φ(u )Φ( v )
(12)
where Φ(τ ) denotes the elementary weight associated to the tree τ and ◦ corresponds to the Butcher’s product of trees (see [8]). Moreover, because of the symplecticity property (5), Φ(τ ) is independent of the type of the root of τ . More in details, define u = [τ1 , τ2 , . . . , τn ]b a generic tree in T P b obtained by connecting the roots of τ1 , τ2 , . . . , τn to a new black root and, similarly, v = [σ1 , σ2 , . . . , σl ] w a generic tree in T P w obtained by connecting the roots of σ1 , σ2 , . . . , σl to a new white root. Moreover, suppose that v u = [σ1 , σ2 , . . . , σl ]b ∈ T P b , and u v = [τ1 , τ2 , . . . , τn ] w ∈ T P w , where v u and u v are obtained from v and u by changing the colour of their root. Then, the following relation holds
Φ( v ) = Φ( v u )Φ(u ) = Φ(u v ).
(13)
A direct consequence of (12) and (13) is the relation
Φ(u ◦ v ) + Φ( v u ◦ u ) = Φ(u )Φ( v u )
(14)
which leads to the following further simplifying assumptions which generalize the results given in [9]. Proposition 3.1. Suppose u , ω1 , ω2 , . . . , ωn , ∈ T P b , v ∈ T P w . Let z be the elementary white tree or z = [ω1 , ω2 , . . . , ωn ] w ∈ T P w . Then, the following relation holds
836
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
Table 1 Nyström trees, elementary weights and density function values up to order five.
τ ◦
Φ(τ ) s+1
bi =
i =1
s
γ (τ )
ai
τ
s+1
1
i =1
s+1
Φ(τ )
γ (τ ) (a)
bi c i
2
i =1
(a) 2
bi c i
i −1 s+1
3
i =1
(b)
bi a j c j
6
i =1 j =1
s+1
(a) 3
bi c i
s+1 i −1 i
4
i =1
(a)
b i a j bk c k
24
i =1 j =1 k=1
s+1
(a) 4
bi c i
j l−1 s+1 i −1
5
i =1
(b)
b i a j bl ak ck
120
i =1 j =1 l=1 k=1
j s+1 i −1
(a)
(a)
b i c i a j bl c l
j s+1 i −1
30
i =1 j =1 l=1
(a) 2
b i a j bl c l
60
i =1 j =1 l=1
Φ z ◦ (u ◦ v ) + Φ z ◦ ( v u ◦ u ) = Φ ( z ◦ u ) ◦ v u
(15)
whenever v u and ( v u ◦ u ) ∈ T P b . Proposition 3.2. Suppose u , v , r ∈ T P b that verify Φ(u ) + Φ( v ) = Φ(r ). Let z be the elementary white tree or z = [ω1 , ω2 , . . . ,
ωn ] w ∈ T P w . Then, the following relation holds Φ( z ◦ u ) + Φ( z ◦ v ) = Φ( z ◦ r ).
(16)
For sake of clarity, we refer the reader to Appendix A for the proofs. By applying these results to our problem we find out that, with respect to the trees which assure the fifth order accuracy to a ESPRKN method (reported in Table 1 together with the elementary weights as well as their density functions), new trees appears and are given in Table 2 (top). In the same table (bottom), the (well-known) additional conditions related to the classical separable method (9) are reported for comparison. Both approaches share the same total number of order conditions i.e. 14, even if they are related to different trees, with different weights. Notice that here we are dealing with a non-separable method (10). Remembering that in the case of a generic non-separable problem we may have up to 40 order conditions [17] for symplectic methods, we proved that the method here analyzed does not perform worse than the classical separable one (9). On the other hand, we observe that the advantages one can get within the splitting and composition framework, cannot be recovered by applying an ESPRKN method as in (10). 4. The method of exact integration and its approximations When the coefficients ai and b i of an ESPRKN method are used for the evaluation of
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i
1 Q i = Q i −1 + hai
(a)
M tn + hc i
,
+ hai s ds P i ,
0
the resulting approximations
p n +1 = p n + h
s +1
(a)
b i f Q i −1 , tn + hc i
,
i =1
qn+1 = qn + h
s i =1
1
(a)
M tn + h c i
ai
+ ai s ds P i ,
(17)
0
inherit the accuracy of the ESPRKN method [6]. Indeed, that the splitting underlying the method (17), satisfies splitting of autonomous second order equations (3). This same of classical ESPRKN methods. In particular, the 11
by performing the order analysis via Lie formalism, it was shown a property on the commutators which characterizes the classical implies that the order conditions for the proposed splitting are the conditions given in Table 1 are sufficient to assure the fifth order
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
837
Table 2 Additional trees, elementary weights and density function values up to order five for method (10) (top) and for method (9) (bottom).
τ
Φ(τ ) s+1 i −1
γ (τ ) (b)
(a)
b i a j c j (c j + α j )
12
i =1 j =1 j s+1 i −1
(a)
(a)
b i a j bl cl (c j + α j )
40
i =1 j =1 l=1 s+1 i −1
(b)
(a)
(a)
b i a j c j ((c j )2 + 2α j c j + α j a j )
20
i =1 j =1
τ
Φ(τ ) s+1 i −1
γ (τ ) (b)
b i a j (c j )2
12
i =1 j =1 j s+1 i −1
(b)
(a)
b i a j c j bl c l
40
i =1 j =1 l=1 s+1 i −1
(b)
b i a j (c j )3
20
i =1 j =1
accuracy of the method (17). However, the method (17) cannot be viewed as the application of an ESPRK method, since it includes the integration of the explicit time dependence. Nevertheless, by replacing the integral according to a quadrature rule in [0, 1] with weights w k and knots sk for k = 1, . . . , r,
qn+1 = qn + h
r s
+ ai sk P i
(a)
ai w k M tn + h c i
i =1 k =1
it is possible to obtain methods of the desired order (see [5]) which, however, do not correspond to any application of an ESPRK method defined by (7) for solving (2). In the followings we propose an alternative formulation of the problem (1) that allows to reproduce the optimal performance of ESPRKN methods. Our approach is based on the main observation that, from the relations
1
(a)
M tn + h c i
hai
(a)
tn +hc i +1
+ ai s ds =
0
M (s) ds, (a)
tn +hc i
(a) tn +hc i +1
(b)
M (s) ds = (a)
M (s) ds + (a)
tn +hc i
(a)
tn +hc i +1
tn +hc i
tn +hc i
M (s) ds,
(18)
(b)
tn +hc i
it results that the method (9) performs the first order rectangular right and left quadrature approximations of the integrals in the right-hand side of (18). Similarly, the method (10), corresponds to the first order rectangular left and right quadrature approximations. Let us generalize this r approach by considering a r-step quadrature rule with weights w r and knots sr . Assuming a symmetric rule, with i =1 w r = 1 and sr −i +1 = 1 − s i for i = 1, . . . , r, the problem (6) can be formulated as the equivalent partitioned system
⎧ q˙ = 1, ⎪ ⎨ 0 r q˙ = w j M ( s j q 0 + s r − j +1 p 0 ) p , ⎪ ⎩ j =1
p˙ 0 = 1, p˙ = f (q0 , q).
(19)
We claim that Proposition 4.1. A 2r order ESPRKN method preserves its accuracy when applied to non-autonomous problems (2) written in the equivalent form of partitioned systems (19) where w j and s j are the weights and knots, respectively, of a 2r order symmetric quadrature rule in [0, 1].
838
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
(b)
Proof. Firstly, consider that the amplitudes of the integrals in (18) are h(c i since
(a)
(a)
(b)
− c i ) = hβi , and h(c i +1 − c i ) = hαi . Moreover,
(a)
(b) (a) + sr − j +1 tn + hc i = tn + hc i + hβi sr − j +1 , (a) (b) = tn + hc i(b) + hαi s j , s j tn + hc i +1 + sr − j +1 tn + hc i s j tn + hc i
it results that the application of the 2r order ESPRKN methods to the partitioned system (19) leads to a method which evolves according to the following approximations
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i Q i = Q i −1 + β i h
r
,
(a)
w j M tn + hc i
+ h β i s r − j +1 P i
j =1 r
+ αi h
(b)
w j M tn + hc i
+ h αi s j P i .
(20)
j =1
Since w j = w r − j +1 for j = 1, . . . , r, in the Q i ’s formula we recognize 2r order quadrature rules applied to the integrals in (18). Hence, the method (20), whose coefficients are provided by a 2rth order ESPRKN method, is a 2rth order approximation of the method (17), therefore preserves its accuracy. 2 5. Numerical tests The numerical examples herein are aimed at showing the preservation of the accuracy of ESPRKN methods when applied as in (20) for solving non-autonomous second order dynamics (2). To fully exploit the results given in Proposition 4.1, we consider methods defined by symmetric coefficients. That way the 2rth order methods are obtained by satisfying the conditions for the 2(r − 1)th order. In particular, we consider the optimal 11-stage symmetric fourth order ESPRKN method b S R K N 11 given in [7] which exhibits a sixth order behaviour on autonomous problems (3). By means of its coefficients (6)
we define as S R K Q N 11 the method that exploits in (20) the approximation provided by the sixth order Gauss–Legendre quadrature rule:
Q i = Q i −1 + h β i
+ h αi +
5−
5 18
√
5 18 15
√ 5 − 15 (a) (b) βi M tn + hc i + hβi + αi M tn + hc i
10
(a)
h βi
h αi
(b)
βi M tn + hc i + + αi M tn + hc i + 2 2 √ √ 5 + 15 5 + 15 (a) (b) Pi. βi M tn + hc i + hβi + αi M tn + hc i + hαi 10
+
4
9
10
10
b We will compare the results of the proposed method with the performances of S R K N 11 when applied according to the classical approach (9). We measure the 2-norm error on the numerical approximations with respect to the solution computed by the Matlab built-in function ODE45 with tolerances set at the machine precision. We will show that the order reduction we can appreciate on the classical implementation, when the effects of the time dependent matrix M (t ) are increased, can (6) be nullified by the proposed method S R K Q N 11 .
5.1. Example 1: Time dependent damped harmonic oscillator: Autonomous damping force linear in the velocity As a test problem we choose the Duffing oscillator (see [12]):
q˙ = M (t ) p ,
p˙ = − M −1 (t ) q3 − q − δ cos(ωt ) .
(21)
with M (t ) = e − t and 0 1. By varying the values of , the dependence of M (t ) on time is amplified or destroyed. Indeed, for values of ≈ 0 the problem behaves as the second order differential equation
q˙ = p ,
p˙ = q3 − q − δ cos(ωt )
that can be accurately integrated with ESPRKN methods. For increasing values of presence of the exponential term. Notice that, from the relation
d dt
˙ M −1 q˙ + M −1 q¨ M −1 q˙ = − M −1 M
their behavior is more affected by the
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
839
˙ M −1 , the problem (1) can be equivalently formulated as and by setting A = M q¨ = A (t )˙q + M (t ) f (q, t )
(22) e − t , hence it results
that is the classical formulation of a damped oscillator. In the case of the Duffing oscillator, M (t ) = that A (t ) = − . In this case, the amplification of the time dependence of M (t ) in (21) corresponds to the amplification of the ‘damped’ structure of the oscillator (22). Again, for increasing values of , ESPRKN methods feel the effects of the autonomous damping term − q˙ and they can lose the order they show when used to integrate the undamped oscillator
q¨ = q3 − q − δ cos(ωt ) that corresponds to = 0. In our test we perform the integration in the interval [0, 10π ] with initial conditions q(0) = 1.75 and p (0) = 0. In order to verify the order of the different methods, the initial step h = π4 is successively halved and the maximum error is plotted against the decreasing stepsizes. In Fig. 1(a), (b), (c) the values of the parameter δ is set at 1/4 and ω = 1, while different choices for parameter are considered. Firstly, we suppose = 1e − 4 in order to neglect the explicit time dependence of M (t ); in this case both the considered methods show a sixth order accuracy as expected and as it can be seen in Fig. 1(a). The cases corresponding to = 0.1, 0.5, are illustrated in Fig. 1(b), (c): we show both the fourth order accuracy of the approximation provided by (6) b and the theoretically predicted sixth order accuracy of the new method S R K Q N 11 . The difference in the method S R K N 11 performance of the two methods is enhanced for increasing values of the parameter. 5.2. Example 2: Time dependent damped harmonic oscillator: Non-autonomous damping force linear in the velocity We consider the one-dimensional time dependent harmonic oscillator given in [20]. The model is defined by M (t ) = t e − L (t ) where L (t ) = 0 sin(τ /π )dτ and 0 < < 1; moreover f (q, t ) = − M (t )−1 ω2 (t )q where ω(t ) = cos(t /2). The equivalent formulation (22) corresponds to
q¨ = − sin(t /π )˙q − ω2 (t )q; Also in this case, the effect of the non-autonomous damping term − sin(t /π )˙q is amplified for increasing values of . The integration is performed in the interval [0, 100π ], with initial step h = π4 and initial conditions q(0) = 1, p (0) = 0. Fig. 1(d), (e), (f) corresponds to different values of parameter : at = 0.018π that defines the model in [20], both the methods show a sixth order accuracy; when the parameter is set at = 0.1π the approximations provided by the method (6) (9) feel the effects of this increment and begin to lose their accuracy. At = 0.5π the S R K Q N 11 outperforms the method b S R K N 11 providing sixth order approximations.
5.3. Quantum damped oscillators: The Walker–Preston model The Walker–Preston model of a diatomic molecule in a strong laser field [21] corresponds to the one-dimensional Schrödinger equation
i
∂ 1 ∂2 ψ(x, t ) = − + V ( x ) + f ( t ) x ψ(x, t ) ∂t 2μ ∂ x2
where h¯ = 1 and ψ(x, 0) = ψ0 (x). In the model, the Morse potential V (x) = D (1 − e −α x )2 and f (t ) = A cos(γ t ) are considered. The parameters, in atomic units a.u, are set at the values reported in [11,19,4] i.e. μ = 1745, D = 0.2251, α = 1.1741, A = 0.011025 and γ = 0.01787. The spatial interval [−0.8, 4.32] is divided in N = 64 parts of length x = 0.08 and periodic boundary conditions are imposed. The initial condition is given by
ψ0 (x) = σ e −(Ω−0.5)α x e −Ω e
−α x
√
with Ω = 2D /ω0 , ω0 = α 2D /μ and σ a normalizing constant. After a space discretization by means of finite differences, we obtain the following system of ODEs:
i u˙ = H (t )u where H is a real time dependent symmetric matrix and the vector u (t ) ∈ C N approximates ψ(x, t ) at the nodal values xn . The vector components of the initial condition u (0) are given by un (0) = ψ0 (xn ) for n = 1, . . . N. In order to formulate the model as a problem described by (2), we use two set of N variables q, p which correspond to the real and the imaginary parts of u:
q˙ = H (t ) p , p˙ = − H (t )q.
840
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
Fig. 1. Left: Duffing oscillator. Efficiency comparison at
= 0.018π (d), = 0.1π (e), = 0.5π (f).
= 1e − 4 (a), = 0.1 (b), = 0.5 (b). Right: damped harmonic oscillator. Order comparison at
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
Fig. 2. Quantum damped oscillator,order comparison. Left:
841
γ = 0.01787. Right: γ = 0.2.
We integrate the system in the time interval [0, T ] with T = 200π and initial stepsize h = 2.5π . We measure the 2-norm error u teor ( T ) − u ( T ) , where the values of u teor are computed, using the same spatial discretization, by the function ODE45. In this as well as in the following case, since the most demanding part of the problem is the function evaluation we plot the efficiency curves (error vs. number of function evaluations). For numerical comparisons, in addition to the (6) (2) (4) proposed method S R K Q N 11 we also consider the methods S R K Q N 11 and S R K Q N 11 corresponding to the second and
b method. In the fourth order Gauss–Legendre approximations. We compare the results with the performances of the S R K N 11
(6)
Fig. 2(left), we can appreciate the sixth order accuracy of the proposed method S R K Q N 11 for the chosen set of parameters. (4)
Notice that the excellent performance of the S R K Q N 11 , which actually should be a fourth order method, last unabated for a wide range of the parameters. In order to discriminate between a fourth and a sixth order effective behaviour we have to force the frequency γ , which weigths the explicit time dependence, to a quite higher γ = 0.2. Results are shown in Fig. 2(right). 5.4. Quantum damped oscillators: The modified Caldirola–Kanai oscillator 1 − sin γ t The modified Caldirola–Kanai oscillator (see [15]) can be considered as a variable mass m(t ) = μ e system and the corresponding the Schrödinger equation (h¯ = 1) is given by
i
∂ m(t ) ∂ 2 ψ(x, t ) = − + V ( x , t ) ψ(x, t ). ∂t 2 ∂ x2
The potential V (x, t ) = 12 μω(t )2 e sin γ t x2 , is defined by means of ω(t ) = A cos(γ t ) with A = 0.001, γ = 0.01, μ = 1745. The second derivative is approximated with the central difference on the same spatial meshgrid considered in the previous example and we integrate the system for the real and imaginary parts of the vector u (t ), starting from the same initial conditions. The time integration is in the interval [0, T ] with T = 10π and large initial stepsize h = 4π . As for the previous example, the proposed method has been implemented by means of quadrature formulae of second, fourth and sixth (6) order, respectively. Fig. 3(left), confirms the sixth order accuracy of the proposed method S R K Q N 11 and the good perfor(4)
mance of the S R K Q N 11 method for the given set of parameters. In Fig. 3(right), with a different choice of the parameters ( A = 0.011025, γ = 0.2 and initial stepsize h = π ), all the methods behaves as fourth order ones, except for the method (6)
S R K Q N 11 which still preserves the sixth order, as expected. 6. Conclusions We have considered explicit symplectic partitioned Runge–Kutta methods for the numerical integration of a specific class of non-autonomous second order differential equations. As concerns splitting and composition methods applied to these problems, in [6] it has been shown that the way the explicit time dependence is treated can be critical. Herein, by performing the P -series analysis, we theoretically demonstrate that the ‘natural’ choices which come straightforward from integrating the time with one or the other set of variables p and q, introduce additional elementary differentials that increase the number of conditions to be imposed at a given order. Hence, by treating the time as two added system’s coordinates, we propose a simple strategy which allows to use ESPRK tailored on dynamics q˙ = p ; p˙ = f (q, t ), as integrators
842
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
Fig. 3. Caldirola–Kanai oscillator, order comparison. Left: A = 0.001,
γ = 0.01, initial stepsize h = 4π . Right: A = 0.011025, γ = 0.2, initial stepsize h = π . (6)
for the class of problems q˙ = M (t ) p ; p˙ = f (q, t ). The performances of the proposed sixth order method S R K Q N 11 have
b been compared with the classical integrator S R K N 11 method on the solution of several test problems.
Acknowledgements We would like to thank the referees for the positive suggestions and comments we believe helped us to greatly improve the presentation of the work as well as to clarify the main results. Appendix A In this appendix we provide the proof of Proposition 3.1 given in Section 3. Proof of Proposition 3.1. Denote with z w the tree ( z ◦ u ) ◦ v u and notice that z w = [ω1 , ω2 , . . . , ωn , u , v u ] w . Then, the elementary weight Φ related to method (7) is provided by
Φ( z w ) =
s
βi Φ2i −1 ( z w ) +
s
i =1
αi Φ2i (z w )
i =1
or, (α0 = βs+1 = 0) equivalently,
Φ( z w ) =
s +1
βi Φ2i −1 ( z w ) +
s +1
i =1
αi−1 Φ2i−2 (z w ),
(23)
i =1
wherewe set Φ0 ( z w ) = 0, Φk ( z w ) = Πk Ψk (u )Ψk ( v u ) for k = 1, . . . , 2s and Πk = 1 when z is an elementary tree, otherwise Πk = ln=1 Ψk (ωl ). From the definition of Ψk (see [13]), it follows that, for all τ ∈ T P b ,
Ψ1 (τ ) = 0, Ψ2i (τ ) =
i
βm Φ2m−1 (τ ) +
m =1
Ψ2i +1 (τ ) = Ψ2i (τ ),
i
αm Φ2m (τ ), i = 1, . . . , s,
m =1
i = 1, . . . , s − 1.
Since u , v u , ωl ∈ T P b , for l = 1, . . . , m, form (24) it follows that Φ( z w ) in (23), is given by
Φ( z w ) =
s +1 (βi + αi −1 )Π2i −1 Ψ2i −1 (u )Ψ2i −1 ( v u ) i =1
=
s +1
b i Π2i −1 Ψ2i −1 (u )Ψ2i −1 ( v u )
i =1
=
s i =1
b i +1 Π2i Ψ2i (u )Ψ2i ( v u ).
(24)
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
843
As in [9], we notice that (see [18]) stopping the method (7) at the 2ith stage with 1 i s, a new ESPRK (inconsistent) method of the same form is obtained. From (24), it results that the elementary weight of the new method, say Φ (2i ) (τ ), is the elementary weight Ψ2i (τ ) associated with the 2i stage of (7). As a consequence, we reformulate Φ( z w ) as
Φ( z w ) =
s +1
b i +1 Π2i Φ (2i ) (u )Φ (2i ) ( v u ).
i =1
Since relation (16) still holds when replacing Φ(τ ) with Φ (2i ) (τ ), relation (16) follows.
2
In a similar way, we are able to demonstrate Proposition 3.2. Proof of Proposition 3.2. By adopting the previous notations, we have that
Φ( z ◦ r ) =
s i =1
b i +1 Π2i Ψ2i (r ) =
s
b i +1 Π2i Φ (2i ) (r ).
i =1
Since Φ(u ) + Φ( v ) = Φ(r ) still holds when replacing Φ(τ ) with Φ (2i ) (τ ), the result follows.
2
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
L. Abia, J.M. Sanz-Serna, Partitioned Runge–Kutta methods for separable Hamiltonian problems, Math. Comput. 60 (1993) 617–634. J. Awrejcewicz, D. Sendkowski, Stick-slip chaos detection in coupled oscillators with friction, Int. J. Sol. Struct. 42 (2005) 5669–5682. S. Blanes, F. Casas, Splitting methods for non-autonomous separable dynamical systems, J. Phys. A: Math. Gen. 39 (2006) 5405–5423. S. Blanes, F. Casas, P.C. Moan, Splitting methods for non-autonomous linear systems, Int. J. Comput. Math. 84 (2007) 713–727. S. Blanes, F. Casas, A. Murua, Splitting methods in the numerical integration of non-autonomous dynamical systems, RACSAM: Rev. R. Acad. Cien. Serie A. Mat. (2010), submitted for publication. S. Blanes, F. Diele, C. Marangi, S. Ragni, Splitting and composition methods for explicit time dependence in separable dynamical systems, J. Comput. Appl. Math. 235 (2010) 646–659. S. Blanes, P.C. Moan, Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyström methods, Appl. Math. 142 (2002) 313–330. J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations. Runge–Kutta and General Linear Methods, Wiley, Chicester, 1982. M.P. Calvo, E. Hairer, Further reduction in the number of independent order conditions for sympletic, explicit partitioned Runge–Kutta and Runge– Kutta–Nyström methods, App. Num. Math. 18 (1995) 107–114. J.M. Cervero, J. Villarroel, On the quantum theory of the damped harmonic oscillator, J. Phys. A: Math. Gen. 17 (1984) 2963–2971. S. Gray, J.M. Verosky, Classical Hamiltonian structures in wave packet dynamics, J. Chem. Phys. 100 (1994) 5011–5022. J. Guckenheimer, P. Holmes, Nonlinear Oscillators, Dynamical Systems, and Bifurcations of Vector Fields, Springer-Verlag, New York, 1983. E. Hairer, Backward analysis of numerical integrators and symplectic methods, Annals of Numerical Mathematics 1 (1994) 107–132. E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations, Springer-Verlag, Berlin, 2002. A.N. Ikot, I.O. Akpabio, M.I. Umo, E.E. Ituen, Quantum damped mechanical oscillator, Int. J. Optics (2010), doi:10.1155/2010/275910. P.G.L. Leach, Harmonic oscillator with variable mass, J. Phys. A: Math. Gen. 16 (1983) 3261–3269. A. Murua, On order conditions for partitioned symplectic methods, SIAM J. Numer. Anal. 34 (1997) 2204–2211. D. Okunbor, R.D. Skeel, Explicit canonical methods for Hamiltonian systems, Math. Comput. 59 (1992) 439–455. J.M. Sanz-Serna, A. Portillo, Classical numerical integrators for wave-packet dynamics, J. Chem. Phys. 104 (1996) 2349–2355. J. Struckmeier, C. Riedel, Canonical transformations and exact invariants for time-dependent Hamiltonian systems, Ann. Phys. 11 (2002) 15–38. R.B. Walker, K. Preston, Quantum versus classical dynamics in treatment of multiple photon excitation of anharmonic-oscillator, J. Chem. Phys. 67 (1977) 2017–2028.
Contents lists available at ScienceDirect
Applied Numerical Mathematics www.elsevier.com/locate/apnum
Explicit symplectic partitioned Runge–Kutta–Nyström methods for non-autonomous dynamics ✩ Fasma Diele, Carmela Marangi ∗ Istituto per le Applicazioni del Calcolo M. Picone, CNR, Via Amendola 122, 70126 Bari, Italy
a r t i c l e
i n f o
a b s t r a c t
Article history: Received 23 June 2010 Received in revised form 26 January 2011 Accepted 5 February 2011 Available online 18 February 2011 Keywords: Symplectic partitioned Runge–Kutta methods Order analysis Nyström methods
We consider explicit symplectic partitioned Runge–Kutta (ESPRK) methods for the numerical integration of non-autonomous dynamical systems. It is known that, in general, the accuracy of a numerical method can diminish considerably whenever an explicit time dependence enters the differential equations and the order reduction can depend on the way the time is treated. In the present paper, we demonstrate that explicit symplectic partitioned Runge–Kutta–Nyström (ESPRKN) methods specifically designed for d ( M −1 q˙ ) = f (q, t ), undergo an order reduction when second order differential equations dt M = M (t ), independently of the way the time is approximated. Furthermore, by means of symmetric quadrature formulae of appropriate order, we propose a different but still equivalent formulation of the original non-autonomous problem that treats the time as two added coordinates of an enlarged differential system. In so doing, the order reduction is avoided as confirmed by the presented numerical tests. © 2011 IMACS. Published by Elsevier B.V. All rights reserved.
1. Introduction Second order differential equations and their numerical solution play a crucial role in a variety of applications in science and engineering. In this paper we focus on the numerical solution of initial value problems for the special second order equations
d dt
M −1 (t )˙q = f (q, t ),
t ∈ [0, T ]
(1)
where the evolution of the variable q = q(t ) ∈ Rd depends on time by means of the time dependent, non-singular, continuously differentiable matrix M : [0, T ] → Rd×d and through the function f : Rd × [0, T ] → Rd sufficiently smooth in order to assure the well-posedness of the problem. By setting p = M −1 (t )˙q the problem can be recast into the following first order differential system
q˙ = M (t ) p ,
p˙ = f (q, t ).
(2)
Both Eqs. (1) and (2) allow to analyze dynamics relevant to nonlinear science whenever time dependent damped or driven harmonic oscillators are considered. Several examples can be found in electronics, in optics and in classical mechanics literature. Single and coupled oscillators with friction (see [10,2,12]) as well as variable mass dynamics (see [16]), can be described by means of the equations above. Moreover, in quantum mechanics they can describe, after a spatial discretization, the evolution of a time dependent Schrödinger equation related to quantum damped oscillators (see [15,21]). ✩
*
This work has been supported in part by National Research Council grant CNR-RSTL id.332/2008. Corresponding author. E-mail addresses: [email protected] (F. Diele), [email protected] (C. Marangi).
0168-9274/$30.00 © 2011 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2011.02.003
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
833
In case of a time independent matrix M (t ) = M 0 , it is easy to see that the original system (1) can be written as an autonomous second order equation Q¨ = F ( Q ) by treating the time in f (q, t ) as an added coordinate q0 = t, setting Q = (q0 , q), and denoting F = (0, M 0 f ( Q )). The corresponding first order differential system, obtained by setting p = q˙ and P = (1, p ) can be written as
Q˙ = P ,
P˙ = F ( Q ).
(3)
A large number of numerical methods for the integration of such equations can be found in literature. They exploit the two main features of such a kind of problems, that are the separability of the variables Q and P , as well as the linear dependence on P . In particular, powerful methods were provided in [7], where an optimal efficiency was obtained through the minimization of the effective error. Unfortunately, in case of a time dependent matrix M (t ), the original non-autonomous problem (1) cannot be reformulated so to recover the useful algebraic structure of (3). In this case, the classical way to face the problem is to further enlarge the system by introducing another time variable p 0 to account for the explicit time dependence in M (t ). By this trick we can gain the separability but we lose the linear dependence on one of the two sets of variables describing the enlarged system. As a consequence, in the Taylor series expansion of the exact solution, some extra elementary differentials appear which determine further order conditions that are not satisfied by an ESPRKN method. Different numerical approaches have been developed in literature in order to design effective numerical methods for non-autonomous dynamics (2) (see [3,4]). Notice that all the methods referred above are symplectic methods. It is known that, for a given problem of a given algebraic structure, symplectic methods turns out to be the most accurate ones with respect to any other non-symplectic PRK method (see e.g. [14]). As a matter of fact, the problems we are dealing with in the present paper, can all be described in an hamiltonian framework by introducing the function
H ( p , q, t ) =
1 2
p T M (t ) p + V (q, t )
where M (t ) is taken as a positive definite symmetric matrix and −∇q V (q, t ) = f (q, t ). In the special case of autonomous systems (3), the Hamiltonian function is conserved and symplectic methods outperform any other method in preserving the qualitative features of the theoretical flow. Whenever non-autonomous systems (2) are considered the corresponding Hamiltonian system loses some of its features as, for example, the conservation of the energy. In this context it is reasonable to wonder about the real advantage of using symplectic methods. The answer lies in the P -series order analysis of PRK methods where the graphical representation of elementary differentials as bicoloured trees is adopted. Indeed, by means of symplecticness one is able to prove useful relations among trees’ elementary weights useful to reduce the number of conditions at a given order [17]. So, even in the context of non-autonomous systems, symplecticness plays a major role in the reduction of the order conditions. In the recent work [6], developed in the framework of splitting and composition techniques, the authors proposed an optimal strategy to deal with the explicit time dependence. Such a strategy, based on the introduction of a single time variable q0 , loses the separability, but regains the linear dependence. Then, starting from those results, one may wonder if there is any special advantage in the linearity on P with respect to separability, that can be exploited to build optimal PRK methods for this class of problems. In order to provide an answer, in this paper we compare the classical approach, where two time coordinates, q0 and p 0 are added to the system, to an alternative one that makes use of a single time variable q0 as in [6]. As a matter of fact, thanks to the symplecticness, in our paper we are able to provide new results that allow us to establish the number of order conditions for the non-separable problem considered. Unfortunately, we find out that the advantages obtained form the linear dependence on P are the same one gets from the separability of the two sets of variables, i.e. the same number of order conditions is obtained in both cases. As a main result of the paper we propose a new strategy to approximate the time, based on the introduction of two time coordinates, as in the classical approach, and the use of symmetric quadrature formulae. We prove that an ESPRKN method applied on a new (still equivalent) formulation of the problem (3), preserves its accuracy avoiding order reduction. Notice that the proposed strategy for the integral approximation is based on symmetric quadrature formulae in order to preserve the symmetry of the resulting ESPRK method. We recall that, for symmetric methods, a 2nth order method can be obtained by satisfying the order conditions for the (2n − 1)th order. In details, in Section 2, the integration of the explicit time dependence in the framework of explicit symplectic partitioned Runge–Kutta methods (ESPRK) is addressed; in Section 3 we perform the P -series analysis of an ESPRK method applied to the non-separable problem obtained by introducing the only variable q0 and we establish the order conditions of a method up to the fifth order. In Section 4 we suggest how to modify the original problem in order to gain the desired order whenever an ESPRKN method is used to approximate the solutions of the non-autonomous system (2). Theoretical results are validated in Section 5 on the solution of perturbed time dependent oscillators as well as on time dependent linear Schrödinger equations chosen as test problems. Finally, in Section 6 we draw our conclusions.
834
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
2. Classical implementation of ESPRK methods for general non-autonomous dynamics ( A)
We recall that a partitioned symplectic k-stages Runge–Kutta method is defined by two RK methods ( A i j , b i
( B i j , b(i B ) , c (i B ) ), which satisfy the symplecticity conditions (see e.g. [14]): ( A)
bi
(B)
(B)
( A) (B )
B i j + b j A i j = bi
bj ,
( A)
, c i ) and
i , j = 1, . . . , k ,
(4)
( A)
i = 1, . . . , k , (5) ( A) (B) where c i = = j B i j are the nodes of the two Runge–Kutta methods. An important feature of PRK methods j A i j and c i bi
= bi ,
is the possibility of having explicit methods for non-autonomous dynamics
q˙ = g ( p , t ),
p˙ = f (q, t ).
(6)
I.e. dynamics which reduce to separable ones whenever there is no explicit time dependence in one or both the equations describing the q’s and p’s evolution. ( A) (B) 0. The These methods are based on two diagonally implicit methods satisfying A ii · B ii = 0 with b i = 0 and b i = classical form of a k = 2s stage method is given in literature by the following Butcher arrays
0
β1 β1 β1 β1 .. . β1 β1
α1 α1 0 α1 β2 α2 α1 β2 α2 0 .. .
.. . α1 β2 α1 β2
.. .
.. . α2 · · · α2 · · ·
..
. βs αs βs αs
β1 β1 β1 β1 .. . β1 β1 β1
0
α1 β2 α1 β2 0 .. .
.. . α1 β2 α1 β2 α1 β2
.. . ··· ··· ···
..
(7)
.
α s −1 β s α s −1 β s 0 α s −1 β s α s
which define the A (left) and B (right) RK methods with αi = 0 and βi = 0 for i = 1, . . . , s. It is easy to verify that the partitioned method (7) satisfies both symplecticity conditions (4)–(5). For non-autonomous dynamics (6), the application of an explicit symplectic Runge–Kutta method, defined by (7), to the partitioned system
q˙ 0 = 1, q˙ = g ( p 0 , p ),
p˙ 0 = 1, p˙ = f (q0 , q)
where (q0 , q T ) ∈ Rd+1 are solved by the A method and ( p 0 , p T ) ∈ Rd+1 variables are solved by the B method, leads to the following steps
P 0 = pn ,
Q 0 = qn ,
Do i = 1, s + 1
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i
(b)
Q i = Q i −1 + hai g P i , tn + hc i enddo p n +1 = P s +1 , where, by setting (b)
and c i
=
i
qn+1 = Q s+1 ,
t n +1 = t n + h , (a)
α0 = βs+1 = 0, we define bi = αi−1 + βi , for i = 1, . . . , s + 1, ai = αi + βi , for i = 1, . . . , s, c i =
k=1 b i .
Moreover, we set as+1 = 0. (a)
Notice that the nodes of the Runge–Kutta formula (7) are different from the coefficients c i (a)
( A)
(a)
( A) = c 2i −2
( A) = c 2i −1 for i
(b)
and c i
(8) i −1 k=1 ai
for i = 1, . . . , s + 1 (b)
(B)
(B)
which appear in (9). They are related as follows: c 1 = c 1 = 0, c i = 2, . . . , s + 1 and c i = c 2i −1 = c 2i for i = 1, . . . , s + 1. This is the classical implementation of non-autonomous problems: the time acts as two additional coordinates treated as p and q respectively. This way the original non-separable problem can be treated as a separable one in an enlarged space. However, the classical partitioning is simply one of the possible forms in which an ESPRK method can be rewritten in case of non-autonomous dynamics. One may as well treat the time as a single additional coordinated evolving with q in both equations of (6) as in [6], where this particular coupling of the time variable is what provides the optimal performance for the splitting and composition methods when applied to non-autonomous problems (1). Since the introduction of a single time variable leads to a non-separable problem, it is not so trivial to infer the performance of an ESPRK method for the case at hand. We will see that, despite what can be proved in the case of splitting and composition methods, this approach does
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
835
not provide advantages with respect to the classical implementation. On the other hand, we are able to demonstrate that the number of order conditions does not increase with respect to the classical partitioning, even if we are dealing with a non-separable problem. 3. A different implementation and relative order analysis Eqs. (2) belongs to the class of non-autonomous dynamics of type (6) where g ( p , t ) = M (t ) p. We search for approximate solution obtained by applying an ESPRK method defined by (7). As we noticed before, different implementations are possible according to how the time variable is approximated. The classical implementation, described in the previous section is provided by the following inner stages:
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i
(b)
Q i = Q i −1 + hai M tn + hc i
,
Pi.
(9)
The implementation we focus on, instead, corresponds to
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i
(a)
Q i = Q i −1 + h βi M tn + hc i
a) + αi M tn + hc i(+ Pi. 1
(10)
The above method describes the application of an ESPRK to the partitioned system
q˙ 0 = 1, q˙ = M (q0 ) p ,
p˙ = f (q0 , q)
(11)
where the time variable is integrated together with the q coordinate in both equations of (6). To perform the order analysis of the partitioned Runge–Kutta (PRK) method (10) applied to the differential system (11), we follow the classical procedure as given in [14]. Notice that, the application of the theory in this context is not straightforward. In order to get the minimal set of order conditions, we will extend the results given in [9], concerning separable systems, to the non-separable problem at hand and provide new relations that allow to properly reduce the number of necessary order conditions. The order of a PRK method can be characterized in terms of its coefficients via bi-coloured rooted trees. We represent with a ‘white’ vertex the elementary tree representing a f = f (q0 , q), and with a ‘black’ vertex the elementary tree representing a g = M (q0 ) p. The rules for constructing trees, for the case at hand, are the following: - ◦ vertices can have only • children. - • vertices can have both ◦ and • children; ◦ children cannot be more than one, • children are elementary trees. As can be seen from the second rule, the introduction of time affects the separability of the problem and is responsible for the peculiar constraint on black children of black vertices. So, unless we are able to prove that such conditions are enough to drastically reduce the number of trees, we should expect a number of order conditions similar to the one that features an ESPRK method when applied to a generic non-separable problem. In order to apply the order theory within the new set of rules, we denote by T P b and T P w all the bi-coloured trees with black and white roots, respectively. For symplectic PRK methods, from condition (4) and for all u ∈ T P b and v ∈ T P w , it follows that [1]
Φ(u ◦ v ) + Φ( v ◦ u ) = Φ(u )Φ( v )
(12)
where Φ(τ ) denotes the elementary weight associated to the tree τ and ◦ corresponds to the Butcher’s product of trees (see [8]). Moreover, because of the symplecticity property (5), Φ(τ ) is independent of the type of the root of τ . More in details, define u = [τ1 , τ2 , . . . , τn ]b a generic tree in T P b obtained by connecting the roots of τ1 , τ2 , . . . , τn to a new black root and, similarly, v = [σ1 , σ2 , . . . , σl ] w a generic tree in T P w obtained by connecting the roots of σ1 , σ2 , . . . , σl to a new white root. Moreover, suppose that v u = [σ1 , σ2 , . . . , σl ]b ∈ T P b , and u v = [τ1 , τ2 , . . . , τn ] w ∈ T P w , where v u and u v are obtained from v and u by changing the colour of their root. Then, the following relation holds
Φ( v ) = Φ( v u )Φ(u ) = Φ(u v ).
(13)
A direct consequence of (12) and (13) is the relation
Φ(u ◦ v ) + Φ( v u ◦ u ) = Φ(u )Φ( v u )
(14)
which leads to the following further simplifying assumptions which generalize the results given in [9]. Proposition 3.1. Suppose u , ω1 , ω2 , . . . , ωn , ∈ T P b , v ∈ T P w . Let z be the elementary white tree or z = [ω1 , ω2 , . . . , ωn ] w ∈ T P w . Then, the following relation holds
836
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
Table 1 Nyström trees, elementary weights and density function values up to order five.
τ ◦
Φ(τ ) s+1
bi =
i =1
s
γ (τ )
ai
τ
s+1
1
i =1
s+1
Φ(τ )
γ (τ ) (a)
bi c i
2
i =1
(a) 2
bi c i
i −1 s+1
3
i =1
(b)
bi a j c j
6
i =1 j =1
s+1
(a) 3
bi c i
s+1 i −1 i
4
i =1
(a)
b i a j bk c k
24
i =1 j =1 k=1
s+1
(a) 4
bi c i
j l−1 s+1 i −1
5
i =1
(b)
b i a j bl ak ck
120
i =1 j =1 l=1 k=1
j s+1 i −1
(a)
(a)
b i c i a j bl c l
j s+1 i −1
30
i =1 j =1 l=1
(a) 2
b i a j bl c l
60
i =1 j =1 l=1
Φ z ◦ (u ◦ v ) + Φ z ◦ ( v u ◦ u ) = Φ ( z ◦ u ) ◦ v u
(15)
whenever v u and ( v u ◦ u ) ∈ T P b . Proposition 3.2. Suppose u , v , r ∈ T P b that verify Φ(u ) + Φ( v ) = Φ(r ). Let z be the elementary white tree or z = [ω1 , ω2 , . . . ,
ωn ] w ∈ T P w . Then, the following relation holds Φ( z ◦ u ) + Φ( z ◦ v ) = Φ( z ◦ r ).
(16)
For sake of clarity, we refer the reader to Appendix A for the proofs. By applying these results to our problem we find out that, with respect to the trees which assure the fifth order accuracy to a ESPRKN method (reported in Table 1 together with the elementary weights as well as their density functions), new trees appears and are given in Table 2 (top). In the same table (bottom), the (well-known) additional conditions related to the classical separable method (9) are reported for comparison. Both approaches share the same total number of order conditions i.e. 14, even if they are related to different trees, with different weights. Notice that here we are dealing with a non-separable method (10). Remembering that in the case of a generic non-separable problem we may have up to 40 order conditions [17] for symplectic methods, we proved that the method here analyzed does not perform worse than the classical separable one (9). On the other hand, we observe that the advantages one can get within the splitting and composition framework, cannot be recovered by applying an ESPRKN method as in (10). 4. The method of exact integration and its approximations When the coefficients ai and b i of an ESPRKN method are used for the evaluation of
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i
1 Q i = Q i −1 + hai
(a)
M tn + hc i
,
+ hai s ds P i ,
0
the resulting approximations
p n +1 = p n + h
s +1
(a)
b i f Q i −1 , tn + hc i
,
i =1
qn+1 = qn + h
s i =1
1
(a)
M tn + h c i
ai
+ ai s ds P i ,
(17)
0
inherit the accuracy of the ESPRKN method [6]. Indeed, that the splitting underlying the method (17), satisfies splitting of autonomous second order equations (3). This same of classical ESPRKN methods. In particular, the 11
by performing the order analysis via Lie formalism, it was shown a property on the commutators which characterizes the classical implies that the order conditions for the proposed splitting are the conditions given in Table 1 are sufficient to assure the fifth order
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
837
Table 2 Additional trees, elementary weights and density function values up to order five for method (10) (top) and for method (9) (bottom).
τ
Φ(τ ) s+1 i −1
γ (τ ) (b)
(a)
b i a j c j (c j + α j )
12
i =1 j =1 j s+1 i −1
(a)
(a)
b i a j bl cl (c j + α j )
40
i =1 j =1 l=1 s+1 i −1
(b)
(a)
(a)
b i a j c j ((c j )2 + 2α j c j + α j a j )
20
i =1 j =1
τ
Φ(τ ) s+1 i −1
γ (τ ) (b)
b i a j (c j )2
12
i =1 j =1 j s+1 i −1
(b)
(a)
b i a j c j bl c l
40
i =1 j =1 l=1 s+1 i −1
(b)
b i a j (c j )3
20
i =1 j =1
accuracy of the method (17). However, the method (17) cannot be viewed as the application of an ESPRK method, since it includes the integration of the explicit time dependence. Nevertheless, by replacing the integral according to a quadrature rule in [0, 1] with weights w k and knots sk for k = 1, . . . , r,
qn+1 = qn + h
r s
+ ai sk P i
(a)
ai w k M tn + h c i
i =1 k =1
it is possible to obtain methods of the desired order (see [5]) which, however, do not correspond to any application of an ESPRK method defined by (7) for solving (2). In the followings we propose an alternative formulation of the problem (1) that allows to reproduce the optimal performance of ESPRKN methods. Our approach is based on the main observation that, from the relations
1
(a)
M tn + h c i
hai
(a)
tn +hc i +1
+ ai s ds =
0
M (s) ds, (a)
tn +hc i
(a) tn +hc i +1
(b)
M (s) ds = (a)
M (s) ds + (a)
tn +hc i
(a)
tn +hc i +1
tn +hc i
tn +hc i
M (s) ds,
(18)
(b)
tn +hc i
it results that the method (9) performs the first order rectangular right and left quadrature approximations of the integrals in the right-hand side of (18). Similarly, the method (10), corresponds to the first order rectangular left and right quadrature approximations. Let us generalize this r approach by considering a r-step quadrature rule with weights w r and knots sr . Assuming a symmetric rule, with i =1 w r = 1 and sr −i +1 = 1 − s i for i = 1, . . . , r, the problem (6) can be formulated as the equivalent partitioned system
⎧ q˙ = 1, ⎪ ⎨ 0 r q˙ = w j M ( s j q 0 + s r − j +1 p 0 ) p , ⎪ ⎩ j =1
p˙ 0 = 1, p˙ = f (q0 , q).
(19)
We claim that Proposition 4.1. A 2r order ESPRKN method preserves its accuracy when applied to non-autonomous problems (2) written in the equivalent form of partitioned systems (19) where w j and s j are the weights and knots, respectively, of a 2r order symmetric quadrature rule in [0, 1].
838
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
(b)
Proof. Firstly, consider that the amplitudes of the integrals in (18) are h(c i since
(a)
(a)
(b)
− c i ) = hβi , and h(c i +1 − c i ) = hαi . Moreover,
(a)
(b) (a) + sr − j +1 tn + hc i = tn + hc i + hβi sr − j +1 , (a) (b) = tn + hc i(b) + hαi s j , s j tn + hc i +1 + sr − j +1 tn + hc i s j tn + hc i
it results that the application of the 2r order ESPRKN methods to the partitioned system (19) leads to a method which evolves according to the following approximations
(a)
P i = P i −1 + hb i f Q i −1 , tn + hc i Q i = Q i −1 + β i h
r
,
(a)
w j M tn + hc i
+ h β i s r − j +1 P i
j =1 r
+ αi h
(b)
w j M tn + hc i
+ h αi s j P i .
(20)
j =1
Since w j = w r − j +1 for j = 1, . . . , r, in the Q i ’s formula we recognize 2r order quadrature rules applied to the integrals in (18). Hence, the method (20), whose coefficients are provided by a 2rth order ESPRKN method, is a 2rth order approximation of the method (17), therefore preserves its accuracy. 2 5. Numerical tests The numerical examples herein are aimed at showing the preservation of the accuracy of ESPRKN methods when applied as in (20) for solving non-autonomous second order dynamics (2). To fully exploit the results given in Proposition 4.1, we consider methods defined by symmetric coefficients. That way the 2rth order methods are obtained by satisfying the conditions for the 2(r − 1)th order. In particular, we consider the optimal 11-stage symmetric fourth order ESPRKN method b S R K N 11 given in [7] which exhibits a sixth order behaviour on autonomous problems (3). By means of its coefficients (6)
we define as S R K Q N 11 the method that exploits in (20) the approximation provided by the sixth order Gauss–Legendre quadrature rule:
Q i = Q i −1 + h β i
+ h αi +
5−
5 18
√
5 18 15
√ 5 − 15 (a) (b) βi M tn + hc i + hβi + αi M tn + hc i
10
(a)
h βi
h αi
(b)
βi M tn + hc i + + αi M tn + hc i + 2 2 √ √ 5 + 15 5 + 15 (a) (b) Pi. βi M tn + hc i + hβi + αi M tn + hc i + hαi 10
+
4
9
10
10
b We will compare the results of the proposed method with the performances of S R K N 11 when applied according to the classical approach (9). We measure the 2-norm error on the numerical approximations with respect to the solution computed by the Matlab built-in function ODE45 with tolerances set at the machine precision. We will show that the order reduction we can appreciate on the classical implementation, when the effects of the time dependent matrix M (t ) are increased, can (6) be nullified by the proposed method S R K Q N 11 .
5.1. Example 1: Time dependent damped harmonic oscillator: Autonomous damping force linear in the velocity As a test problem we choose the Duffing oscillator (see [12]):
q˙ = M (t ) p ,
p˙ = − M −1 (t ) q3 − q − δ cos(ωt ) .
(21)
with M (t ) = e − t and 0 1. By varying the values of , the dependence of M (t ) on time is amplified or destroyed. Indeed, for values of ≈ 0 the problem behaves as the second order differential equation
q˙ = p ,
p˙ = q3 − q − δ cos(ωt )
that can be accurately integrated with ESPRKN methods. For increasing values of presence of the exponential term. Notice that, from the relation
d dt
˙ M −1 q˙ + M −1 q¨ M −1 q˙ = − M −1 M
their behavior is more affected by the
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
839
˙ M −1 , the problem (1) can be equivalently formulated as and by setting A = M q¨ = A (t )˙q + M (t ) f (q, t )
(22) e − t , hence it results
that is the classical formulation of a damped oscillator. In the case of the Duffing oscillator, M (t ) = that A (t ) = − . In this case, the amplification of the time dependence of M (t ) in (21) corresponds to the amplification of the ‘damped’ structure of the oscillator (22). Again, for increasing values of , ESPRKN methods feel the effects of the autonomous damping term − q˙ and they can lose the order they show when used to integrate the undamped oscillator
q¨ = q3 − q − δ cos(ωt ) that corresponds to = 0. In our test we perform the integration in the interval [0, 10π ] with initial conditions q(0) = 1.75 and p (0) = 0. In order to verify the order of the different methods, the initial step h = π4 is successively halved and the maximum error is plotted against the decreasing stepsizes. In Fig. 1(a), (b), (c) the values of the parameter δ is set at 1/4 and ω = 1, while different choices for parameter are considered. Firstly, we suppose = 1e − 4 in order to neglect the explicit time dependence of M (t ); in this case both the considered methods show a sixth order accuracy as expected and as it can be seen in Fig. 1(a). The cases corresponding to = 0.1, 0.5, are illustrated in Fig. 1(b), (c): we show both the fourth order accuracy of the approximation provided by (6) b and the theoretically predicted sixth order accuracy of the new method S R K Q N 11 . The difference in the method S R K N 11 performance of the two methods is enhanced for increasing values of the parameter. 5.2. Example 2: Time dependent damped harmonic oscillator: Non-autonomous damping force linear in the velocity We consider the one-dimensional time dependent harmonic oscillator given in [20]. The model is defined by M (t ) = t e − L (t ) where L (t ) = 0 sin(τ /π )dτ and 0 < < 1; moreover f (q, t ) = − M (t )−1 ω2 (t )q where ω(t ) = cos(t /2). The equivalent formulation (22) corresponds to
q¨ = − sin(t /π )˙q − ω2 (t )q; Also in this case, the effect of the non-autonomous damping term − sin(t /π )˙q is amplified for increasing values of . The integration is performed in the interval [0, 100π ], with initial step h = π4 and initial conditions q(0) = 1, p (0) = 0. Fig. 1(d), (e), (f) corresponds to different values of parameter : at = 0.018π that defines the model in [20], both the methods show a sixth order accuracy; when the parameter is set at = 0.1π the approximations provided by the method (6) (9) feel the effects of this increment and begin to lose their accuracy. At = 0.5π the S R K Q N 11 outperforms the method b S R K N 11 providing sixth order approximations.
5.3. Quantum damped oscillators: The Walker–Preston model The Walker–Preston model of a diatomic molecule in a strong laser field [21] corresponds to the one-dimensional Schrödinger equation
i
∂ 1 ∂2 ψ(x, t ) = − + V ( x ) + f ( t ) x ψ(x, t ) ∂t 2μ ∂ x2
where h¯ = 1 and ψ(x, 0) = ψ0 (x). In the model, the Morse potential V (x) = D (1 − e −α x )2 and f (t ) = A cos(γ t ) are considered. The parameters, in atomic units a.u, are set at the values reported in [11,19,4] i.e. μ = 1745, D = 0.2251, α = 1.1741, A = 0.011025 and γ = 0.01787. The spatial interval [−0.8, 4.32] is divided in N = 64 parts of length x = 0.08 and periodic boundary conditions are imposed. The initial condition is given by
ψ0 (x) = σ e −(Ω−0.5)α x e −Ω e
−α x
√
with Ω = 2D /ω0 , ω0 = α 2D /μ and σ a normalizing constant. After a space discretization by means of finite differences, we obtain the following system of ODEs:
i u˙ = H (t )u where H is a real time dependent symmetric matrix and the vector u (t ) ∈ C N approximates ψ(x, t ) at the nodal values xn . The vector components of the initial condition u (0) are given by un (0) = ψ0 (xn ) for n = 1, . . . N. In order to formulate the model as a problem described by (2), we use two set of N variables q, p which correspond to the real and the imaginary parts of u:
q˙ = H (t ) p , p˙ = − H (t )q.
840
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
Fig. 1. Left: Duffing oscillator. Efficiency comparison at
= 0.018π (d), = 0.1π (e), = 0.5π (f).
= 1e − 4 (a), = 0.1 (b), = 0.5 (b). Right: damped harmonic oscillator. Order comparison at
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
Fig. 2. Quantum damped oscillator,order comparison. Left:
841
γ = 0.01787. Right: γ = 0.2.
We integrate the system in the time interval [0, T ] with T = 200π and initial stepsize h = 2.5π . We measure the 2-norm error u teor ( T ) − u ( T ) , where the values of u teor are computed, using the same spatial discretization, by the function ODE45. In this as well as in the following case, since the most demanding part of the problem is the function evaluation we plot the efficiency curves (error vs. number of function evaluations). For numerical comparisons, in addition to the (6) (2) (4) proposed method S R K Q N 11 we also consider the methods S R K Q N 11 and S R K Q N 11 corresponding to the second and
b method. In the fourth order Gauss–Legendre approximations. We compare the results with the performances of the S R K N 11
(6)
Fig. 2(left), we can appreciate the sixth order accuracy of the proposed method S R K Q N 11 for the chosen set of parameters. (4)
Notice that the excellent performance of the S R K Q N 11 , which actually should be a fourth order method, last unabated for a wide range of the parameters. In order to discriminate between a fourth and a sixth order effective behaviour we have to force the frequency γ , which weigths the explicit time dependence, to a quite higher γ = 0.2. Results are shown in Fig. 2(right). 5.4. Quantum damped oscillators: The modified Caldirola–Kanai oscillator 1 − sin γ t The modified Caldirola–Kanai oscillator (see [15]) can be considered as a variable mass m(t ) = μ e system and the corresponding the Schrödinger equation (h¯ = 1) is given by
i
∂ m(t ) ∂ 2 ψ(x, t ) = − + V ( x , t ) ψ(x, t ). ∂t 2 ∂ x2
The potential V (x, t ) = 12 μω(t )2 e sin γ t x2 , is defined by means of ω(t ) = A cos(γ t ) with A = 0.001, γ = 0.01, μ = 1745. The second derivative is approximated with the central difference on the same spatial meshgrid considered in the previous example and we integrate the system for the real and imaginary parts of the vector u (t ), starting from the same initial conditions. The time integration is in the interval [0, T ] with T = 10π and large initial stepsize h = 4π . As for the previous example, the proposed method has been implemented by means of quadrature formulae of second, fourth and sixth (6) order, respectively. Fig. 3(left), confirms the sixth order accuracy of the proposed method S R K Q N 11 and the good perfor(4)
mance of the S R K Q N 11 method for the given set of parameters. In Fig. 3(right), with a different choice of the parameters ( A = 0.011025, γ = 0.2 and initial stepsize h = π ), all the methods behaves as fourth order ones, except for the method (6)
S R K Q N 11 which still preserves the sixth order, as expected. 6. Conclusions We have considered explicit symplectic partitioned Runge–Kutta methods for the numerical integration of a specific class of non-autonomous second order differential equations. As concerns splitting and composition methods applied to these problems, in [6] it has been shown that the way the explicit time dependence is treated can be critical. Herein, by performing the P -series analysis, we theoretically demonstrate that the ‘natural’ choices which come straightforward from integrating the time with one or the other set of variables p and q, introduce additional elementary differentials that increase the number of conditions to be imposed at a given order. Hence, by treating the time as two added system’s coordinates, we propose a simple strategy which allows to use ESPRK tailored on dynamics q˙ = p ; p˙ = f (q, t ), as integrators
842
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
Fig. 3. Caldirola–Kanai oscillator, order comparison. Left: A = 0.001,
γ = 0.01, initial stepsize h = 4π . Right: A = 0.011025, γ = 0.2, initial stepsize h = π . (6)
for the class of problems q˙ = M (t ) p ; p˙ = f (q, t ). The performances of the proposed sixth order method S R K Q N 11 have
b been compared with the classical integrator S R K N 11 method on the solution of several test problems.
Acknowledgements We would like to thank the referees for the positive suggestions and comments we believe helped us to greatly improve the presentation of the work as well as to clarify the main results. Appendix A In this appendix we provide the proof of Proposition 3.1 given in Section 3. Proof of Proposition 3.1. Denote with z w the tree ( z ◦ u ) ◦ v u and notice that z w = [ω1 , ω2 , . . . , ωn , u , v u ] w . Then, the elementary weight Φ related to method (7) is provided by
Φ( z w ) =
s
βi Φ2i −1 ( z w ) +
s
i =1
αi Φ2i (z w )
i =1
or, (α0 = βs+1 = 0) equivalently,
Φ( z w ) =
s +1
βi Φ2i −1 ( z w ) +
s +1
i =1
αi−1 Φ2i−2 (z w ),
(23)
i =1
wherewe set Φ0 ( z w ) = 0, Φk ( z w ) = Πk Ψk (u )Ψk ( v u ) for k = 1, . . . , 2s and Πk = 1 when z is an elementary tree, otherwise Πk = ln=1 Ψk (ωl ). From the definition of Ψk (see [13]), it follows that, for all τ ∈ T P b ,
Ψ1 (τ ) = 0, Ψ2i (τ ) =
i
βm Φ2m−1 (τ ) +
m =1
Ψ2i +1 (τ ) = Ψ2i (τ ),
i
αm Φ2m (τ ), i = 1, . . . , s,
m =1
i = 1, . . . , s − 1.
Since u , v u , ωl ∈ T P b , for l = 1, . . . , m, form (24) it follows that Φ( z w ) in (23), is given by
Φ( z w ) =
s +1 (βi + αi −1 )Π2i −1 Ψ2i −1 (u )Ψ2i −1 ( v u ) i =1
=
s +1
b i Π2i −1 Ψ2i −1 (u )Ψ2i −1 ( v u )
i =1
=
s i =1
b i +1 Π2i Ψ2i (u )Ψ2i ( v u ).
(24)
F. Diele, C. Marangi / Applied Numerical Mathematics 61 (2011) 832–843
843
As in [9], we notice that (see [18]) stopping the method (7) at the 2ith stage with 1 i s, a new ESPRK (inconsistent) method of the same form is obtained. From (24), it results that the elementary weight of the new method, say Φ (2i ) (τ ), is the elementary weight Ψ2i (τ ) associated with the 2i stage of (7). As a consequence, we reformulate Φ( z w ) as
Φ( z w ) =
s +1
b i +1 Π2i Φ (2i ) (u )Φ (2i ) ( v u ).
i =1
Since relation (16) still holds when replacing Φ(τ ) with Φ (2i ) (τ ), relation (16) follows.
2
In a similar way, we are able to demonstrate Proposition 3.2. Proof of Proposition 3.2. By adopting the previous notations, we have that
Φ( z ◦ r ) =
s i =1
b i +1 Π2i Ψ2i (r ) =
s
b i +1 Π2i Φ (2i ) (r ).
i =1
Since Φ(u ) + Φ( v ) = Φ(r ) still holds when replacing Φ(τ ) with Φ (2i ) (τ ), the result follows.
2
References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21]
L. Abia, J.M. Sanz-Serna, Partitioned Runge–Kutta methods for separable Hamiltonian problems, Math. Comput. 60 (1993) 617–634. J. Awrejcewicz, D. Sendkowski, Stick-slip chaos detection in coupled oscillators with friction, Int. J. Sol. Struct. 42 (2005) 5669–5682. S. Blanes, F. Casas, Splitting methods for non-autonomous separable dynamical systems, J. Phys. A: Math. Gen. 39 (2006) 5405–5423. S. Blanes, F. Casas, P.C. Moan, Splitting methods for non-autonomous linear systems, Int. J. Comput. Math. 84 (2007) 713–727. S. Blanes, F. Casas, A. Murua, Splitting methods in the numerical integration of non-autonomous dynamical systems, RACSAM: Rev. R. Acad. Cien. Serie A. Mat. (2010), submitted for publication. S. Blanes, F. Diele, C. Marangi, S. Ragni, Splitting and composition methods for explicit time dependence in separable dynamical systems, J. Comput. Appl. Math. 235 (2010) 646–659. S. Blanes, P.C. Moan, Practical symplectic partitioned Runge–Kutta and Runge–Kutta–Nyström methods, Appl. Math. 142 (2002) 313–330. J.C. Butcher, The Numerical Analysis of Ordinary Differential Equations. Runge–Kutta and General Linear Methods, Wiley, Chicester, 1982. M.P. Calvo, E. Hairer, Further reduction in the number of independent order conditions for sympletic, explicit partitioned Runge–Kutta and Runge– Kutta–Nyström methods, App. Num. Math. 18 (1995) 107–114. J.M. Cervero, J. Villarroel, On the quantum theory of the damped harmonic oscillator, J. Phys. A: Math. Gen. 17 (1984) 2963–2971. S. Gray, J.M. Verosky, Classical Hamiltonian structures in wave packet dynamics, J. Chem. Phys. 100 (1994) 5011–5022. J. Guckenheimer, P. Holmes, Nonlinear Oscillators, Dynamical Systems, and Bifurcations of Vector Fields, Springer-Verlag, New York, 1983. E. Hairer, Backward analysis of numerical integrators and symplectic methods, Annals of Numerical Mathematics 1 (1994) 107–132. E. Hairer, C. Lubich, G. Wanner, Geometric Numerical Integration. Structure-Preserving Algorithms for Ordinary Differential Equations, Springer-Verlag, Berlin, 2002. A.N. Ikot, I.O. Akpabio, M.I. Umo, E.E. Ituen, Quantum damped mechanical oscillator, Int. J. Optics (2010), doi:10.1155/2010/275910. P.G.L. Leach, Harmonic oscillator with variable mass, J. Phys. A: Math. Gen. 16 (1983) 3261–3269. A. Murua, On order conditions for partitioned symplectic methods, SIAM J. Numer. Anal. 34 (1997) 2204–2211. D. Okunbor, R.D. Skeel, Explicit canonical methods for Hamiltonian systems, Math. Comput. 59 (1992) 439–455. J.M. Sanz-Serna, A. Portillo, Classical numerical integrators for wave-packet dynamics, J. Chem. Phys. 104 (1996) 2349–2355. J. Struckmeier, C. Riedel, Canonical transformations and exact invariants for time-dependent Hamiltonian systems, Ann. Phys. 11 (2002) 15–38. R.B. Walker, K. Preston, Quantum versus classical dynamics in treatment of multiple photon excitation of anharmonic-oscillator, J. Chem. Phys. 67 (1977) 2017–2028.