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

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.

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 Section 2.(ii)In Section 3, 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 Section 4, the conversion of the closed Newton-Cotes differential methods into multilayer symplectic structures is presented.(iv)Numerical results are presented in Section 5.

2. Basic Theory on Symplectic Schemes and Numerical Methods

Based on Zhu et al. [3] and on the division of the interval [𝑎,𝑏] with 𝑁 points, we have the following discrete scheme for the n-step approximation to the solution:𝑝_𝑛+1𝑞_𝑛+1=𝑀_𝑛+1𝑝_𝑛𝑞_𝑛,𝑀_𝑛+1=𝑤_𝑛+1𝑦_𝑛+1𝑧_𝑛+1𝑔_𝑛+1. (2.1) Based on the above we can write the n-step approximation to the solution as𝑝_𝑛𝑞_𝑛=𝑤_𝑛𝑦_𝑛𝑧_𝑛𝑔_𝑛𝑤_𝑛−1𝑦_𝑛−1𝑧_𝑛−1𝑔_𝑛−1⋯𝑤₁𝑦₁𝑧₁𝑔₁𝑝₀𝑞₀=𝑀_𝑛𝑀_𝑛−1⋯𝑀₁𝑝₀𝑞₀. (2.2) Defining𝑆=𝑀_𝑛𝑀_𝑛−1⋯𝑀₁=𝑊_𝑛𝑌_𝑛𝑍_𝑛𝐺_𝑛, (2.3) the discrete transformation can be written as𝑝_𝑛𝑞_𝑛𝑝=𝑆⁰𝑞₀. (2.4)

A discrete scheme (2.1) is a symplectic scheme if the transformation matrix 𝑆 is symplectic. A matrix 𝐴 is symplectic if 𝐴^𝑇𝐽𝐴=𝐽, where𝐽=01−10. (2.5) The product of symplectic matrices is also symplectic. Hence, if each matrix 𝑀_𝑛 is symplectic the transformation matrix 𝑆 is symplectic. Consequently, the discrete scheme (3.5) is symplectic if each matrix 𝑀_𝑛 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:^𝑏_𝑎𝑓(𝑥)𝑑𝑥≈𝑧ℎ^𝑘_𝑖=0𝑡_𝑖𝑓𝑥_𝑖, (3.1) whereℎ=𝑏−𝑎𝑁,𝑥_𝑖=𝑎+𝑖ℎ,𝑖=0,1,2,…,𝑁. (3.2) The coefficient 𝑧 as well as the weights 𝑡_𝑖 are given in Table 1.

Table 1 

Closed Newton-Cotes integral rules.

Remark 3.1. It is easy for one to see that the coefficients 𝑡_𝑖 in the Table 1 are symmetric, that is, one has the following relation: 𝑡_𝑖=𝑡_𝑘−𝑖𝑘,𝑖=0,1,…,2. (3.3)

The closed Newton-Cotes differential methods are produced from the integral rules. From Table 1 we have the following differential methods: 𝑘=1𝑦_𝑛+1−𝑦_𝑛=ℎ2𝑓_𝑛+1+𝑓_𝑛,𝑘=2𝑦_𝑛+1−𝑦_𝑛−1=ℎ3𝑓_𝑛−1+4𝑓_𝑛+𝑓_𝑛+1,𝑘=3𝑦_𝑛+1−𝑦_𝑛−2=3ℎ8𝑓_𝑛−2+3𝑓_𝑛−1+3𝑓_𝑛+𝑓_𝑛+1,𝑘=4𝑦_𝑛+2−𝑦_𝑛−2=2ℎ457𝑓_𝑛−2+32𝑓_𝑛−1+12𝑓_𝑛+32𝑓_𝑛+1+7𝑓_𝑛+1,𝑘=5𝑦_𝑛+2−𝑦_𝑛−3=5ℎ28819𝑓_𝑛−3+75𝑓_𝑛−2+50𝑓_𝑛−1+50𝑓_𝑛+75𝑓_𝑛+1+19𝑓_𝑛+2,𝑘=6𝑦_𝑛+3−𝑦_𝑛−3=ℎ14041𝑓_𝑛−3+216𝑓_𝑛−2+27𝑓_𝑛−1+272𝑓_𝑛+27𝑓_𝑛+1+216𝑓_𝑛+2+41𝑓_𝑛+3,𝑘=7𝑦_𝑛+3−𝑦_𝑛−4=7ℎ17280751𝑓_𝑛−4+3577𝑓_𝑛−3+1323𝑓_𝑛−2+2989𝑓_𝑛−1+2989𝑓_𝑛+1323𝑓_𝑛+1+3577𝑓_𝑛+2+751𝑓_𝑛+3,𝑘=8𝑦_𝑛+4−𝑦_𝑛−4=4ℎ14175989𝑓_𝑛−4+5888𝑓_𝑛−3−928𝑓_𝑛−2+10496𝑓_𝑛−1−4540𝑓_𝑛+10496𝑓_𝑛+1−928𝑓_𝑛+2+5888𝑓_𝑛+3+989𝑓_𝑛+4. (3.4)

