# An Energy-Work Relationship Integration Scheme for Nonconservative Hamiltonian Systems

Authors: Fu Jingli; Chen Benyong; Xie Fengping
Journal: Research Letters in Physics
 (2008)

---

## Abstract

This letter focuses on studying a new energy-work relationship numerical integration scheme of nonconservative Hamiltonian systems. The signal-stage, multistage, and parallel composition numerical integration schemes are presented for this system. The high-order energy-work relation scheme of the system is constructed by a parallel connection ofn multistage scheme of order 2 which its order of accuracy is 2n. The connection, which is discrete analog of usual case, between the change of energy and work of nonconservative force is obtained for nonconservative Hamiltonian systems.This letter also shows that the more the stages of the schemes are, the less the error rate of the scheme is for nonconservative Hamiltonian systems. Finally, an applied example is discussed to illustrate these results.


---

## Body

## 1. Introduction

Recently, there have been a great number of studies on
the so-called geometric numerical integration scheme which preserve the structure
of systems [1–3]. Leimkuher and Reich pointed out that the geometric numerical
integrators are time-stepping methods, designed such that they exactly satisfy
conservation laws, symmetries or symplectic properties of a system of
differential equations [1]. Hairer et al. presented the symplectic integration of
separable Hamiltonian ordinary and partial differential equations. In this way,
the symplectic scheme is performed prior to the spatial step as opposed to the
standard approach of spatially discrediting the PDE to form a system of
Hamiltonian ODEs to which a symplectic integrator can be applied [2]. An energy-conserving scheme is one of such geometric numerical
integration scheme [4–8]. It is very known that a high-order scheme can be
constructed by connecting low-order scheme in series (hereafter we will call
it series composition) [1–3]. Now, the high-order energy-conserving scheme has been
constructed with method [9]. In this letter, we will present a new numerical
integration scheme, which is energy-work relation integration scheme, of nonconservative
Hamiltonian systems. This work also study that a high-order energy-work
relation scheme, which it has a structure connecting the order 2 multistage
scheme in parallel scheme (hereafter we will call it parallel composition
scheme), can be constructed by connecting low-order scheme in series.


## 2. Numerical Integration for Nonconservative Hamiltonian Systems

