# New Stable Closed Newton-Cotes Trigonometrically Fitted Formulae for Long-Time Integration

Authors: T. E. Simos
Journal: Abstract and Applied Analysis
 (2012)

---

## Abstract

The closed Newton-Cotes differential methods of high algebraic order for small number of function evaluations are unstable. In this work, we propose a new closed Newton-Cotes trigonometrically fitted differential method of high algebraic order which gives much more efficient results than the well-know ones.


---

## Body

## 1. Introduction

In the recent years, there is a great interest in the construction of numerical methods for ordinary differential equations that preserve qualitative properties of the analytic solution.Symplectic integrators are necessary in the case that we wish to preserve the characteristics of the Hamiltonian system in the approximate solution. Much research has been done recently mainly on the development of one-step symplectic integrators (see [1, 2]). In their work, Zhu et. al [3] and Chiou and Wu [4] constructed multistep symplectic integrators by writing open Newton-Cotes differential schemes as multilayer symplectic structures.Last decades much work has been done on trigonometrically fitting and the numerical solution of periodic initial value problems (see [5–20] and references therein).In this paper, we follow the steps described below.(i)
The new condition is described.(ii)
The trigonometrically fitted method is developed.(iii)
The closed Newton-Cotes differential methods are presented as multilayer symplectic integrators.(iv)
The closed Newton-Cotes methods are applied to nonlinear problems and the efficiency of the new methods is presented.We note that the aim of this paper is to generate methods that can be used for nonlinear differential equations as well as linear ones.The construction of the paper is given below.(i)
The theory for the symplectic schemes is presented in Section2.(ii)
In Section3, we present the closed Newton-Cotes differential methods and the new condition for the development of the methods. We also develop the new trigonometrically-fitted methods.(iii)
In Section4, the conversion of the closed Newton-Cotes differential methods into multilayer symplectic structures is presented.(iv)
Numerical results are presented in Section5.


## 2. Basic Theory on Symplectic Schemes and Numerical Methods

Based on Zhu et al. [3] and on the division of the interval [a,b] with N points, we have the following discrete scheme for the n-step approximation to the solution:(2.1)(pn+1qn+1)=Mn+1(pnqn),Mn+1=(wn+1yn+1zn+1gn+1).
Based on the above we can write the n-step approximation to the solution as(2.2)(pnqn)=(wnynzngn)(wn-1yn-1zn-1gn-1)⋯(w1y1z1g1)(p0q0)=MnMn-1⋯M1(p0q0).
Defining(2.3)S=MnMn-1⋯M1=(WnYnZnGn),
the discrete transformation can be written as(2.4)(pnqn)=S(p0q0).A discrete scheme (2.1) is a symplectic scheme if the transformation matrix S is symplectic. A matrix A is symplectic if ATJA=J, where(2.5)J=(01-10).
The product of symplectic matrices is also symplectic. Hence, if each matrix Mn is symplectic the transformation matrix S is symplectic. Consequently, the discrete scheme (3.5) is symplectic if each matrix Mn is symplectic.Remark 2.1.
The proposed methods can be used for nonlinear differential equations as well as linear ones.


## 3. Trigonometrically Fitted Closed Newton-Cotes Differential Methods

### 3.1. General Closed Newton-Cotes Formulae

