New Eight-Step Symmetric Embedded [603941]

New Eight-Step Symmetric Embedded
Predictor-Corrector (EPC2M) Method with
Vanished Phase-lag and its First and Second
Derivatives for the Numerical Integration of
the Schr¨ odinger equation
P.I. Stasinos
Laboratory of Computational Sciences, Department of Infor matics and
Telecommunications, University of Peloponnese, GR-221 00 Tripolis, Greece
T.E. Simos1
Department of Mathematics, College of Sciences, King Saud Un iversity, P.
O. Box 2455, Riyadh 11451, Saudi Arabia and Laboratory of Com putational
Sciences, Department of Informatics and Telecommunicatio ns, University of
Peloponnese, GR-221 00 Tripolis, Greece. Tel: 0030 210 94 21 510, Fax:
0030 210 94 20 091, e-mail: [anonimizat]
Abstract
Weconstruct a new embedded predictor-corrector method, pha se-fitted with
infinite order of phase-lag . It is an eight-step symmetric multistep m ethod
related to the multistep method of Quinlan and Tremaine. We introduc e
a new pair of embedded predictor-corrector methods that can be used to
solve numerically IVPs with oscillatory solutions , the Schr¨ odinger eq uation
and orbital problems. It is about a new method whith vanished phase -lag,
tenth algebraic order and it has vanished its first and second deriva tives.
After testing the new method , according to the numerical results , this new
method is more efficient than other methods that were compared.
1Highly Cited Researcher ( http://isihighlycited.com/ ),Active Member of the Eu-
ropean Academy of Sciences and Arts. Active Member of the Europ ean Academy of
Sciences. Corresponding Member of European Academy of Arts, S ciences and Humani-
ties, Please use the following address for all correspondence: Dr. T.E. Simos, 10 Konitsis
Street, Amfithea – Paleon Faliron, GR-175 64 Athens, Greece.
Preprint submitted to Computer Physics Communications Sep tember 20, 2018

Keywords: multistep , IVPs, predictor-correcto , symmetric, phase-fitted ,
phase-lag, embedded, Schr¨ odinger equation, P-stable, interva l of periodicity
2000 MSC: 65L05
PACS:0.00, 0.00.E
1. Introduction
We are dealing with the numerical solution of the second-order IVPs :
q′′=f(x,y),q(x0) =η,q′(x0) =η′(1.0.1)
for which the first derivative q′is missing. Dahlquist [ ?] suggests for the
equation of the form (1.0.1) the following method:
λ/summationdisplay
i=0aiyn+i=h2λ′/summationdisplay
i=0bifn+i (1.0.2)
This method is related to the characteristic polynomials [ ?] :
ρ(w) =λ/summationdisplay
i=0aiwi(1.0.3)
and
σ(w) =λ′/summationdisplay
i=0biwi. (1.0.4)
sowecanrefertothismethodas( ρ,σ)method. Themethodiscalledimplicit
ifλ=λ′and ifλ′=λ−1 explicit. Regarding the accuracy , the stability
and the convergence of linear multistep methods Dahlquist [ ?] presented
a complete work. The method has order q, if for any test function z(x) we
have:
λ/summationdisplay
i=0aiz(x+ih)−h2λ′/summationdisplay
i=0biz′′(x+ih) =Cq+2hq+2z(q+2)(x)+O(hq+3) (1.0.5)
The error constant Cq+2is given by :
Cp=1
p!λ/summationdisplay
j=0jp−2(j2aj−p(p−1)bj)−λ′/summationdisplay
λ+1jp−2
(p−2)!bj,p>2 (1.0.6)
We know that this method satisfies [ ?] the following :
2

• |a0|+|b0| /negationslash= 0,aλ= 1
•ρandσshould not have any common factors
•ρ(1) = 0,ρ′(1) = 0,ρ′′(1) = 2σ(1)
•is zero-stable method
Also , if :
ai=aλ−i,i= 0(1)λ,bi=bλ−i,i= 0(1)λ (1.0.7)
then the method is called symmetric. The following two definitions link th e
concepts P-stability and interval of periodicity , concerning form p roblems :
y′′+q2y= 0. (1.0.8)
DEFINITION 1. [ ?] The (ρ,σ) method , is said to have an interval of
periodicity (0 ,∞) , if the roots of :
Π(w,H2) =ρ(w)+H2σ(w) = 0,H=qh (1.0.9)
satisfy
w1=eiθ(H),w2=e−iθ(H),|wk| ≤1,k>2 (1.0.10)
for allH2in the interval and θ(H) is a real function.
DEFINITION 2. [ ?] If (ρ,σ) method has an interval of periodicity of the
form (0,∞) , then the method is said to be P-stable.
Lambert and Watson [ ?] proved that if ( ρ,σ) method is symmetric then
it has a non-vanishing interval of periodicity and for P-stability, the order
is at most two. The following theorem gives the result of Fukushima [ ?],
where given the necessary condition for which a symmetric multistep method
has a non-vanishing interval of periodicity.
Theorem 1. [?] Consider a symmetric multistep method for which we
define a function :
d(θ) =−ρ(eiθ)
σ(eiθ)(1.0.11)
The method may have a non-vanishing interval of periodicity if and on ly if :
•d(θ) in the interval [0 ,π] has no nonzero double roots , or
3

•all the nonzero roots of d(θ) , of multiplicity two , satisfy the condition
d′′(θ)>0 in the interval [0 ,π].
Quinlan and Tremaine [ ?] continued the work of Lambert and Watson to
achieve high-order symmetrical methods for planetary integratio ns. They in-
troducedtheirresultsontheincredibleperformanceoftheirmeth odsforlong-
timeintegrations. ThecontributionofQuinlanandTremaineonthepr operty
of symmetric multistep methods where the longtitude error increas es linearly
over time , is also significant.
The association of multistep symmetric methods with St¨ ormer symm etric
methods was done by Quinlan and Tremaine [ ?]. They present a method
of twelfth order, symmetric and more efficient than the symmetric S t¨ ormer
method of thirteenth order. With regard to the bibliography for th e above
category of problems (1.0.1) you can refer to [ ?]-[?],[?]-[?].
WeareparticularlydealingwiththeSchr¨ odinger equation.TheSchr¨ odinger
equation is:
−/planckover2pi12
2m∂2ψ(x)
∂x2+V(x)ψ(x) =Eψ(x) (1.0.12)
whereψ(x) is a wavefunction that represents the energy self-condition, alo ng
with the normalization condition:
ˆH=−/planckover2pi12
2m∂2
∂x2+V(ˆx) (1.0.13)
If we define as υ(˜x) = 2mL2V(x)//planckover2pi12= 2mL2V(˜xL)//planckover2pi12and as
ǫ= 2mL2E//planckover2pi12and change symbols as ˜ x→x,˜ψ→ψ, we will have :
ψ′′(x) =−(ǫ−υ(x))ψ(x) (1.0.14)
We observe that the equation (1.0.12) is transformed into equation (1.0.14)
which refers to the initial value problems (1.0.1).
2. Phase-lag of symmetric multistep methods
Regarding problems :
q′′=f(x,y),q(x0) =η,q′(x0) =η′(2.0.1)
we apply an λ−stepmethod of the form :
λ/summationdisplay
i=0aizn+i=h2λ/summationdisplay
i=0bif(xn+i,zn+i) (2.0.2)
4

withh=|xi+1−xi|,i= 0(1)m−1, and|a0|+|b0| /negationslash= 0.
Applying the scalar test equation :
z′′=−t2z (2.0.3)
in a 2λ-step symmetric method, for i=−λ(1)λ, resulting a difference equa-
tion :
λ/summationdisplay
i=1Ai(v)(zn+i+zn−i)+A0(v)zn= 0 (2.0.4)
whereA0(v),A1(v,…,Aλ(v) are polynomials and v=th,his the step
length.
Then from the form (2.0.4) we get the characteristic equation :
λ/summationdisplay
i=1Ai(v)(wi+w−i)+A0(v) = 0 (2.0.5)
which relates to DEFINITION 1.
T.E.Simos and P.S.Williams [ ?] in 1997 presented a direct formula of
calculating the phase-lag for a symmetric four-step method. This f ormula
helpedtoconstruct afour-stepmethodwithinfinite orderofphas e-lag. Their
theorem is :
Theorem 2. [?]A symmetric 2λ-step method , with corresponding char-
acteristic equation (2.0.5) has phase-lag constant θand phase-lag order p
determined by the formula :
−θvp+2+O(vp+4) =2λ/summationtext
j=1Aj(v)cos(jv)+A0(v)
2λ/summationtext
j=1j2Aj(v)(2.0.6)
So for the new method we have :
−θvp+2+O(vp+4) =2A4(v)cos(4v)+2A3(v)cos(3v)+2A2(v)cos(2v)+2A1(v)cos(v)+A0(v)
32A4(v)+18A3(v)+8A2(v)+2A1(v)
(2.0.7)
as it is an eight-step symmetric method.
5

3. New Embedded Predictor-Corrector
In this section we present the theoretical background on which we relied
for the construction of the new method. From Lambert [ ?] we have the
mainλ−steppredictor-corrector pair :
λ/summationtext
i=0a∗
izn+i=hλ−1/summationtext
i=0b∗
ifn+i
λ/summationtext
i=0aizn+i=hλ/summationtext
i=0bifn+i

(3.0.1)
If we denote q∗as the predictor order of accuracy and qas the corrector
order accuracy, then the predictor-corrector method has tot al accuracy q∗+r
ifq∗<qandr<=q−q∗−1 .
For the purposes of our research we consider the pair :
λ/summationtext
i=0a∗
izn+i=h2λ−1/summationtext
i=0b∗
ifn+i
λ/summationtext
i=0aizn+i=h2λ/summationtext
i=0bifn+i

(3.0.2)
of multistep linear methods where |a∗
0|+|b∗
0| /negationslash= 0,|a0|+|b0| /negationslash= 0,b∗
λ= 0
andbλ/negationslash= 0.
If we assume a∗
λ= 1 andaλ= 1 then we can take:
zn+λ+λ−1/summationtext
i=0a∗
izn+i=h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)
zn+λ+λ−1/summationtext
i=0aizn+i=h2/parenleftBig
bλf(xn+λ,zn+λ)+λ−1/summationtext
i=0bif(xn+i,zn+i)/parenrightBig

(3.0.3)
and so we can have :
zn+λ=−λ−1/summationtext
i=0a∗
izn+i+h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)
zn+λ=−λ−1/summationtext
i=0aizn+i+h2bλf(xn+λ,zn+λ)+h2λ−1/summationtext
i=0bif(xn+i,zn+i)