In the present paper, we will investigate the case 𝑘=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 case 𝑘=8, we can produce the well known coefficients: 𝑎₄=3956/14175, 𝑎₃=23552/14175, 𝑎₂=−3712/14175, 𝑎₁=41984/14175, and 𝑎₀=−3632/2835 (see the closed Newton-Cotes differential scheme for 𝑘=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(∑^𝑛_𝑖=1|𝑎_𝑖|+|𝑎₀|)/2𝑛=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: 𝑛=4: 𝑎₀=−19672/945+70𝑎₄, 𝑎₁=1952/105−56𝑎₄, 𝑎₂=−848/105+28𝑎₄, 𝑎₃=736/189−8𝑎₄ and to the condition 9836/33075≤𝑎₄≤244/735. We choose the value: 𝑎₄=3/10, which satisfies the above condition.

3.3. Exponentially Fitted Closed Newton-Cotes Differential Method

Requiring the differential scheme: 𝑦_𝑛+4−𝑦_𝑛−4=ℎ𝑎₄𝑓_𝑛−4+𝑎₃𝑓_𝑛−3+𝑎₂𝑓_𝑛−2+𝑎₁𝑓_𝑛−1+𝑎₀𝑓_𝑛+𝑎₁𝑓_𝑛+1+𝑎₂𝑓_𝑛+2+𝑎₃𝑓_𝑛+3+𝑎₄𝑓_𝑛+4 (3.5) to be accurate for the following set of functions (we note that 𝑓_𝑖=𝑦^′_𝑖,𝑖=𝑛−1,𝑛,𝑛+1): 1,𝑥,𝑥²,𝑥³,𝑥⁴,𝑥⁵,𝑥⁶,cos(v𝑥),sin(v𝑥), (3.6) the following set of equations is obtained:8cos(𝑤)sin(𝑤)2(cos(𝑤))²−1=𝑤2𝑎₄−2𝑎₂+𝑎₀−6𝑎₃cos(𝑤)+4𝑎₂(cos(𝑤))²+2𝑎₁cos(𝑤)+8𝑎₃(cos(𝑤))³+16𝑎₄(cos(𝑤))⁴−16𝑎₄(cos(𝑤))²2𝑎₄+2𝑎₃+2𝑎₂+2𝑎₁+𝑎₀=896𝑎₄+54𝑎₃+24𝑎₂+6𝑎₁=1282560𝑎₄+810𝑎₃+160𝑎₂+10𝑎₁=2048. (3.7) Requesting that 𝑎₄=3/10 and solving the above system of equations, we obtain𝑎₀=1735𝑤cos(𝑤)+1353𝑤cos(3𝑤)+242𝑤cos(2𝑤)+270𝑤cos(4𝑤)−900sin(4𝑤)(𝑤)()(),𝑎675𝑤cos+45𝑤cos3𝑤−450𝑤−270𝑤cos2𝑤₁=−1856𝑤cos(3𝑤)−1735𝑤−1404𝑤cos(2𝑤)−405𝑤cos(4𝑤)+1350sin(4𝑤)(𝑤)()(),𝑎1350𝑤cos+90𝑤cos3𝑤−900𝑤−540𝑤cos2𝑤₂=702𝑤cos(𝑤)+418𝑤cos(3𝑤)−121𝑤+81𝑤cos(4𝑤)−270sin(4𝑤)(𝑤)()(),𝑎675𝑤cos+45𝑤cos3𝑤−450𝑤−270𝑤cos2𝑤₃=−27𝑤cos(4𝑤)−836𝑤cos(2𝑤)−1353𝑤+1856𝑤cos(𝑤)+90sin(4𝑤)(𝑤)()(),1350𝑤cos+90𝑤cos3𝑤−900𝑤−540𝑤cos2𝑤 (3.8) where 𝑤=vℎ.

For small values of v, the above formulae are subject to heavy cancellations. In this case the following Taylor series expansions must be used:𝑎₀=4661−37804073𝑤5670²+40193𝑤249480⁴−6980443𝑤681080400⁶+9455989𝑤49037788800⁸−19551709𝑤8336424096000¹⁰[phone]𝑤1900704693888000¹²[phone]𝑤4390627842881280000¹⁴[phone]𝑤1615751046180311040000¹⁶[phone]𝑤75617148961238556672000¹⁸𝑎+⋯,₁=499+5254073𝑤7560²−40193𝑤332640⁴+6980443𝑤908107200⁶−9455989𝑤65383718400⁸[phone]𝑤11115232128000¹⁰−69457813𝑤2534272925184000¹²−3310479379𝑤5854170457175040000¹⁴−58254816773𝑤2154334728240414720000¹⁶−97232951747𝑤100822865281651408896000¹⁸𝑎+⋯,₂=781−10504073𝑤18900²+40193𝑤831600⁴−6980443𝑤2270268000⁶+9455989𝑤163459296000⁸−19551709𝑤27788080320000¹⁰[phone]𝑤6335682312960000¹²[phone]𝑤14635426142937600000¹⁴[phone]𝑤5385836820601036800000¹⁶[phone]𝑤252057163204128522240000¹⁸𝑎+⋯,₃=6493+47254073𝑤113400²−40193𝑤4989600⁴+6980443𝑤13621608000⁶−9455989𝑤980755776000⁸[phone]𝑤166728481920000¹⁰−69457813𝑤38014093877760000¹²−3310479379𝑤87812556857625600000¹⁴−58254816773𝑤32315020923606220800000¹⁶−97232951747𝑤1512342979224771133440000¹⁸+⋯. (3.9)

The behaviour of the coefficients is given in the following Figure 1.

(a)

(b)

(c)

(d)

Figure 1 

Behavior of the coefficients of the new proposed method given by (3.8) for several values of 𝑤.

The local truncation error for the above differential method is given byL.T.E(ℎ)=−593ℎ⁹𝑦28350^𝑛(9)+v²𝑦_𝑛⁽⁷⁾. (3.10) The L.T.E is obtained expanding the terms 𝑦_𝑛±𝑗 and 𝑓_𝑛±𝑗,𝑗=1(1)4 in (3.5) into Taylor series expansions and substituting the Taylor series expansions of the coefficients of the method.

In Figure 2, we present the behaviour of the quantity ST=(2∑^𝑛_𝑖=1|𝑎_𝑖|+|𝑎₀|)/2𝑛 for several values of v.

Figure 2 

Behaviour of the quantity ST for several values of v.

So, we have the following theorem.

Theorem 3.2. The method (3.5) with coefficients 𝑎_𝑖,𝑖=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) where 𝑚 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: 𝑤−𝑧𝑢𝑧𝑤𝑦^𝑛+1_𝑛+1=𝑤𝑧𝑢−𝑧𝑤𝑦^𝑛_𝑛 (4.2) is symplectic.