The closed Newton-Cotes integral rules can be presented with the formula:(3.1)∫abf(x)dx≈zh∑i=0ktif(xi),
where(3.2)h=b-aN,xi=a+ih,i=0,1,2,…,N.
The coefficient z as well as the weights ti are given in Table 1.Table 1
Closed Newton-Cotes integral rules.
kzt0t1t2t3t4t5t6t7t801111/21121/314133/8133142/457321232755/28819755050751961/1404121627272272164177/1728075135771323298929891323357775184/141759895888−92810496−454010496−9285888989Remark 3.1.
It is easy for one to see that the coefficientsti in the Table 1 are symmetric, that is, one has the following relation:
(3.3)ti=tk-i,i=0,1,…,k2.The closed Newton-Cotes differential methods are produced from the integral rules. From Table1 we have the following differential methods: (3.4)k=1yn+1-yn=h2(fn+1+fn),k=2yn+1-yn-1=h3(fn-1+4fn+fn+1),k=3yn+1-yn-2=3h8(fn-2+3fn-1+3fn+fn+1),k=4yn+2-yn-2=2h45(7fn-2+32fn-1+12fn+32fn+1+7fn+1),k=5yn+2-yn-3=5h288(19fn-3+75fn-2+50fn-1+50fn+75fn+1+19fn+2),k=6yn+3-yn-3=h140(41fn-3+216fn-2+27fn-1+272fn+27fn+1+216fn+2+41fn+3),k=7yn+3-yn-4=7h17280(751fn-4+3577fn-3+1323fn-2+2989fn-1+2989fn+1323fn+1+3577fn+2+751fn+3),k=8yn+4-yn-4=4h14175(989fn-4+5888fn-3-928fn-2+10496fn-1-4540fn+10496fn+1-928fn+2+5888fn+3+989fn+4).In the present paper, we will investigate the casek=8 and we will produce trigonometrically fitted differential methods of order 1.


### 3.2. Development of Closed Newton-Cotes Differential Schemes

For the development of a Newton-Cotes differential method of the above form, two procedures can be applied.(i)
The procedure which is based on the minimization of the local truncation error. Based on this procedure and for the casek=8, we can produce the well known coefficients:  a4=3956/14175, a3=23552/14175, a2=-3712/14175, a1=41984/14175, and a0=-3632/2835 (see the closed Newton-Cotes differential scheme for k=8 presented above).(ii)
The procedure which is based on
(1)
the minimization of the local truncation error;(2)
the satisfaction of the condition:2(∑i=1n|ai|+|a0|)/2n=1. This condition is produced by application of least squares method to the production of the differential method (see more details in [21] (about stable quadrature rules) and [22]).The above procedure leads to the following coefficients for the case:n=4: a0=-19672/945+70a4, a1=1952/105-56a4, a2=-848/105+28a4, a3=736/189-8a4 and to the condition 9836/33075≤a4≤244/735. We choose the value: a4=3/10, which satisfies the above condition.


### 3.3. Exponentially Fitted Closed Newton-Cotes Differential Method