(3.0.4)
6

Thus the numerical integration for second-order initial-value prob lems takes
the form :
z∗
n+λ=−λ−1/summationtext
i=0a∗
izn+i+h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)
zn+λ=−λ−1/summationtext
i=0a∗
izn+i+h2bλf(xn+λ,z∗
n+λ)+h2λ−1/summationtext
i=0bif(xn+i,zn+i)

(3.0.5)
where|a∗
0|+|b∗
0| /negationslash= 0 ,|a0|+|b0| /negationslash= 0 andbλ/negationslash= 0.
Ifa∗
i=a∗
λ−i,b∗
i=b∗
λ−i,ai=aλ−iandbi=bλ−i,i= 0(1)⌊λ
2⌋we said that
the method is symmetric.
Now if we put λ= 8 into (3.0.5) we get :
z∗
4=−(z−4+a∗
3(z3+z−3)+a∗
2(z2+z−2)+a∗
1(z1+z−1)+a∗
0z0)
+h2(b∗
3(f3+f−3)+b∗
2(f2+f−2)+b∗
1(f1+f−1)+b∗
0f0)
z4=−(z−4+a3(z3+z−3)+a2(z2+z−2)+a1(z1+z−1)+a0z0)
+h2(b4(f4+f−4)+b3(f3+f−3)+b2(f2+f−2)+
+b1(f1+f−1)+b0f0)