Let the configuration of a mechanical system is
described byn generalized coordinates qs(s=1,…,n) and n gener-alized momentums ps(s=1,…,n).
Suppose the system is subjected to n nonpotential generalized forces Qs′.
The gener-alized Hamiltonian canonical equations of system as (1)q˙s=∂H∂ps,p˙s=−∂H∂qs+Qs′,s=1,2,…,n, here, the Hamiltonian (2)H=H(q1,…,qn,p1,…,pn), which represents the total energy. The relationship
between the change of energy and the power of nonconservative force is easily
verified as (3)dHdt=∑s=1n(∂H∂qsq˙s+∂H∂psp˙s)=∑s=1n(∂H∂qs∂H∂ps−∂H∂ps∂H∂qs+∂H∂psQs')=∑s=1n∂H∂psQs'. The numerical integration is considered as the
discretization of (4)qs(tk+1)=qs(tk)+∫tktk+1∂H(q1,…,qn,p1,…,pn)∂psdt,ps(tk+1)=ps(tk)−∫tktk+1∂H(q1,…,qn,p1,…,pn)∂qsdt+Qs'Δt,tk=kΔt,k=0,1,2,…, which are obtained by integrating both sides of (1) on
the interval [tk,tk+1], where Δt is the step size.


## 3. Second-Order Schemes of Numerical Integration for Nonconservative Hamiltonian Systems

### 3.1. Single-Stage Scheme

Letpsk and qsk be the numerical approximations of ps(tk) and qs(tk),
respectively. Then a 1-stage scheme is given by(5)qsk+1=qsk+Ips1,0,psk+1=psk−Iqs1,0+Qs'Δt,s=1,2,…,n, with(6)Iqsa,b=(a−b)Δtδqsa,bμ¯qsa,bHk,Ipsa,b=(a−b)Δtδpsa,bμ¯psa,bHk,Hk=H(q1k,…,qnk,p1k,…,pnk). The notations δpsa,b and δqsa,b denote the partial difference quotient
operators with respect to ps and qs,
respectively, which are defined as(7)δqsa,bF(q1l1,…,qnln,p1k1,…,pnkn)=(Eqsa−Eqsb)F(Eqsa−Eqsb)qsls,δpsa,bF(q1l1,…,qnln,p1k1,…,pnkn)=(Epsa−Epsb)F(Epsa−Epsb)psks, where Eqsa and Epsa are the shift operators defined as(8)EqsaF(q1l1,…,qs−1ls−1,qsls,qs+1ls+1,…,qnln,p1k1,…,pnkn)=F(q1l1,…,qs−1ls−1,qsls+a,qs+1ls+1,…,qnln,p1k1,…,pnkn),EpsaF(q1l1,…,qnln,p1k1,…,ps−1ks−1,psks,ps+1ks+1,…,pnkn)=F(q1l1,…,qnln,p1k1,…,ps−1ks−1,psks+a,ps+1ks+1,…,pnkn). The notations μ¯qsa,b and μ¯psa,b denote the mean difference operators with
respect to all variables except for qs and ps,
respectively, which are defined as(9)μ¯qsa,bF(q1l1,…,qnln,p1k1,…,pnkn)=Ma,b(Eq1,…,Eqs−1,Eqs+1,…,Eqn,Ep1,…,Epn),μ¯qsa,bF(q1l1,…,qnln,p1k1,…,pnkn)=Ma,b(Eq1,…,Eqn,Ep1,…,Eps−1,Eps+1,…,Epn), with(10)Ma,b(x1,x2,…,xr−1)=1r!∑l=1rper(x1ax2a⋯xr−1a⋮⋮⋮x1ax2a⋯xr−1ax1bx2b⋯xr−1b⋮⋮⋮x1bx2b⋯xr−1b)}r−1}l−1, where per(A)
denotes the permanent or plus determinant of a matrix A [10]. For example, in the case d = 1, we have(11)μ¯q1a,b=Ma,b(Eq1)=12(Eq1a+Eq1b),μ¯p1a,b=Ma,b(Ep1)=12(Ep1a+Ep1b). The operators δqsa,b,δpsa,b,μ¯qsa,b,
and μ¯psa,b have symmetry expressed as(12)δqsa,b=δqsb,a,δpsa,b=δpsb,a,μ¯qsa,b=μ¯qsb,a,μ¯psa,b=μ¯psb,a.


### 3.2. Relation between the Energy and Work  of Nonconservative Force for Nonconservative Hamiltonian Systems