Requiring the differential scheme:(3.5)yn+4-yn-4=h(a4fn-4+a3fn-3+a2fn-2+a1fn-1+a0fn+a1fn+1+a2fn+2+a3fn+3+a4fn+4)
to be accurate for the following set of functions (we note that fi=yi′,i=n-1,n,n+1): (3.6){1,x,x2,x3,x4,x5,x6,cos⁡(vx),sin⁡(vx)},
the following set of equations is obtained:(3.7)8cos⁡(w)sin⁡(w)(2(cos⁡(w))2-1)=w[(cos⁡(w))42a4-2a2+a0-6a3cos⁡(w)+4a2(cos⁡(w))2+2a1cos⁡(w)+8a3(cos⁡(w))3+16a4(cos⁡(w))4-16a4(cos⁡(w))2]2a4+2a3+2a2+2a1+a0=896a4+54a3+24a2+6a1=1282560a4+810a3+160a2+10a1=2048.
Requesting that a4=3/10 and solving the above system of equations, we obtain(3.8)a0=1735wcos⁡(w)+1353wcos⁡(3w)+242wcos⁡(2w)+270wcos⁡(4w)-900sin⁡(4w)675wcos⁡(w)+45wcos⁡(3w)-450w-270wcos⁡(2w),a1=-1856wcos⁡(3w)-1735w-1404wcos⁡(2w)-405wcos⁡(4w)+1350sin⁡(4w)1350wcos⁡(w)+90wcos⁡(3w)-900w-540wcos⁡(2w),a2=702wcos⁡(w)+418wcos⁡(3w)-121w+81wcos⁡(4w)-270sin⁡(4w)675wcos⁡(w)+45wcos⁡(3w)-450w-270wcos⁡(2w),a3=-27wcos⁡(4w)-836wcos⁡(2w)-1353w+1856wcos⁡(w)+90sin⁡(4w)1350wcos⁡(w)+90wcos⁡(3w)-900w-540wcos⁡(2w),
where w=vh.For small values ofv, the above formulae are subject to heavy cancellations. In this case the following Taylor series expansions must be used:(3.9)a0=46613780-40735670w2+40193249480w4-6980443681080400w6+945598949037788800w8-195517098336424096000w10+694578131900704693888000w12+33104793794390627842881280000w14+582548167731615751046180311040000w16+9723295174775617148961238556672000w18+⋯,a1=499525+40737560w2-40193332640w4+6980443908107200w6-945598965383718400w8+1955170911115232128000w10-694578132534272925184000w12-33104793795854170457175040000w14-582548167732154334728240414720000w16-97232951747100822865281651408896000w18+⋯,a2=7811050-407318900w2+40193831600w4-69804432270268000w6+9455989163459296000w8-1955170927788080320000w10+694578136335682312960000w12+331047937914635426142937600000w14+582548167735385836820601036800000w16+97232951747252057163204128522240000w18+⋯,a3=64934725+4073113400w2-401934989600w4+698044313621608000w6-9455989980755776000w8+19551709166728481920000w10-6945781338014093877760000w12-331047937987812556857625600000w14-5825481677332315020923606220800000w16-972329517471512342979224771133440000w18+⋯.The behaviour of the coefficients is given in the following Figure1.Behavior of the coefficients of the new proposed method given by (3.8) for several values of w.
(a)(b)(c)(d)The local truncation error for the above differential method is given by(3.10)L.T.E(h)=-593h928350(yn(9)+v2yn(7)).
The L.T.E is obtained expanding the terms yn±j and fn±j,j=1(1)4 in (3.5) into Taylor series expansions and substituting the Taylor series expansions of the coefficients of the method.In Figure2, we present the behaviour of the quantity ST=(2∑i=1n|ai|+|a0|)/2n for several values of v.Figure 2
Behaviour of the quantity ST for several values ofv.So, we have the following theorem.Theorem 3.2.
The method (3.5) with coefficients ai,i=0(1)4, obtained by the solution of the system (3.7) is accurate for the set of functions (3.6) and is of eighth algebraic order.


## 3.1. General Closed Newton-Cotes Formulae

The closed Newton-Cotes integral rules can be presented with the formula:(3.1)∫abf(x)dx≈zh∑i=0ktif(xi),
where(3.2)h=b-aN,xi=a+ih,i=0,1,2,…,N.
The coefficient z as well as the weights ti are given in Table 1.Table 1
Closed Newton-Cotes integral rules.
kzt0t1t2t3t4t5t6t7t801111/21121/314133/8133142/457321232755/28819755050751961/1404121627272272164177/1728075135771323298929891323357775184/141759895888−92810496−454010496−9285888989Remark 3.1.
It is easy for one to see that the coefficientsti in the Table 1 are symmetric, that is, one has the following relation:
(3.3)ti=tk-i,i=0,1,…,k2.The closed Newton-Cotes differential methods are produced from the integral rules. From Table1 we have the following differential methods: (3.4)k=1yn+1-yn=h2(fn+1+fn),k=2yn+1-yn-1=h3(fn-1+4fn+fn+1),k=3yn+1-yn-2=3h8(fn-2+3fn-1+3fn+fn+1),k=4yn+2-yn-2=2h45(7fn-2+32fn-1+12fn+32fn+1+7fn+1),k=5yn+2-yn-3=5h288(19fn-3+75fn-2+50fn-1+50fn+75fn+1+19fn+2),k=6yn+3-yn-3=h140(41fn-3+216fn-2+27fn-1+272fn+27fn+1+216fn+2+41fn+3),k=7yn+3-yn-4=7h17280(751fn-4+3577fn-3+1323fn-2+2989fn-1+2989fn+1323fn+1+3577fn+2+751fn+3),k=8yn+4-yn-4=4h14175(989fn-4+5888fn-3-928fn-2+10496fn-1-4540fn+10496fn+1-928fn+2+5888fn+3+989fn+4).In the present paper, we will investigate the casek=8 and we will produce trigonometrically fitted differential methods of order 1.