(3.0.6)
which is the form of the predictor-corrector eight-step symmetr ic method ,
where
zi=z(x+ih),fi=f(x+ih,z(x+ih)),i=−4(1)3 andf4=f(x+4h,z∗
4).
The above method corresponds to the characteristic equation :
4/summationdisplay
i=1Ai(v)(wi+w−i)+A0(v) = 0 (3.0.7)
whereAi(v) =ai+(bi−a∗
ib4)v2−b∗
ib4v4, i= 0(1)4, a4=a∗
4= 1,b∗
4= 0.
3.1. Embedded predictor-corrector
Replacing as ai=a∗
i, i= 0(1)λ−1 to (3.0.5) we get :
z∗
n+λ=−λ−1/summationtext
i=0a∗
izn+i+h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)
zn+λ=−λ−1/summationtext
i=0a∗
izn+i+h2bλf(xn+λ,z∗
n+λ)+h2λ−1/summationtext
i=0bif(xn+i,zn+i)

(3.1.1)
7

where|a∗
0|+|b∗
0| /negationslash= 0 ,|a0|+|b0| /negationslash= 0 andbλ/negationslash= 0.
If we put :
bi=bi+0 =bi−b∗
i+b∗
i= (bi−b∗
i)+b∗
i,
in the (3.1.1) and βi=bi−b∗
i,i= 0(1)λ−1, then:
bi=βi+b∗
i
and so we get :
h2λ−1/summationdisplay
i=0bif(xn+i,zn+i) =h2λ−1/summationdisplay
i=0(βi+b∗
i)f(xn+i,zn+i) =
h2λ−1/summationdisplay
i=0βif(xn+i,zn+i)+h2λ−1/summationdisplay
i=0b∗
if(xn+i,zn+i) (3.1.2)
so it is written as:
z∗
n+λ=−λ−1/summationtext
i=0a∗
izn+i+h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)
zn+λ=−λ−1/summationtext
i=0a∗
izn+i+h2bλf(xn+λ,z∗
n+λ)
+h2λ−1/summationtext
i=0βif(xn+i,zn+i)+h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)

(3.1.3)
or
z∗
n+λ=−λ−1/summationtext
i=0a∗
izn+i+h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)
qzn+λ=−λ−1/summationtext
i=0a∗
izn+i+h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)
+h2bλf(xn+λ,z∗
n+λ)+h2λ−1/summationtext
i=0βif(xn+i,zn+i)

(3.1.4)
where|a∗
0|+|b∗
0| /negationslash= 0 ,|a∗
0|+|β0| /negationslash= 0 andbλ/negationslash= 0.
8

Because it is b∗
λ= 0 andbλ/negationslash= 0 and as βλ=bλ−b∗
λ, then we get:
βλ=bλ−b∗
λ=bλ−0 =bλ/negationslash= 0
and as:
βi=bi−b∗
i,i= 0(1)λ (3.1.5)
the (3.1.4) can take the new pair form :
z∗
n+λ=−λ−1/summationtext
i=0a∗
izn+i+h2λ−1/summationtext
i=0b∗
if(xn+i,zn+i)
zn+λ=z∗
n+λ+h2λ/summationtext
i=0βif(xn+i,zn+i)

(3.1.6)
where|a∗
0|+|b∗
0| /negationslash= 0,|a0|+|β0| /negationslash= 0. One can see that on this new pair of
methods ,the predictor is totally enclosed in corrector , so we call it a s the
Embedded Predictor-Corrector Method (EPCM) .
4. The new embedded (EPC2M) method
Foraλ= 1 to (2.0.2) we take :
zn+λ+λ−1/summationdisplay
i=0aizn+i=h2λ/summationdisplay
i=0bif(xn+i,zn+i),
and :
zn+λ=−λ−1/summationdisplay
i=0aizn+i+h2λ/summationdisplay
i=0bif(xn+i,zn+i) (4.0.1)
and ifai=aλ−iandbi=bλ−i,i= 0(1)⌊λ
2⌋we have a symmetric method.
4.1. The predictor (explicit method)
For an eight-step symmetric method , if we put b8= 0 to the form (4.0.1)
then we get the explicit method :
z4=−(z−4+a3(z3+z−3)+a2(z2+z−2)+a1(z1+z−1)+a0z0)
+h2(b3(f3+f−3)+b2(f2+f−2)+b1(f1+f−1)+b0f0)(4.1.1)
9

wherezi=z(x+ih) andfi=f(x+ih,z(x+ih)).
The characteristic equation (2.0.5) now becomes :
4/summationdisplay
i=1Ai(v)(wi+w−i)+A0(v) = 0 (4.1.2)
whereAi(v) =ai+v2bi, i= 0(1)4, a 4= 1, b 4= 0.
If we have as :
a3=−2, a 2= 2, a 1=−1, a 0= 0,
b3=17671
12096,b2=−3937
2016,b1=20483
4032,b0=−12629
3024(4.1.3)
in the form (4.1.1) , this is Quinlan-Tremaine [ ?] multistep symmetric
method with eight steps, eight order of phase-lag, eight algebraic o rder and
interval of periodicity (0 ,u2
0), whereu2
0= 0.52.
Takingaicoefficients from (4.1.3) and by wiping out the first and the second
derivative of the phase-lag , we get from (4.1.1):
a∗
3=−2, a∗
2= 2, a∗
1=−1, a∗
0= 0,
b∗
1=5
2−b∗
3−b∗
2−1
2b∗
0, b∗
0=A
B, b∗
2=C
D, b∗
3=E
F(4.1.4)
where
A=−6+140·sin(v)·v·cos(v)4+160·sin(v)·cos(v)5·v
−96·sin(v)·cos(v)7·v+2·v·cos(v)·sin(v)
−60·v·cos(v)3·sin(v)−134·v·cos(v)2·sin(v)
+32·cos(v)8·v2+20·sin(v)·v−36·cos(v)·v2
+88·cos(v)3·v2+24·cos(v)2·v2+15·cos(v)·v4
+20·cos(v)4·v2+20·cos(v)4·v4+30·cos(v)2·v4
+25·cos(v)3·v4−68·cos(v)5·v2+16·cos(v)7·v2
−64·cos(v)6·v2+18·cos(v)−120·cos(v)4+36·cos(v)5
+192·cos(v)6−12·v2+30·cos(v)2−54·cos(v)3−96·cos(v)8
+10·cos(v)5·v4−32·sin(v)·cos(v)6·v
B= 2v4·(cos(v)5−cos(v)4−2·cos(v)3+2·cos(v)2+cos(v)−1)
10