Proof. We rewrite (4.2) as 𝑢_𝑛+1𝑦_𝑛+1=𝑤−𝑧𝑧𝑤⁻¹𝑤𝑧𝑢−𝑧𝑤𝑦^𝑛_𝑛. (4.3) Defining 𝑀=𝑤−𝑧𝑧𝑤⁻¹=𝑤𝑧1−𝑧𝑤𝑤²+𝑧²𝑤²−𝑧²2𝑤𝑧−2𝑤𝑧𝑤²−𝑤², (4.4) it can easily be proved that 𝑀^𝑇𝐽𝑀=𝐽. (4.5) Thus, the matrix 𝑀 is symplectic.
In [3], Zhu et al. have studied the well-known second-order differential scheme (SOD). They have proved that the scheme: 𝑞_𝑛+𝑖−𝑞_𝑛−𝑖=2𝑖ℎ𝑓_𝑛,𝑖=1(1)4 (4.6) 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 for 𝑛=4 to the linear Hamiltonian system (4.1) gives 𝑢_𝑛+4−𝑢_𝑛−4=𝑠𝑎₀𝑦_𝑛−4+𝑎₁𝑦_𝑛−3+𝑎₂𝑦_𝑛−2+𝑎₃𝑦_𝑛−1+𝑎₄𝑦_𝑛+𝑎₅𝑦_𝑛+1+𝑎₆𝑦_𝑛+2+𝑎₇𝑦_𝑛+3+𝑎₈𝑦_𝑛+4,𝑦_𝑛+4−𝑦_𝑛−4=−𝑠𝑎₀𝑢_𝑛−4+𝑎₁𝑢_𝑛−3+𝑎₂𝑢_𝑛−2+𝑎₃𝑢_𝑛−1+𝑎₄𝑢_𝑛+𝑎₅𝑢_𝑛+1+𝑎₆𝑢_𝑛+2+𝑎₇𝑢_𝑛+3+𝑎₈𝑢_𝑛+4, (4.7) where 𝑠=𝑚ℎ, where 𝑚 is defined in (4.1).
From (4.6), we have that𝑢_𝑛+𝑖−𝑢_𝑛−𝑖=2𝑖𝑠𝑦_𝑛,𝑦_𝑛+𝑖−𝑦_𝑛−𝑖=−2𝑖𝑠𝑢_𝑛(1)1,𝑖=14or𝑖=2(1)52. (4.8)
We now consider the approximation based on the first formula of (4.8) for (𝑛+1)-step gives (taking into account the second formula of (4.8))𝑢_𝑛+𝑖+𝑢_𝑛−𝑖=𝑢_𝑛+𝑠𝑦_{𝑛+𝑖−1/2}+𝑢_𝑛−𝑠𝑦_{𝑛−𝑖+1/2}=𝑢_{𝑛+𝑖−1}+𝑢_{𝑛−𝑖+1}+𝑠𝑦_{𝑛+𝑖−1/2}−𝑦_{𝑛−𝑖+1/2}=2−𝑖²𝑠²𝑢_𝑛,𝑖=1(1)3. (4.9)
Substituting (4.9) into (4.7) and considering that 𝑎₀=𝑎₈, 𝑎₁=𝑎₇, 𝑎₂=𝑎₆, and 𝑎₃=𝑎₅, we have: 𝑢_𝑛+4−𝑢_𝑛−4=𝑠𝑎₀𝑦_𝑛−4+𝑦_𝑛+4+𝑎₁2−9𝑠²+2𝑎₂1−2𝑠²+𝑎₃2−𝑠²+𝑎₄𝑦_𝑛𝑦_𝑛+4−𝑦_𝑛−4=𝑠𝑎₀𝑢_𝑛−4+𝑢_𝑛+4+𝑎₁2−9𝑠²+2𝑎₂1−2𝑠²+𝑎₃2−𝑠²+𝑎₄𝑢_𝑛, (4.10) and with (4.8) we have 𝑢_𝑛+4−𝑢_𝑛−4=𝑠𝑎₀𝑦_𝑛−4+𝑦_𝑛+4+𝑎₁2−9𝑠²+2𝑎₂1−2𝑠²+𝑎₃2−𝑠²+𝑎₄𝑢_𝑛+4−𝑢_𝑛−4,𝑦8𝑠_𝑛+4−𝑦_𝑛−4=𝑠𝑎₀𝑢_𝑛−4+𝑢_𝑛+4+𝑎₁2−9𝑠²+2𝑎₂1−2𝑠²+𝑎₃2−𝑠²+𝑎₄−𝑦_𝑛+4−𝑦_𝑛−4,8𝑠 (4.11) which gives: 𝑢_𝑛+4−𝑢_𝑛−4𝑎1−¹2−9𝑠²+2𝑎₂1−2𝑠²+𝑎₃2−𝑠²+𝑎₄8=𝑠𝑎₀𝑦_𝑛+4+𝑦_𝑛−4𝑦_𝑛+4−𝑦_𝑛−4𝑎1−¹2−9𝑠²+2𝑎₂1−2𝑠²+𝑎₃2−𝑠²+𝑎₄8=−𝑠𝑎₀𝑢_𝑛+4+𝑢_𝑛−4. (4.12)
The above formula in matrix form can be written as𝑄(𝑠)−𝑠𝑎₀𝑠𝑎₀𝑄(𝑠)𝑢_𝑛+4𝑦_𝑛+4=𝑄(𝑠)𝑠𝑎₀−𝑠𝑎₀𝑄(𝑠)𝑢_𝑛−4𝑦_𝑛−4, (4.13) where 𝑄(𝑠)𝑎=1−¹2−9𝑠²+2𝑎₂1−2𝑠²+𝑎₃2−𝑠²+𝑎₄8, (4.14) 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: 𝑢^′′+𝜔²()𝑢=2𝑢𝑣−sin2𝜔𝑥𝑢²+𝑣^23/2,𝑢(0)=1,𝑢^′(0)𝑣=0,_′′+𝜔²𝑢𝑣=²−𝑣²−cos(2𝜔𝑥)𝑢²+𝑣^23/2,𝑣(0)=0,𝑣^′(0)=𝜔. (5.1)