## 3.2. Development of Closed Newton-Cotes Differential Schemes

For the development of a Newton-Cotes differential method of the above form, two procedures can be applied.(i)
The procedure which is based on the minimization of the local truncation error. Based on this procedure and for the casek=8, we can produce the well known coefficients:  a4=3956/14175, a3=23552/14175, a2=-3712/14175, a1=41984/14175, and a0=-3632/2835 (see the closed Newton-Cotes differential scheme for k=8 presented above).(ii)
The procedure which is based on
(1)
the minimization of the local truncation error;(2)
the satisfaction of the condition:2(∑i=1n|ai|+|a0|)/2n=1. This condition is produced by application of least squares method to the production of the differential method (see more details in [21] (about stable quadrature rules) and [22]).The above procedure leads to the following coefficients for the case:n=4: a0=-19672/945+70a4, a1=1952/105-56a4, a2=-848/105+28a4, a3=736/189-8a4 and to the condition 9836/33075≤a4≤244/735. We choose the value: a4=3/10, which satisfies the above condition.


## 3.3. Exponentially Fitted Closed Newton-Cotes Differential Method

Requiring the differential scheme:(3.5)yn+4-yn-4=h(a4fn-4+a3fn-3+a2fn-2+a1fn-1+a0fn+a1fn+1+a2fn+2+a3fn+3+a4fn+4)
to be accurate for the following set of functions (we note that fi=yi′,i=n-1,n,n+1): (3.6){1,x,x2,x3,x4,x5,x6,cos⁡(vx),sin⁡(vx)},
the following set of equations is obtained:(3.7)8cos⁡(w)sin⁡(w)(2(cos⁡(w))2-1)=w[(cos⁡(w))42a4-2a2+a0-6a3cos⁡(w)+4a2(cos⁡(w))2+2a1cos⁡(w)+8a3(cos⁡(w))3+16a4(cos⁡(w))4-16a4(cos⁡(w))2]2a4+2a3+2a2+2a1+a0=896a4+54a3+24a2+6a1=1282560a4+810a3+160a2+10a1=2048.
Requesting that a4=3/10 and solving the above system of equations, we obtain(3.8)a0=1735wcos⁡(w)+1353wcos⁡(3w)+242wcos⁡(2w)+270wcos⁡(4w)-900sin⁡(4w)675wcos⁡(w)+45wcos⁡(3w)-450w-270wcos⁡(2w),a1=-1856wcos⁡(3w)-1735w-1404wcos⁡(2w)-405wcos⁡(4w)+1350sin⁡(4w)1350wcos⁡(w)+90wcos⁡(3w)-900w-540wcos⁡(2w),a2=702wcos⁡(w)+418wcos⁡(3w)-121w+81wcos⁡(4w)-270sin⁡(4w)675wcos⁡(w)+45wcos⁡(3w)-450w-270wcos⁡(2w),a3=-27wcos⁡(4w)-836wcos⁡(2w)-1353w+1856wcos⁡(w)+90sin⁡(4w)1350wcos⁡(w)+90wcos⁡(3w)-900w-540wcos⁡(2w),
where w=vh.For small values ofv, the above formulae are subject to heavy cancellations. In this case the following Taylor series expansions must be used:(3.9)a0=46613780-40735670w2+40193249480w4-6980443681080400w6+945598949037788800w8-195517098336424096000w10+694578131900704693888000w12+33104793794390627842881280000w14+582548167731615751046180311040000w16+9723295174775617148961238556672000w18+⋯,a1=499525+40737560w2-40193332640w4+6980443908107200w6-945598965383718400w8+1955170911115232128000w10-694578132534272925184000w12-33104793795854170457175040000w14-582548167732154334728240414720000w16-97232951747100822865281651408896000w18+⋯,a2=7811050-407318900w2+40193831600w4-69804432270268000w6+9455989163459296000w8-1955170927788080320000w10+694578136335682312960000w12+331047937914635426142937600000w14+582548167735385836820601036800000w16+97232951747252057163204128522240000w18+⋯,a3=64934725+4073113400w2-401934989600w4+698044313621608000w6-9455989980755776000w8+19551709166728481920000w10-6945781338014093877760000w12-331047937987812556857625600000w14-5825481677332315020923606220800000w16-972329517471512342979224771133440000w18+⋯.The behaviour of the coefficients is given in the following Figure1.Behavior of the coefficients of the new proposed method given by (3.8) for several values of w.
(a)(b)(c)(d)The local truncation error for the above differential method is given by(3.10)L.T.E(h)=-593h928350(yn(9)+v2yn(7)).
The L.T.E is obtained expanding the terms yn±j and fn±j,j=1(1)4 in (3.5) into Taylor series expansions and substituting the Taylor series expansions of the coefficients of the method.In Figure2, we present the behaviour of the quantity ST=(2∑i=1n|ai|+|a0|)/2n for several values of v.Figure 2
Behaviour of the quantity ST for several values ofv.So, we have the following theorem.Theorem 3.2.
The method (3.5) with coefficients ai,i=0(1)4, obtained by the solution of the system (3.7) is accurate for the set of functions (3.6) and is of eighth algebraic order.