C=−6+48·sin(v)·cos(v)5·v+240·sin(v)·v·cos(v)4
−48·v·cos(v)3·sin(v)−126·v·cos(v)2·sin(v)
−6·v·cos(v)·sin(v)+20·sin(v)·v−40·cos(v)·v2
+128·cos(v)3·v2+16·cos(v)2·v2+15·cos(v)·v4
−8·cos(v)4·v2+30·cos(v)2·v4+15·cos(v)3·v4
−136·cos(v)5·v2+18·cos(v)−84·cos(v)4+192·cos(v)5
+48·cos(v)6−8·v2+42·cos(v)2−114·cos(v)3−96·cos(v)7
+48·cos(v)7·v2−128·sin(v)·cos(v)6·v
D= 4·v4·sin(v)4·(cos(v)−1)
E=−48·cos(v)6+48·cos(v)6·v2−80·sin(v)·cos(v)5·v
−48·cos(v)5·v2+48·cos(v)5+80·sin(v)·v·cos(v)4
+72·cos(v)4−96·cos(v)4·v2−78·cos(v)3+96·cos(v)3·v2
+104·v·cos(v)3·sin(v)−102·v·cos(v)2·sin(v)
+48·cos(v)2·v2−18·cos(v)2+5·cos(v)2·v4
−18·v·cos(v)·sin(v)+30·cos(v)−48·cos(v)·v2
+10·cos(v)·v4−6+5·v4+16·sin(v)·v
F=−8·v4·sin(v)4·(cos(v)−1)
By using the Taylor expansion on the previous formulas , for small va lues of
v,we accept cancellations and so we have :
b∗
0=−12629
3024+45767
12096v2−9837221
7983360v4+153204313
653837184v6−2356782689
87178291200v8
+20347993339
9700566220800v10−8744186458121
77410518441984000v12+133502728560739
28100018194440192000v14
−2016098025337469
15511210043330985984000v16+456680883838857389
84691206836587183472640000v18
b∗
2=−3937
2016+45767
40320v2−8607
39424v4+51408821
2724321600v6−35318011
34871316480v8
+3348191339
118562476032000v10−56104711163
43667471941632000v12−1538176483573
31222242438266880000v14
−1555777699603
202760915599097856000v16−14727745335969683
16606118987566114406400000v18
and
b∗
3=17671/12096
−45767
241920v2+22153
4561920v4−41092123
130767436800v6−7321421
348713164800v8
−5642643317
2134124568576000v10−210863655707
681212562289459200v12−364884558191
10035720783728640000v14
−264125909808473
62044840173323943936000v16−840723413884952309
1693824136731743669452800000v18
11

From (4.1.4) we get the explicit method .
The behavior of the coefficients b∗
0,b∗
2andb∗
4for several values of v=th
is presented in Figure 1.
Figure 1: The behavior of b∗
0,b∗
2andb∗
4for several values of v=th.
4.2. The corrector (implicit method)
For the construction of the corrector method , we set λ= 8 into (4.0.1)
and we take :
z4=−z−4−a3(z3+z−3)−a2(z2+z−2)−a1(z1+z−1)+
h2(b4(f∗
4+f−4)+b3(f3+f−3)+b2(f2+f−2)+b1(f1+f−1)+b0f0)(4.2.1)
wherezi=q(x+ih),fi=f(x+ih,q(x+ih)) andf∗
4=f(x+4h,q∗
4).
Characteristic equation (2.0.5) now becomes:
4/summationdisplay
i=1Ai(v)(wi+w−i)+A0(v) = 0 (4.2.2)
where
Ai(v) =αi+v2βi, i= 0(1)4, α 4= 1.
12

Keeping the same aicoefficients (4.1.3) and by wiping out the first and the
second derivative of the phase-lag , from the form (4.1.1) we get:
a3=−2, a 2= 2, a 1=−1, a 0= 0,
b1=20483
4032+4b4+5
3b3−2
3b0,
b2=−3937
2016−5b4−8
3b3+1
6b0
b0=G
H, b 3=I
J, b 4=K
L(4.2.3)
where
G= 15876−196560·sin(v)·v·cos(v)4−296352·sin(v)·cos(v)5·v
−152029·cos(v)3·v4−35251·cos(v)6·v4−212998·cos(v)5·v4
−289829·cos(v)4·v4+72576·sin(v)·cos(v)7·v
−24192·sin(v)·cos(v)6·v+72576·cos(v)7+217728 ·cos(v)8
+241920 ·v·cos(v)3·sin(v)−24948·v·cos(v)·sin(v)
+271404 ·v·cos(v)2·sin(v)−43092·cos(v)+362880 ·cos(v)4
+56700·cos(v)·v2+18900·cos(v)3·v2+79380·cos(v)2·v2
−31873·cos(v)·v4−60480·cos(v)4·v2−91064·cos(v)2·v4
−75600·cos(v)5·v2−43848·sin(v)·v−208656·cos(v)5
−489888·cos(v)6−106596·cos(v)2+179172 ·cos(v)3+19244·v4
−30240·cos(v)6·v2+11340·v2
H= 3024·v4·(cos(v)6−2·cos(v)5−cos(v)4+4·cos(v)3−cos(v)2
−2·cos(v)+1)
I= 18144−145152·sin(v)·v·cos(v)4+96768·sin(v)·cos(v)5·v
−59158·cos(v)3·v4+17671·cos(v)5·v4−53013·cos(v)4·v4
−193536·v·cos(v)3·sin(v)+42336 ·v·cos(v)·sin(v)
+241920 ·v·cos(v)2·sin(v)−90720·cos(v)−217728·cos(v)4
+117936 ·cos(v)·v2−93744·cos(v)3·v2−63504·cos(v)2·v2
−49233·cos(v)·v4+72576·cos(v)4·v2−56638·cos(v)2·v4
−24192·cos(v)5·v2−42336·sin(v)·v−145152·cos(v)5
+145152 ·cos(v)6+54432·cos(v)2+235872 ·cos(v)3+18931·v4
−9072·v2
J=−12096·v4·(cos(v)5−3·cos(v)4+2·cos(v)3
+2·cos(v)2−3·cos(v)+1)
13