Proposition 1.
The relation between the energy and work of nonconservative
force for nonconservative Hamiltonian system holds:
(13)Hk+a−Hk+b=∑s=1n[(δqsa,bμ¯qsa,bHk)(qsk+a−qsk+b)+(δpsa,bμ¯psa,bHk)(psk+a−psk+b)].Proof.
For simplicity, we set
(14)E¯a=(Eq1a,…,Eqna,Ep1a,…,Epna)T. We first note the identity
(15)∏snEqsaEpsa=1(2n)!per(E→a,…,E→a︸2n). It follows that (16)Hk+a−Hk+b=(∏s=1nEqsaEpsa−∏s=1nEqsbEpsb)Hk=∑s=1n[(Eqsa−Eqsb)Ma,b(Eq1,…,Eqn,Ep1,…,Eps−1,Eps+1,…,Epn)+(Epsa−Epsb)Ma,b(Eq1,…,Eqs−1,Eqs+1,…,Eqn,Ep1,…,Epn)]Hk=∑s=1n[(δqsa,bμ¯qsa,bHk)(qsk+a−qsk+b)+(δpsa,bμ¯psa,bHk)(psk+a−psk+b)], where we have used the properties of the permanent and
the definitions  of operators (7)–(10) [10].Proposition 2.
The scheme{(1),…,(6)} satisfies the relation between the change of energy
and the work of nonconservative force for the system.Proof.
We see from the chain rule (13) that the change
of energy Hk is equivalent to the work of nonconservative
force:
(17)Hk+1−Hk=∑s=1n[(δqs1,0μ¯qs1,0Hk)(qsk+1−qsk)+(δps1,0μ¯ps1,0Hk)(psk+1−psk)]=Δt∑s=1n[(δqs1,0μ¯qs1,0Hk)(δps1,0μ¯ps1,0Hk)+(δps1,0μ¯ps1,0Hk)(−δqs1,0μ¯qs1,0Hk+Qs')]=∑s=1nΔt(δps1,0μ¯ps1,0Hk)Qs', which is a discrete analog of (3).


### 3.3. Order of Accuracy

The local errors involved in the determination of{psk+1,qsk+1}s=1n from {psk,qsk}s=1n are O(Δt3),
that is, Iqs1,0 and Ips1,0 in the scheme (5) are the second-order
approximations of the integrals in (4), respectively. Although this can be
proved by the Taylor
expansions, it is obvious because the scheme is symmetric (see Section 4.3).


### 3.4. Multistage Scheme

Anc-stage
scheme is constructed by connecting the second-order scheme with small
integration interval of length Δt/c in series:(18)Psk+m/c=Psk+(m−1)/c−IQsm/c,(m−1)/c+Qs'(Q1k,…,Qnk,P1k,…,Pnk)Δtc,Qsk+m/c=Qsk+(m−1)/c+IPsm/c,(m−1)/cPsk+1psk+1,Psk=psk,Qsk+1=qsk+1,Qsk=qsk,s=1,2,…,n,m=1,2,…,c, with(19)IQsa,b=(a−b)ΔtδQsa,bμ¯Qsa,bH(Q1k,…,Qnk,P1k,…,Pnk),IPsa,b=(a−b)ΔtδPsa,bμ¯Psa,bH(Q1k,…,Qnk,P1k,…,Pnk), where Psk+m/c and Qsk+m/c are the internal  stage variables. It should be
noted that the above scheme is equivalent to the scheme:(20)psk+1=psk−∑l=1cIQsl/c,(l−1)/c+Qs'(q1k,…,qnk,p1k,…,pnk)Δtc,qsk+1=qsk−∑l=1cIPsl/c,(l−1)/c,Psk+m/c=c−mc(psk−∑l=1mIQsl/c,(l−1)/c+Qs'Δtc)+mc(psk+1+∑l=m+1cIQsl/c,(l−1)/c+Qs'Δtc),Qsk+m/c=c−mc(qsk+∑l=1mIPsl/c,(l−1)/c)+mc(qsk+1−∑l=m+1cIPsl/c,(l−1)/c),s=1,2,…,n,m=1,2,…,c−1. The latter scheme (20) will be used in the next
section to construct a higher-order scheme.It is obvious for thec-stage that the relationship
between the change of energy and the work of nonconservative force is exactly equivalent
and that the order of accuracy is 2. We point out here that the local error is expressed as c×O[(Δt/c)3]=c−2O(Δt)3.


## 3.1. Single-Stage Scheme

Letpsk and qsk be the numerical approximations of ps(tk) and qs(tk),
respectively. Then a 1-stage scheme is given by(5)qsk+1=qsk+Ips1,0,psk+1=psk−Iqs1,0+Qs'Δt,s=1,2,…,n, with(6)Iqsa,b=(a−b)Δtδqsa,bμ¯qsa,bHk,Ipsa,b=(a−b)Δtδpsa,bμ¯psa,bHk,Hk=H(q1k,…,qnk,p1k,…,pnk). The notations δpsa,b and δqsa,b denote the partial difference quotient
operators with respect to ps and qs,
respectively, which are defined as(7)δqsa,bF(q1l1,…,qnln,p1k1,…,pnkn)=(Eqsa−Eqsb)F(Eqsa−Eqsb)qsls,δpsa,bF(q1l1,…,qnln,p1k1,…,pnkn)=(Epsa−Epsb)F(Epsa−Epsb)psks, where Eqsa and Epsa are the shift operators defined as(8)EqsaF(q1l1,…,qs−1ls−1,qsls,qs+1ls+1,…,qnln,p1k1,…,pnkn)=F(q1l1,…,qs−1ls−1,qsls+a,qs+1ls+1,…,qnln,p1k1,…,pnkn),EpsaF(q1l1,…,qnln,p1k1,…,ps−1ks−1,psks,ps+1ks+1,…,pnkn)=F(q1l1,…,qnln,p1k1,…,ps−1ks−1,psks+a,ps+1ks+1,…,pnkn). The notations μ¯qsa,b and μ¯psa,b denote the mean difference operators with
respect to all variables except for qs and ps,
respectively, which are defined as(9)μ¯qsa,bF(q1l1,…,qnln,p1k1,…,pnkn)=Ma,b(Eq1,…,Eqs−1,Eqs+1,…,Eqn,Ep1,…,Epn),μ¯qsa,bF(q1l1,…,qnln,p1k1,…,pnkn)=Ma,b(Eq1,…,Eqn,Ep1,…,Eps−1,Eps+1,…,Epn), with(10)Ma,b(x1,x2,…,xr−1)=1r!∑l=1rper(x1ax2a⋯xr−1a⋮⋮⋮x1ax2a⋯xr−1ax1bx2b⋯xr−1b⋮⋮⋮x1bx2b⋯xr−1b)}r−1}l−1, where per(A)
denotes the permanent or plus determinant of a matrix A [10]. For example, in the case d = 1, we have(11)μ¯q1a,b=Ma,b(Eq1)=12(Eq1a+Eq1b),μ¯p1a,b=Ma,b(Ep1)=12(Ep1a+Ep1b). The operators δqsa,b,δpsa,b,μ¯qsa,b,
and μ¯psa,b have symmetry expressed as(12)δqsa,b=δqsb,a,δpsa,b=δpsb,a,μ¯qsa,b=μ¯qsb,a,μ¯psa,b=μ¯psb,a.