## 4. Closed Newton-Cotes Can Be Expressed as Symplectic Integrators

Let consider Hamilton’s equations of motion:(4.1)u̇=my,ẏ=-mu,
where m is a constant scalar or matrix. It is well known that (4.1) is important in the fields of physics, chemistry, material sciences, and so forth.Theorem 4.1.
A discrete scheme of the form:(4.2)(w-zzw)(un+1yn+1)=(wz-zw)(unyn)
is symplectic.Proof.
We rewrite (4.2) as
(4.3)(un+1yn+1)=(w-zzw)-1(wz-zw)(unyn).
Defining
(4.4)M=(w-zzw)-1(wz-zw)=1w2+z2(w2-z22wz-2wzw2-w2),
it can easily be proved that
(4.5)MTJM=J.
Thus, the matrix M is symplectic.
In [3], Zhu et al. have studied the well-known second-order differential scheme (SOD). They have proved that the scheme:
(4.6)qn+i-qn-i=2ihfn,i=1(1)4
has a symplectic structure.
The above methods have been produced by the simplest Open Newton-Cotes integral formula.
Based on [4, 7], the Closed Newton-Cotes differential schemes will be written as multilayer symplectic structures.
Application of the Newton-Cotes differential formula forn=4 to the linear Hamiltonian system (4.1) gives (4.7)un+4-un-4=s(a0yn-4+a1yn-3+a2yn-2+a3yn-1+a4yn+a5yn+1+a6yn+2+a7yn+3+a8yn+4),yn+4-yn-4=-s(a0un-4+a1un-3+a2un-2+a3un-1+a4un+a5un+1+a6un+2+a7un+3+a8un+4),
where s=mh, where m is defined in (4.1).
From (4.6), we have that(4.8)un+i-un-i=2isyn,yn+i-yn-i=-2isun,i=1(1)4ori=12(1)52.
We now consider the approximation based on the first formula of (4.8) for (n+1)-step gives (taking into account the second formula of (4.8))(4.9)un+i+un-i=(un+syn+i-1/2)+(un-syn-i+1/2)=un+i-1+un-i+1+s(yn+i-1/2-yn-i+1/2)=(2-i2s2)un,i=1(1)3.
Substituting (4.9) into (4.7) and considering that a0=a8, a1=a7, a2=a6, and a3=a5, we have: (4.10)un+4-un-4=s[a0(yn-4+yn+4)+(a1(2-9s2)+2a2(1-2s2)+a3(2-s2)+a4)yn]yn+4-yn-4=s[a0(un-4+un+4)+(a1(2-9s2)+2a2(1-2s2)+a3(2-s2)+a4)un],
and with (4.8) we have
(4.11)un+4-un-4=s[un+4-un-48sa0(yn-4+yn+4)+(a1(2-9s2)+2a2(1-2s2)+a3(2-s2)+a4)un+4-un-48s],yn+4-yn-4=s[yn+4-yn-48sa0(un-4+un+4)+(a1(2-9s2)+2a2(1-2s2)+a3(2-s2)+a4)[-yn+4-yn-48s]],
which gives:
(4.12)(un+4-un-4)[1-a1(2-9s2)+2a2(1-2s2)+a3(2-s2)+a48]=sa0(yn+4+yn-4)(yn+4-yn-4)[1-a1(2-9s2)+2a2(1-2s2)+a3(2-s2)+a48]=-sa0(un+4+un-4).
The above formula in matrix form can be written as(4.13)(Q(s)-sa0sa0Q(s))(un+4yn+4)=(Q(s)sa0-sa0Q(s))(un-4yn-4),
where
(4.14)Q(s)=1-a1(2-9s2)+2a2(1-2s2)+a3(2-s2)+a48,
which is a discrete scheme of the form (4.2) and hence it is symplectic.