K= 288·cos(v)6−96·cos(v)6·v2+288·sin(v)·cos(v)5·v
−288·cos(v)5+192·cos(v)5·v2−432·cos(v)4+96·cos(v)4·v2
−480·sin(v)·v·cos(v)4+468·cos(v)3−444·cos(v)3·v2
−336·v·cos(v)3·sin(v)−125·cos(v)3·v4
+588·v·cos(v)2·sin(v)+108·cos(v)2+36·cos(v)2·v2
−215·cos(v)2·v4−180·cos(v)+252·cos(v)·v2
+12·v·cos(v)·sin(v)−55·cos(v)·v4+36+35 ·v4
−36·v2−72·sin(v)·v
L= 96·v4·sin(v)4·(cos(v)2−2·cos(v)+1)
wherev=th,qis the frequency.
By using the Taylor expansion on the previous formulas , for small va lues of
v,we accept cancellations and so we have :
b0=45767
10368+58061
152064v2−182872531
2490808320v4+219498427
31384184832v6
−31166250649
106706228428800v8+562198948603
48658040163532800v10+34124493178829
224800145555521536000v12
+88156974516427
1938901255416373248000v14+24877817270533589
4839497533519267627008000v16
+1264420372774109977
2032588964078092403343360000v18+1935985026451356263
26197813314784302087536640000v20
+552400667876709828689
63474543997162912978934169600000v22+157465715679401807715941593
154997219495792173944999770062848000000v24
b3=−45767
90720−58061
1330560v2+60053897
43589145600v4−5734501
156920924160v6
−651395527
266765571072000v8−11665883797
42575785143091200v10−14989282592857
562000363888803840000v12
−5795535181309
2215887149047283712000v14−231843167482133
891486387753549299712000v16
−531006438529013
20325889640780924033433600v18−15707159026192169
5954048480632795928985600000v20
−11269380349983157813217
42210571758113337130991222784000000v22−455577298477649913331231
16847523858238279776630409789440000000v24
14

b4=45767
725760+58061
10644480v2+88852949
174356582400v4+81007601
1569209241600v6
+857181503
152437469184000v8+2185407102427
3406062811447296000v10+168261172258691
2248001455555215360000v12
+170514753691237
19389012554163732480000v14+69983291279922121
67752965469269746778112000v16
+13136856210243949
108694597009523657932800000v18+3685532027088797737
261978133147843020875366400000v20
+137450787313695211229171
84421143516226674261982445568000000v22+290371063645984514944193479
1549972194957921739449997700628480000000v24
From (4.2.3) we get the implicit method.
The behavior of the coefficients b0,b3andb4for several values of v=this
presented in Figure 2.
Figure 2: The behavior of b0,b3andb4for several values of v=th.
4.3. The development of our new (EPC2M) method
In this subsection we present the methodology we have followed to a chieve
the development of our new method.
Initially we suppose that coefficients b∗
i,bi,i= 0(1)λof (3.0.2), depend on
v,(b∗
i=b∗
i(v),bi=bi(v)).
15

Thenfrom(3.1.5)each βibecomes:βi=bi−b∗
i=bi(v)−b∗
i(v) =βi(v), i=
0(1)λ.
Replacing the above coefficients to (3.1.6) we get:
q∗
n+λ=−λ−1/summationtext
i=0a∗
izn+i+h2λ−1/summationtext
i=0b∗
i(v)f(xn+i,zn+i)
zn+λ=q∗
n+λ+h2λ/summationtext
i=0βi(v)f(xn+i,zn+i)

(4.3.1)
where
|a∗
0|+|b∗
0(v)| /negationslash= 0,|a0|+|β0(v)| /negationslash= 0
andβi(v) =bi(v)−b∗
i(v), i= 0(1)λwithb∗
λ(v) = 0.
Ifa∗
i=a∗
λ−i,b∗
i(v) =b∗
λ−i(v) andβi(v) =βλ−i(v),i= 0(1)⌊λ
2⌋the method
is symmetric.
The symmetric method (EPC2M) of eight-steps resulting form (4.3.1 ) for
λ= 8 is :
q∗
4=−(z−4+a∗
3(z3+z−3)+a∗
2(z2+z−2)+a∗
1(z1+z−1)
+a∗
0z0)+h2(b∗
3(v)(f3+f−3)+b∗
2(v)(f2+f−2)+b∗
1(v)(f1+f−1)
+b∗
0(v)f0)
z4=q∗
4+h2/parenleftBig
β4(v)(f4+f−4)+β3(v)(f3+f−3)
+β2(v)(f2+f−2)+β1(v)(f1+f−1)+β0(v)f0/parenrightBig