## 3.2. Relation between the Energy and Work  of Nonconservative Force for Nonconservative Hamiltonian Systems

Proposition 1.
The relation between the energy and work of nonconservative
force for nonconservative Hamiltonian system holds:
(13)Hk+a−Hk+b=∑s=1n[(δqsa,bμ¯qsa,bHk)(qsk+a−qsk+b)+(δpsa,bμ¯psa,bHk)(psk+a−psk+b)].Proof.
For simplicity, we set
(14)E¯a=(Eq1a,…,Eqna,Ep1a,…,Epna)T. We first note the identity
(15)∏snEqsaEpsa=1(2n)!per(E→a,…,E→a︸2n). It follows that (16)Hk+a−Hk+b=(∏s=1nEqsaEpsa−∏s=1nEqsbEpsb)Hk=∑s=1n[(Eqsa−Eqsb)Ma,b(Eq1,…,Eqn,Ep1,…,Eps−1,Eps+1,…,Epn)+(Epsa−Epsb)Ma,b(Eq1,…,Eqs−1,Eqs+1,…,Eqn,Ep1,…,Epn)]Hk=∑s=1n[(δqsa,bμ¯qsa,bHk)(qsk+a−qsk+b)+(δpsa,bμ¯psa,bHk)(psk+a−psk+b)], where we have used the properties of the permanent and
the definitions  of operators (7)–(10) [10].Proposition 2.
The scheme{(1),…,(6)} satisfies the relation between the change of energy
and the work of nonconservative force for the system.Proof.
We see from the chain rule (13) that the change
of energy Hk is equivalent to the work of nonconservative
force:
(17)Hk+1−Hk=∑s=1n[(δqs1,0μ¯qs1,0Hk)(qsk+1−qsk)+(δps1,0μ¯ps1,0Hk)(psk+1−psk)]=Δt∑s=1n[(δqs1,0μ¯qs1,0Hk)(δps1,0μ¯ps1,0Hk)+(δps1,0μ¯ps1,0Hk)(−δqs1,0μ¯qs1,0Hk+Qs')]=∑s=1nΔt(δps1,0μ¯ps1,0Hk)Qs', which is a discrete analog of (3).


## 3.3. Order of Accuracy

The local errors involved in the determination of{psk+1,qsk+1}s=1n from {psk,qsk}s=1n are O(Δt3),
that is, Iqs1,0 and Ips1,0 in the scheme (5) are the second-order
approximations of the integrals in (4), respectively. Although this can be
proved by the Taylor
expansions, it is obvious because the scheme is symmetric (see Section 4.3).


## 3.4. Multistage Scheme

Anc-stage
scheme is constructed by connecting the second-order scheme with small
integration interval of length Δt/c in series:(18)Psk+m/c=Psk+(m−1)/c−IQsm/c,(m−1)/c+Qs'(Q1k,…,Qnk,P1k,…,Pnk)Δtc,Qsk+m/c=Qsk+(m−1)/c+IPsm/c,(m−1)/cPsk+1psk+1,Psk=psk,Qsk+1=qsk+1,Qsk=qsk,s=1,2,…,n,m=1,2,…,c, with(19)IQsa,b=(a−b)ΔtδQsa,bμ¯Qsa,bH(Q1k,…,Qnk,P1k,…,Pnk),IPsa,b=(a−b)ΔtδPsa,bμ¯Psa,bH(Q1k,…,Qnk,P1k,…,Pnk), where Psk+m/c and Qsk+m/c are the internal  stage variables. It should be
noted that the above scheme is equivalent to the scheme:(20)psk+1=psk−∑l=1cIQsl/c,(l−1)/c+Qs'(q1k,…,qnk,p1k,…,pnk)Δtc,qsk+1=qsk−∑l=1cIPsl/c,(l−1)/c,Psk+m/c=c−mc(psk−∑l=1mIQsl/c,(l−1)/c+Qs'Δtc)+mc(psk+1+∑l=m+1cIQsl/c,(l−1)/c+Qs'Δtc),Qsk+m/c=c−mc(qsk+∑l=1mIPsl/c,(l−1)/c)+mc(qsk+1−∑l=m+1cIPsl/c,(l−1)/c),s=1,2,…,n,m=1,2,…,c−1. The latter scheme (20) will be used in the next
section to construct a higher-order scheme.It is obvious for thec-stage that the relationship
between the change of energy and the work of nonconservative force is exactly equivalent
and that the order of accuracy is 2. We point out here that the local error is expressed as c×O[(Δt/c)3]=c−2O(Δt)3.


## 4. Higher-Order Schemes of Numerical Integration for Nonconservative Hamiltonian Systems

### 4.1. Parallel Composition Scheme

Letc1,c2,…,cn be arbitrary positive integers satisfying(21)c1<c2<⋯<cn, then a new scheme is constructed by connection c1-stage,c2-stage,…,cn-stage schemes of order 2 in parallel: (22)psk+1=psk−∑j=1udj∑l=1cjIQsjl/cj,(l−1)/cj+Qs'(q1k,…,qnk,p1k,…,pnk)Δt,qsk+1=qsk+∑j=1udj∑l=1cjIPijl/cj,(l−1)/cj,Psjk+m/cj=cj−mcj(psk−∑l=1mIQsjl/cj,(l−1)/cj+Qsj'Δt)+mcj(psk+1−∑l=m+1mIQsjl/cj,(l−1)/cj+Qsj'Δt),Qsjk+m/cj=cj−mcj(qsk+∑l=1mIQsjl/cj,(l−1)/cj)+mcj(qsk+1+∑l=m+1mIQsjl/cj,(l−1)/cj),Psjk+1=psk+1,Psjk=psk,Qsjk+1=qsk+1,Qsjk=qsk,s=1,2,…,n,j=1,2,…,u,m=1,2,…,cj−1, with the weights(23)dj={1foru=1,cj2u−2∏l=1,l≠ju(cj2−cl2)foru≥2,j=1,2,…,u, where(24)IPsja,b=(a−b)ΔtδPsja,bμ¯Psja,bHjk,IQsja,b=(a−b)ΔtδQsja,bμ¯Qsja,bHjk,Hjk=H(Q1jk,…,Qnjk,P1jk,…,Pnjk),Λsjk=Λs(Q1jk,…,Qnjk,P1jk,…,Pnjk).


### 4.2. Relation between the Change of Energy  and Work of Nonconservative Force

Proposition 3.
The scheme (23) with the condition
(25)∑j=1udj=1 satisfies relationship between the change of energy
and work of nonconservative force for nonconservative Hamiltonian systemsProof.
We first note
(26)Hk+1=Hjk+1,Hk=Hjk,j=1,2,…,u. We see from Proposition 1 that
(27)Hjk+a−Hjk+b=∑s=1n(δPsja,bμ¯Psja,bHjk)(Psjk+a−Psjk+b)+(δQsja,bμ¯Qsja,bHjk)(Qsjk+a−Qsjk+b). It follows from (25)–(27) that
(28)Hk+1−Hk=∑j=1udj(Hjk+1−Hjk)=∑j=1udj∑m=1cj(Hjk+m/cj−Hjk+(m−1)/cj)=∑j=1udjcjΔt∑m=1cj∑s=1n[IPsjk+m/cj,k+(m−1)/cj×(Psjk+m/cj−Psjk+(m−1)/cj)+IQsjk+m/cj,k+(m−1)/cj×(Qsjk+m/cj−Qsjk+(m−1)/cj)]. We obtain from (22)
(29)Psjk+m/cj−Psjk+(m−1)/cj=1cj(∑r=1ucr∑l=1cjIPsrl/cr,(l−1)/cr−∑l=1cjIPsjl/cj,(l−1)/cj)+IPsjm/cj,(m−1)/cj Substituting (29) into (28) yields
(30)Hk+1−Hk=∑j=1udj∑m=1cjcj∑s=1nIPsjm/cj,(m−1)/cjQsj', which is a discrete analog of that relation between
the change of energy and work of nonconservative force for the systems (17).


### 4.3. Order of a Symmetric Scheme

Proposition 4.
Consider the scheme (22) as
mapping
(31)ϕΔt:(q1k,…,qnk,p1k,…,pnk)=(q1k+1,…,qnk+1,p1k+1,…,pnk+1), and let ϕΔt−1 be the inverse mapping of ϕΔt. Then,
one has
(32)ϕ−Δt−1=ϕΔt. That is, the scheme is symmetricProof.
The inverseϕΔt−1 is obtained by exchanging (psk,qsk) and (psk+1,qsk+1).
Replacing Δt by −Δt and rearranging terms in ϕΔt−1 leads to the mapping ϕ−Δt−1.
For this ϕ−Δt−1,
setting
(33)P¯sjk+m/cj=Psjk+1−m/cj,s=1,2,…,n,j=1,2,…,u,m=0,1,…,cj, and omitting the tilde, we can obtain ϕΔt. Therefore,
form (31) holds.Proposition 5.
If one chooses the weightsd1,d2,…,dn as (23), the accuracy of the scheme (22) is at
least of order 2n.Proof.
It is known that if a one-step scheme is
symmetric, its order of accuracy is even [1, 2]. Therefore, the local error of
the scheme ϕΔt is O(Δt2r+1) with a positive integer r. We first choose {dj}j=1u such that
(34)∑j=1udj=1. Since the error of IPsjm/cj,(m−1)/cj, IQsjm/cj,(m−1)/cj,
and Qsj'k are O[(Δt/cj)3],
the error of ϕΔt is expressed as
(35)∑j=1udj×cj×O[(Δtcj)3]=∑j=1ucj−2djO(Δt3). If one chooses {dj}j=1u such that
(36)∑j=1ncj−2dj=0, then the O(Δt3)-term in the error of ϕΔt vanishes. Since the error of ϕΔt is of odd order, it becomes O(Δt5).
The O(Δt5)-term in the error {dj}j=1u-term in the error of ϕΔt is expressed as
(37)∑j=1udj×cj×O[(Δtcj)5]=∑j=1ucj−4djO(Δt5).  These procedures can be repeated. The final condition for {dj}j=1u is
(38)∑j=1ucj−2(u−1)dj=0. Therefore, if one chooses {dj}j=1u such that they satisfy the n simultaneous linear equations:
(39)∑j=1ucj−2ldj={1forl=0,0forl=1,2,…,u−1, then the error of ϕΔt is O(Δt2u+1).
Since the solution of (39) is given by (23), the order of accuracy is 2u.


## 4.1. Parallel Composition Scheme

Letc1,c2,…,cn be arbitrary positive integers satisfying(21)c1<c2<⋯<cn, then a new scheme is constructed by connection c1-stage,c2-stage,…,cn-stage schemes of order 2 in parallel: (22)psk+1=psk−∑j=1udj∑l=1cjIQsjl/cj,(l−1)/cj+Qs'(q1k,…,qnk,p1k,…,pnk)Δt,qsk+1=qsk+∑j=1udj∑l=1cjIPijl/cj,(l−1)/cj,Psjk+m/cj=cj−mcj(psk−∑l=1mIQsjl/cj,(l−1)/cj+Qsj'Δt)+mcj(psk+1−∑l=m+1mIQsjl/cj,(l−1)/cj+Qsj'Δt),Qsjk+m/cj=cj−mcj(qsk+∑l=1mIQsjl/cj,(l−1)/cj)+mcj(qsk+1+∑l=m+1mIQsjl/cj,(l−1)/cj),Psjk+1=psk+1,Psjk=psk,Qsjk+1=qsk+1,Qsjk=qsk,s=1,2,…,n,j=1,2,…,u,m=1,2,…,cj−1, with the weights(23)dj={1foru=1,cj2u−2∏l=1,l≠ju(cj2−cl2)foru≥2,j=1,2,…,u, where(24)IPsja,b=(a−b)ΔtδPsja,bμ¯Psja,bHjk,IQsja,b=(a−b)ΔtδQsja,bμ¯Qsja,bHjk,Hjk=H(Q1jk,…,Qnjk,P1jk,…,Pnjk),Λsjk=Λs(Q1jk,…,Qnjk,P1jk,…,Pnjk).


## 4.2. Relation between the Change of Energy  and Work of Nonconservative Force

Proposition 3.
The scheme (23) with the condition
(25)∑j=1udj=1 satisfies relationship between the change of energy
and work of nonconservative force for nonconservative Hamiltonian systemsProof.
We first note
(26)Hk+1=Hjk+1,Hk=Hjk,j=1,2,…,u. We see from Proposition 1 that
(27)Hjk+a−Hjk+b=∑s=1n(δPsja,bμ¯Psja,bHjk)(Psjk+a−Psjk+b)+(δQsja,bμ¯Qsja,bHjk)(Qsjk+a−Qsjk+b). It follows from (25)–(27) that
(28)Hk+1−Hk=∑j=1udj(Hjk+1−Hjk)=∑j=1udj∑m=1cj(Hjk+m/cj−Hjk+(m−1)/cj)=∑j=1udjcjΔt∑m=1cj∑s=1n[IPsjk+m/cj,k+(m−1)/cj×(Psjk+m/cj−Psjk+(m−1)/cj)+IQsjk+m/cj,k+(m−1)/cj×(Qsjk+m/cj−Qsjk+(m−1)/cj)]. We obtain from (22)
(29)Psjk+m/cj−Psjk+(m−1)/cj=1cj(∑r=1ucr∑l=1cjIPsrl/cr,(l−1)/cr−∑l=1cjIPsjl/cj,(l−1)/cj)+IPsjm/cj,(m−1)/cj Substituting (29) into (28) yields
(30)Hk+1−Hk=∑j=1udj∑m=1cjcj∑s=1nIPsjm/cj,(m−1)/cjQsj', which is a discrete analog of that relation between
the change of energy and work of nonconservative force for the systems (17).


## 4.3. Order of a Symmetric Scheme

Proposition 4.
Consider the scheme (22) as
mapping
(31)ϕΔt:(q1k,…,qnk,p1k,…,pnk)=(q1k+1,…,qnk+1,p1k+1,…,pnk+1), and let ϕΔt−1 be the inverse mapping of ϕΔt. Then,
one has
(32)ϕ−Δt−1=ϕΔt. That is, the scheme is symmetricProof.
The inverseϕΔt−1 is obtained by exchanging (psk,qsk) and (psk+1,qsk+1).
Replacing Δt by −Δt and rearranging terms in ϕΔt−1 leads to the mapping ϕ−Δt−1.
For this ϕ−Δt−1,
setting
(33)P¯sjk+m/cj=Psjk+1−m/cj,s=1,2,…,n,j=1,2,…,u,m=0,1,…,cj, and omitting the tilde, we can obtain ϕΔt. Therefore,
form (31) holds.Proposition 5.
If one chooses the weightsd1,d2,…,dn as (23), the accuracy of the scheme (22) is at
least of order 2n.Proof.
It is known that if a one-step scheme is
symmetric, its order of accuracy is even [1, 2]. Therefore, the local error of
the scheme ϕΔt is O(Δt2r+1) with a positive integer r. We first choose {dj}j=1u such that
(34)∑j=1udj=1. Since the error of IPsjm/cj,(m−1)/cj, IQsjm/cj,(m−1)/cj,
and Qsj'k are O[(Δt/cj)3],
the error of ϕΔt is expressed as
(35)∑j=1udj×cj×O[(Δtcj)3]=∑j=1ucj−2djO(Δt3). If one chooses {dj}j=1u such that
(36)∑j=1ncj−2dj=0, then the O(Δt3)-term in the error of ϕΔt vanishes. Since the error of ϕΔt is of odd order, it becomes O(Δt5).
The O(Δt5)-term in the error {dj}j=1u-term in the error of ϕΔt is expressed as
(37)∑j=1udj×cj×O[(Δtcj)5]=∑j=1ucj−4djO(Δt5).  These procedures can be repeated. The final condition for {dj}j=1u is
(38)∑j=1ucj−2(u−1)dj=0. Therefore, if one chooses {dj}j=1u such that they satisfy the n simultaneous linear equations:
(39)∑j=1ucj−2ldj={1forl=0,0forl=1,2,…,u−1, then the error of ϕΔt is O(Δt2u+1).
Since the solution of (39) is given by (23), the order of accuracy is 2u.


## 5. A Numerical Example

Consider the motion of a particle
with unit mass whose Hamiltonian is(40)H(q,p)=12p2, and the motion of the system is
subjected to nonpotential force(41)Qs'=cosbt, where b is a constant.The
equation of motion of nonconservative particle is(42)dpdt=cosbt,dqdt=p. The analytic solution of (42) is
given by(43)p=Absinbt,q=−Ab2cos(bt+φ), which have the period T=2π/b.
We take the initial conditions:(44)p(π2b)=1,q(π2)=1, and the calculation time t=T.
The parallel composition scheme with(45)cj=j,j=1,2,…,u was used. We calculated the
global error given by(46)e(t)=(pK−p(T))2+(qK−q(T))2, where K=T/Δt.
Since the global error e(t) is about T/Δt times the local error, e(t) is expressed as e(T)=O(Δt2n).We should
point out that the local error of the parallel composition is expressed as(47)∑j=1ucj−2udjO(Δt2n+1)=1c12c22⋯cn2O(Δt2n+1),  the more the stages of the schemes are,
the smaller the error of the scheme for nonconservative Hamiltonian systems.


## 6. Conclusion

In this paper, the new numerical integration schemes
of nonconservative Hamilton
systems are established. This study has given that the numerical connection
between energy of system and work of nonconservative force is an analog of
usual energy-work connection, and the numerical connection between the
high-order energy-work is also contented. Numerical results showed  that the more the stages of the schemes are, the smaller the error of the scheme for nonconservative Hamiltonian systems.



---

*Source: 168231-2008-07-14.xml*