## 5. Numerical Example

### 5.1. A Nonlinear Orbital Problem

Consider the nonlinear system of equations:(5.1)u′′+ω2u=2uv-sin⁡(2ωx)(u2+v2)3/2,u(0)=1,u′(0)=0,v′′+ω2v=u2-v2-cos⁡(2ωx)(u2+v2)3/2,v(0)=0,v′(0)=ω.The analytical solution of the problem is the following:(5.2)u(x)=cos⁡(ωx),v(x)=sin⁡(ωx).The system of (5.1) has been solved for 0≤x≤1000 and ω=10 using the methods(i)
The eighth-order multistep method developed by Quinlan and Tremaine [23] (which is indicated as Method I).(ii)
The tenth-order multistep method developed by Quinlan and Tremaine [23] (which is indicated as Method II).(iii)
The twelfth-order multistep method developed by Quinlan and Tremaine [23] (which is indicated as Method III).(iv)
The Newton-Cotes classical tenth-algebraic-order differential method (which is indicated as Method IV), (with the term classical we mean the closed Newton-Cotes differential method with constant coefficients).(v)
The Newton-Cotes eight-algebraic-order differential method with constant coefficient which corresponds to the New Developed Method VII (which is indicated as Method V).(vi)
The Newton-Cotes tenth-algebraic-order differential method developed in [8] (which is indicated as Method VI).(vii)
The stable Newton-Cotes eight-algebraic-order trigonometrically fitted differential method (which is indicated as Method VII).For this problem, we havew=10. The numerical results obtained for the seven methods mentioned above were compared with the analytical solution. Figure 3 shows the absolute errors  Errmax  defined by(5.3)Errmax⁡=|log10[max⁡(‖u(x)calculated-u(x)theoretical‖,‖v(x)calculated-v(x)theoretical‖)]|,x∈[0,1000],
for several values of the number of function evaluations (NFEs).Figure 3
Errmax for several values of the number of function evaluations (NFE) for the Methods I–VII for the nonlinear orbital problem. The nonexistrnce of a value of Errmax indicates that for these values Errmax is negative.


### 5.2. Duffing’s Equation

Consider the nonlinear initial value problem:(5.4)y′′=-y-y3+0.002cos⁡(1.01t),y(0)=0.20042672806,u′(0)=0.The analytical solution of the problem is the following:(5.5)y(t)=0.200179477536cos⁡(1.01t)+2.4694614310-4cos⁡(3.03t)+3.0401410-7cos⁡(5.05t)+3.7410-10cos⁡(7.07t).The above equation (5.4) has been solved for 0≤x≤1000 using the methods mentioned above.For this problem, we havew=1. The numerical results obtained for the seven methods mentioned above were compared with the analytical solution. Figure 4 shows the absolute errors Errmax defined by (5.6)Errmax⁡=|log10[max⁡(‖y(x)calculated-y(x)theoretical‖)]|,x∈[0,1000],
for several values of the number of function evaluations (NFEs).Figure 4
Errmax for several values of the number of function evaluations (NFEs) for the Methods I–IV for Duffing’s. The nonexistence of a value of Errmax indicates that for these values Errmax is negative.We note here that analogous results for both problems are obtained for interval of integration[0,10000] or [0,1000000].


## 5.1. A Nonlinear Orbital Problem