(4.3.2)
wherezi=q(x+ih),fi=f(x+ih,q(x+ih)),i=−4(1)3, and f4=
f(x+4h,q∗
4) and the step length is h.
Characteristic equation (2.0.5) becomes :
4/summationdisplay
i=1Ai(v)(wi+w−i)+A0(v) = 0 (4.3.3)
whereAi(v) =a∗
i+v2(βi(v)−a∗
iβ4(v))−v4b∗
i(v)β4(v),i= 0(1)4,a∗
4= 1,
b∗
4(v) = 0.
Now from (4.1.4) and (4.2.3),we can take the coefficients :
16

βi(v) =bi(v)−b∗
i(v), i= 0(1)4 which are specifically:
β0(v) =b0(v)−b∗
0(v) =b0−b∗
0,
β1(v) =b1(v)−b∗
1(v) =10403
4032+4b4+5
3b3−2
3b0+b∗
3+b∗
2+1
2b∗
0,
β2(v) =b2(v)−b∗
2(v) =b2−b∗
2=−3937
2016−5b4−8
3b3+1
6b0−b∗
2,
β3(v) =b3(v)−b∗
3(v) =b3−b∗
3
β4(v) =b4(v)−b∗
4(v) =b4(v) =b4(4.3.4)
Thus from (4.3.2), (4.1.4) and (4.3.4) we have the new eight-step sym metric
(EPC2M) method:
z∗
4=−z−4+2(z3+z−3)−2(z2+z−2)+(z1+z−1)
+h2/parenleftBig
b∗
3(v)(f3+f−3)+b∗
2(v)(f2+f−2)
+(5
2−b∗
3(v)−b∗
2(v)−1
2b∗
0(v))(f1+f−1)+b∗
0(v)f0/parenrightBig
z4=z∗
4+h2/parenleftBig
(b4)(f∗
4+f−4)+(b3−b∗
3)(f3+f−3)
+(−3937
2016−5b4−8
3b3+1
6b0−b∗
2)(f2+f−2)
+(10403
4032+4b4+5
3b3−2
3b0+b∗
3+b∗
2+1
2b∗
0)(f1+f−1)
+(b0−b∗
0)f0/parenrightBig

(4.3.5)
where
zi=q(x+ih),fi=f(x+ih,q(x+ih)),f∗
4=f(x+ 4h,q∗
4) , andb∗
0,b∗
2,b∗
3
and alsob0,b3,b4are appropriately defined by the foregoing formulas (4.1.4),
(4.1.5), (4.2.3) and (4.2.4) respectively.
To compute the (LTE) local truncation error we have to express y±λ, λ=
1(1)4 andf±κ, κ= 0(1)4 via Taylor series and with coefficients of (4.1.4) we
received the following expansion :
L.T.E.=12506213339
5794003353600y(12)
nh12+O(h14) (4.3.6)
Eventually we managed and achieved the construction of the new me thod
(4.3.5). This method (4.3.5) is an eight-step symmetric (EPC2M) meth od
phase-fitted of algebraic order ten.
Because of the necessity we had to know the first approach steps , we used
for that purpose the method Runge-Kutta Fehlberg.
17

4.4. Interval of Periodicity -Stability
As mentioned in a previous section Lambert and Watson [ ?] proved
that the method (1.0.2) has a non-vanishing interval of periodicity o nly if
the order of P-stability cannot exceed 2, if it is symmetric,and furth er more
the method is implicit. Fukushima [ ?] gave a necessary condition,Theorem
1, about when a method will have a non-vanishing interval of periodic ity.
To determine the interval of periodicity and for our stability analysis we
considered the problem :
z′′=−δ2z (4.4.1)
whereδ/negationslash=q. Applying the eight-step symmetric EPCM (4.3.2) to the model
equation (4.4.1), we take the following difference equation:
F4(w,n) (zn+4+zn−4)+F3(w,n) (zn+3+zn−3)+
F2(w,n) (zn+2+zn−2)+F1(w,n) (zn+1+zn−1)+F0(w,n)zn= 0(4.4.2)
Now the associated characteristic equation is:
F4(w,n)/parenleftbig
λ8+1/parenrightbig
+F3(w,n)/parenleftbig
λ7+λ/parenrightbig
+F2(w,n)/parenleftbig
λ6+λ2/parenrightbig
+F1(w,n)/parenleftbig
λ5+λ3/parenrightbig
+F0(w,n)λ4= 0 (4.4.3)
where
F4(w,n) = 1, F3(w,n) =a∗
3+w2(β3−a∗
3β4)−w4b∗
3β4,
F2(w,n) =a∗
2+w2(β2−a∗
2β4)−w4b∗
2β4,
F1(w,n) =a∗
1+w2(β1−a∗
1β4)−w4b∗
1β4,
F0(w,n) =w2β0−w4b∗
0β4(4.4.4)
andw=δh,n=qh. From Definition 1 and Definition 2 we have a Remark
when a (ρ,σ) method , associated with (4.3.2) method, has a nonempty
stability area.
Remark 1. From the characteristic equation (4.4.3) if the roots satis fy :
|λj| ≤1, j= 1(1)6. (4.4.5)
then the EPCM (4.3.2) has a nonempty stability area. This are a is also called
w−ndomain.
For the calculation and the graph of the s−νdomain we followed the
following sequence of steps :
18

•Calculation a loop for w(1≤w≤M) andn(1≤n≤L) respectively,
whereM,Lis defined by the user.
•Solution of the characteristic equation (4.4.3) for each set ( w,n) deter-
mined by the above loops.
•Checking any solution of the equation (4.4.3) ,the case of the satisfa c-
tion of the conditions (4.4.5).
•If specific set ( w,n) satisfy conditions (4.4.5), the corresponding point
(w,n) is plotted.
In Figure 3 where the w−ndomain is presented ,we take following the
steps above.
Figure 3: The w−vdomain of the (EPC2M) (4.3.2)
We note that in Figure 3 about the w−ndomain area :
in coloured subdomain: the numerical pair characterized as stable.
in white subdomain: the numerical pair characterized as unstable.
Also we note that : The coloured subarea presented above will be us ed
for real problems where δ/negationslash=q, ignoring the area which is around the first
diagonal. The area around the first diagonal will be used for real pr oblems
withδ=qlike the Schr¨ odinger equation.
For real problems with δ=q, to calculate the interval of periodicity, we
follow these steps:
19