The analytical solution of the problem is the following: 𝑢(𝑥)=cos(𝜔𝑥),𝑣(𝑥)=sin(𝜔𝑥). (5.2)

The system of (5.1) has been solved for 0≤𝑥≤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 have 𝑤=10. The numerical results obtained for the seven methods mentioned above were compared with the analytical solution. Figure 3 shows the absolute errors Err_max defined byErr_max=log₁₀max𝑢(𝑥)_calculated−𝑢(𝑥)_theoretical,𝑣(𝑥)_calculated−𝑣(𝑥)_theoretical,𝑥∈[0,1000], (5.3) 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: 𝑦^′′=−𝑦−𝑦³+0.002cos(1.01𝑡),𝑦(0)=0.20042672806,𝑢^′(0)=0. (5.4)

The analytical solution of the problem is the following: 𝑦(𝑡)=0.200179477536cos(1.01𝑡)+2.4694614310⁻⁴cos(3.03𝑡)+3.0401410⁻⁷cos(5.05𝑡)+3.7410⁻¹⁰cos(7.07𝑡). (5.5)

The above equation (5.4) has been solved for 0≤𝑥≤1000 using the methods mentioned above.

For this problem, we have 𝑤=1. The numerical results obtained for the seven methods mentioned above were compared with the analytical solution. Figure 4 shows the absolute errors Err_max defined by Err_max=log₁₀max𝑦(𝑥)_calculated−𝑦(𝑥)_theoretical,𝑥∈[0,1000], (5.6) 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: 1,𝑥,𝑥²,𝑥³,…,𝑥^𝑚,cos(𝑤𝑥),sin(𝑤𝑥); (6.1)(ii)Satisfaction of the condition 2(∑^𝑛_𝑖=1|𝑎_𝑖|+|𝑎₀|)/2𝑛=1, where 𝑎_𝑖,𝑖=0(1)𝑛 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.