Consider the nonlinear system of equations:(5.1)u′′+ω2u=2uv-sin⁡(2ωx)(u2+v2)3/2,u(0)=1,u′(0)=0,v′′+ω2v=u2-v2-cos⁡(2ωx)(u2+v2)3/2,v(0)=0,v′(0)=ω.The analytical solution of the problem is the following:(5.2)u(x)=cos⁡(ωx),v(x)=sin⁡(ωx).The system of (5.1) has been solved for 0≤x≤1000 and ω=10 using the methods(i)
The eighth-order multistep method developed by Quinlan and Tremaine [23] (which is indicated as Method I).(ii)
The tenth-order multistep method developed by Quinlan and Tremaine [23] (which is indicated as Method II).(iii)
The twelfth-order multistep method developed by Quinlan and Tremaine [23] (which is indicated as Method III).(iv)
The Newton-Cotes classical tenth-algebraic-order differential method (which is indicated as Method IV), (with the term classical we mean the closed Newton-Cotes differential method with constant coefficients).(v)
The Newton-Cotes eight-algebraic-order differential method with constant coefficient which corresponds to the New Developed Method VII (which is indicated as Method V).(vi)
The Newton-Cotes tenth-algebraic-order differential method developed in [8] (which is indicated as Method VI).(vii)
The stable Newton-Cotes eight-algebraic-order trigonometrically fitted differential method (which is indicated as Method VII).For this problem, we havew=10. The numerical results obtained for the seven methods mentioned above were compared with the analytical solution. Figure 3 shows the absolute errors  Errmax  defined by(5.3)Errmax⁡=|log10[max⁡(‖u(x)calculated-u(x)theoretical‖,‖v(x)calculated-v(x)theoretical‖)]|,x∈[0,1000],
for several values of the number of function evaluations (NFEs).Figure 3
Errmax for several values of the number of function evaluations (NFE) for the Methods I–VII for the nonlinear orbital problem. The nonexistrnce of a value of Errmax indicates that for these values Errmax is negative.


## 5.2. Duffing’s Equation

Consider the nonlinear initial value problem:(5.4)y′′=-y-y3+0.002cos⁡(1.01t),y(0)=0.20042672806,u′(0)=0.The analytical solution of the problem is the following:(5.5)y(t)=0.200179477536cos⁡(1.01t)+2.4694614310-4cos⁡(3.03t)+3.0401410-7cos⁡(5.05t)+3.7410-10cos⁡(7.07t).The above equation (5.4) has been solved for 0≤x≤1000 using the methods mentioned above.For this problem, we havew=1. The numerical results obtained for the seven methods mentioned above were compared with the analytical solution. Figure 4 shows the absolute errors Errmax defined by (5.6)Errmax⁡=|log10[max⁡(‖y(x)calculated-y(x)theoretical‖)]|,x∈[0,1000],
for several values of the number of function evaluations (NFEs).Figure 4
Errmax for several values of the number of function evaluations (NFEs) for the Methods I–IV for Duffing’s. The nonexistence of a value of Errmax indicates that for these values Errmax is negative.We note here that analogous results for both problems are obtained for interval of integration[0,10000] or [0,1000000].


## 6. Conclusions

In this paper, we have introduced a new procedure for the development of Newton-Cotes differential schemes. The new procedure consists from the following steps:(i)
requirement the Newton-Cotes differential scheme to be accurate for the following set of functions:(6.1){1,x,x2,x3,…,xm,cos⁡(wx),sin⁡(wx)};(ii)
Satisfaction of the condition2(∑i=1n|ai|+|a0|)/2n=1, where ai,i=0(1)n are the coefficients of the Newton-Cotes differential scheme;(iii)
Expression of the Newton-Cotes differential scheme as multilayer symplectic integrators.We applied the new developed methods to several problems. We presented in this paper the application to a nonlinear orbital problem and to Duffing’s equation and we compared them with well-known integrators from the literature. Based on these illustrations, we conclude that the new procedure produces much more efficient methods than well-known methods of the literature.



---

*Source: 182536-2012-05-08.xml*