1. Examining where δ=q⇔w=nso we determined the stability
polynomials Fj(w,n), j= 0(1)4
2. Exploration the area which is around of the first diagonal of the w−n
domain.
Sotheintervalofperiodicityoftheproposed(EPC2M)(4.3.2)is(0 ,n2
0)where
n2
0= 6.235009.
5. Numerical results
5.1. Problems
We compare our method (4.3.2) with eight initial value problems with
oscillating solutions to measure the efficiency of our new method. In o rder
to estimate the frequency of those problems , we bring each proble m into the
formy′′=My+NwhereM,Nare matrices. As starting point when we
use integration process we use x0≈0.1.
5.1.1. Problem by Bettis and Stiefel
Bettis and Stiefel problem can be written in the desired form as :
z′′+z= 0.001eit, z(0) = 1, z′(0) = 0.9995i, zq∈ C,(5.1.1)
or respectively :
p′′+p= 0.001 cos(t), p(0) = 1, p′(0) = 0,
r′′+r= 0.001 sin(t), r(0) = 0, r′(0) = 0.9995(5.1.2)
with theoretical solution :
z(t) =p(t)+ir(t), p,r∈ R
p(t) = cos(t)+0.0005tsin(t),
r(t) = sin(t)−0.0005tcos(t).
We solve the system of equations (5.1.2) for t∈[0,1000π].
The estimated frequency is: w= 1 (see [ ?] and [?]).
20

5.1.2. Duffing’s Equation
Duffing’s equation is given by :
¨x+δ˙x+αx+βx3=γcos(ωt) (5.1.3)
For this problem we use :
z′′=−z−z3+0.002cos(1.01t), z(0) = 0.200426728067 , z′(0) = 0
(5.1.4)
witht∈[0,1000π].
The above equation has theoretical solution :
z(t) = 0.200179477536cos(1 .01t) + 2.46946143 ·10−4cos(3.03t) + 3.04014·
10−7cos(5.05t)+3.74·10−10cos(7.07t)+….
The estimated frequency is: w= 1 (see [ ?]).
5.1.3. The N-body problem
Thisproblemstudiesthemotionofaclusterofcelestialbodiesthatin teract
with each other due to gravity.
The following formula gives us the acceleration that one body will have
due to the others bodies :
¨qi=−Gn/summationdisplay
j=1,j/negationslash=imj
q3
jiqji (5.1.5)
whereGis the gravitational constant, the mass of body jismjand the
position of body iis a vector qi.
Until now there has been no analytical solution of the n-body equat ion.
Numerical integration is the most popular method to approach the s olution
to the above problem. In order to predict the velocity and the posit ion at a
new time, we can integrate the formula (5.1.5) having regard to the p osition
and the initial conditions for each body.
We test the new (EPC2M) method for the interval t∈[0,105] .We use
value 0.00145044732989 (see [ ?]) for the frequency.
5.1.4. Problem by Franco and Palacios
Franco and Palacios studied the′almost′periodic orbit problem :
f′′+ω2f=−εf (5.1.6)
21

whereε=λ2−ω2.
This periodic orbital problem [ ?] can be written as :
f′′+f=εeiϕt, f(0) = 1, f′(0) =i, f∈ C,(5.1.7)
or respectively :
k′′+k=εcos(ϕt), k(0) = 1, k′(0) = 0,
l′′+l=εsin(ϕt), l(0) = 0, l′(0) = 1,(5.1.8)
whereε= 0.001 andϕ= 0.01.
The theoretical solution of (5.1.7) is :
f(t) =k(t)+il(t), k,l∈ R
k(t) =1−ε−ϕ2
1−ϕ2cos(t)+ε
1−ϕ2cos(ϕt)
l(t) =1−εϕ−ϕ2
1−ϕ2sin(t)+ε
1−ϕ2sin(ϕt)
We can solve that problem as a single equation in complex arithmetic or a s
a pair of coupled equations. We solve it as a system of equations (5.1.8 ) for
t∈[0,1000π].
The estimated frequency is: w= 1 (see [ ?] and [?]).
5.1.5. Schr¨ odinger’s equation
An important problem in nuclear physics isthe movement offree elect rons.
We know that these electrons move in well-defined orbits, around th e nucleus
and in a potential field created by positively charged ions as well as th e other
electrons. Every potential field is determined by parameters such as width,
depth and slope of the potential.The Woods-Saxon potential is given as :
V(r) =−V0
1+e(r−R0
a)(5.1.9)
whereR0is the potential width, V0is the depth of the potential and ais
the surface thickness. The form (5.1.9) of the potential is a funct ion of the
distance r from the center of the nucleus. In the spherical coord inates the
Schr¨ odinger equation is :
(−/planckover2pi12
2m▽2+V(r))ψ(r) =Eψ(r) (5.1.10)
22

or the radial part of it , it is :
R′′(r)+2m
/planckover2pi12[E+V0
1+qe2ar]R(r) = 0 (5.1.11)
where the radial wave function ψ(r) =R(r)/rand inserting a real constant
q we can make conversions as r−R0≡r,1
a≡2a. Here we assume that
ψ(r) = (1
r)R(r) is bounded as r→0. Formula (5.1.11) can be written as:
z′′(x) =/parenleftbiggl(l+1)
x2+V(x)−E/parenrightbigg
z(x) (5.1.12)
whereW(x) =l(l+1)
x2+V(x) is the effective potential. For our convenience,we
assume thatlim
x→∞V(x) = 0 and therefore lim
x→∞W(x) = 0.
Dividing [0,∞) into subintervals [ ai,bi] so thatW(x) is a constant Wiand
ifE >0 then we get the approximation for problem (5.1.11) :
z′′
i= (W−E)zi (5.1.13)
with solution :
zi(x) =Cie/parenleftbigg√
W−Ex/parenrightbigg
+Die/parenleftbigg
−√
W−Ex/parenrightbigg
,
Ci, Di∈R.(5.1.14)
For the needs of our work we are using another form of the Woods- Saxon
potential:
V(x) =k0
1+p+k1p
(1+p)2, p= exp/parenleftbiggx−x0
a/parenrightbigg
, (5.1.15)
and ifl= 0 at the interval [0 ,15] then we integrate problem (5.1.12) using
ask0=−50, a= 0.6, x0= 7 and k1=−k0
a. AsV(x) decays faster
thanl(l+1)
x2, the Schr¨ odinger equation (5.1.12) for large xbecomes :
z′′(x) =/parenleftbiggl(l+1)
x2−E/parenrightbigg
z(x) (5.1.16)
The radial part of the equation looks tough but the above equation has two
solutions ,which is a combination of the spherical Bessel functions an d the
23

spherical Neumann functions , qjbj(qx) andqjnj(qx) [?] .
Ifx→ ∞then the solution have the asymptotic form :
z(x)≈qjbj(qx)−qjnj(qx)
≈D[sin(qx−πl/2)+tan(sl) cos(qx−πl/2)](5.1.17)
where
tan(sl) =z(xi)S(xi+1)−z(xi+1)S(xi)
z(xi+1)C(xi)−z(xi)C(xi+1)(5.1.18)
andslis the scattering phase shift as S(x) =kjbj(kx),C(x) =kjnj(kx) and
xi<xi+1.
5.1.6. Inhomogeneous Equation
We can write the above problem in the desired form as :
a2¨x(t)+a1˙x(t)+a0x(t) =y(t) (5.1.19)
with general solution :
kx1(t)+Lx2(t)+particular integral (5.1.20)
whereKandLobtained from initial conditions and particular integral cal-
culated by special methods. We test our new (EPC2M) in problem:
z′′=−100q+99 sin(t), z(0) = 1, z′(0) = 11t∈[0,1000π] (5.1.21)
with theoretical solution: z(t) = sin(t)+sin(10t)+cos(10t).
The estimated frequency is: w= 10 (see [ ?]).
5.1.7. Kepler problem
Kepler orbit is the motion of one body relative to another, as a parab ola ,
ellipse or hyperbola , which forms a two-dimensional orbital plane in th ree-
dimensional space. In the two-dimensional Kepler problem we consid er only
the point-like gravitational attraction of two bodies and no other f orces. We
try to determine the paths of two bodies that are mutually depende nt on
each other only because of gravity. This problem can be transform ed in the
desired form as :
z′′=−z
(z2+q2)3
2, q′′=−q
(z2+q2)3
2, (5.1.22)
24

withz(0) = 1−g, z′(0) = 0, q(0) = 0, q′(0) =/radicalBig
1+g
1−g, t∈[0,1000π],where
g is the eccentricity.
Theoretical solution is :
z(t) = cos(s)−g, q(t) =/radicalbig
1−g2sin(s) (5.1.23)
wheresis the solution of the equation s−g sin(s)−t= 0.
For the problem (see [ ?] and [?]) we use as frequency the assessment
w=1
(z2+q2)3
4.
5.1.8. Nonlinear Equation
We can transform second order equations into first order equatio ns only
when the function z is missing or when the independent variable t is miss ing.
In our problem the independent variable t is missing so we have :
z′′=−100z+ sin(z), z(0) = 0, z′(0) = 1t∈[0,20π].(5.1.24)
Because we have no information on the theoretical solution we use z(20π) =
3.92823991 ·10−4.The estimated frequency is: w= 10 (see [ ?]).
5.2. The methods
To achieve the numerical integration, we used the following methods :
•The symmetric 8-step method of Quinlan-Tremaine [ ?] of algebraic
order eight: ”Q-T 8step (QT8)”
•The symmetric 8-step EPCM (phase-fitted) of Panopoulos and Simo s
[?] of algebraic order ten: ”EPCM 8-step PF (PS)”
•The symmetric 8-step EPCM of Stasinos and Simos [ ?] of algebraic
order ten with vanished phase-lag and its first derivative : ”EPCM SS
(SS)”
•The new symmetric 8-step EPC2M (4.3.5) of algebraic order ten with
vanished phase-lag and its first and second derivatives: ”New EPCM
(New)”
25

5.3. Comparison
We studied the accuracy of the new (EPC2M) method with regard to
whetherwehadknowledgeoftheoreticalsolutioninadvanceornot inrelation
to CPU time. About the methods we tested, the accuracy was expr essed by
-log10.
Table 1 shows the comparison of the above four methods for all the prob-
lems mentioned in this section.
From that table , one can understand how the new symmetric embed ded
predictor-corrector method (New EPCM (New)) (4.3.5)ismore effic ient than
all other methods.
6. Conclusions
In this work, we studied the development process of the new pair of
predictor-corrector (4.3.5). We also presented the analysis we ma de in this
new method (4.3.5) compared to the other three methods we applied to the
eight problems we mentioned in the corresponding section.
A feature of this new method is how the predictor method is totally en –
closed in the corrector, so we called it as Embedded Predictor-Corr ector
Method. The new method can be applied to IVPs with oscillatory solutio ns,
for which the frequency is known from the beginning or the frequen cy can be
estimated. It can also be used for the numerical solution for proble ms such
as the Schr¨ odinger equation or orbital problems.
After the analysis of this new method, it is easy to see how the (4.3.5)
method is more effective than the methods with which it was compared and
also that the computational cost has been reduced considerably.
In conclusion, we mention that the new (4.3.5) method is an eight step
symmetric method, withinfinite order ofphase-lag , oftenth algebr aicorder,
that has eliminated phase-lag and its first and second derivatives an d has
(0,6.235009) as interval of periodicity.
26

Table 1: Comparison for all the problems solved for Step Length and CPU time
27