A New Explicit Four-Step Symmetric Method for Solving Schrödinger’s Equation

Abstract

In this article we have developed a new explicit four-step linear method of fourth algebraic order with vanished phase-lag and its first derivative. The efficiency of the method is tested by solving effectively the one-dimensional time independent Schrödinger’s equation. The error and stability analysis are studied. Also, the new method is compared with other methods in the literature. It is found that this method is more efficient than these methods.

1. Introduction

In this article the numerical solution of certain second-order initial-value problems with periodical and/or oscillatory solutions as

$${\eta }^{″}\left(x\right)=f(x,\eta ),\eta \left(x_0\right)={\eta }_0\mathrm{and}{\eta }^{\prime }\left(x_0\right)={\eta }_0^{\prime },$$

(1)

is investigated. As it is seen in Equation (1), this model does not involve the first derivative explicitly. In practice, there are several mathematical models in the form of Equation (1), characterizing certain scientific problems in the fields of chemistry, quantum chemistry, physics, quantum mechanics, etc., [1].

In the literature, there are various algorithms addressing the numerical solution of boundary and initial value problems of this type. These methods can be classified in two main categories. The first category involves methods with constant coefficients, and the second one involves methods using variable coefficients. Methods of the first type are developed in Gautschi [2], Jain and et al. [3], and Steifel and Bettis in [4], while methods of the second kind are presented in Chawla and et al. [5,6], Dahlquist [7], Franco [8], Lambert and Watson [9], Krishnaiah [10], Simos and et al. [11,12], Saldanha and Achar [13], Achar [14], and Daele and Vanden Berghe [15].

A well known example of Equation (1), in science, is Schrödinger’s equation. Much research has been conducted on numerical methods for solving such a model numerically, (see for example [16,17,18,19,20,21] and the references therein).

The purpose of this article is to develop a new efficient method, with variable coefficients, for solving the time independent Schrödinger’s equation. The methodology in developing this method is based on the symmetric multistep methods introduced by Quinlan and Tremaine [22]. Its efficiency in handling numerically initial value problems with oscillatory solutions comes from the requirements imposed on its phase-lag and some of its derivatives, as well as its high algebraic order.

The article is organized as follows: Section 2 involves the basic concepts from the literature required for developing the new method. In Section 3 we present the development methodology of this method. The analysis of the error and stability of the new method is considered in Section 4. Finally, in Section 5 we present the numerical results.

2. Preliminaries

Numerical solution of the initial value problem Equation (1), can be derived by using multistep methods with k steps as:

$$\sum _{i=0}^k{\alpha }_i{y}_{n+i}=h^2\sum _{i=0}^k{\beta }_{i}f(x_{n+i},y_{n+i}),$$

(2)

using a set of equally spaced points ${\left\{x_i\right\}}_{i=0}^k$ discretizing the interval $[a,b]$ into k subintervals, with a step size h, where $h=|x_{i+1}-x_i|,\mathrm{for}i=0,1,\dots ,k-1$. When k is an even integer, ${\alpha }_{\nu }={\alpha }_{k-\nu }$ and ${\beta }_{\nu }={\beta }_{k-\nu }$ for $\nu =0,1,\dots ,\frac{k}{2}$, then the method Equation (2) is said to be symmetric method.

The linear multistep method given in Equation (2) is associated with an operator as

$$\mathcal{L}\left(x\right)=\sum _{\nu =0}^k{\alpha }_{\nu }u(x+\nu h)-h^2\sum _{\nu =0}^k{\beta }_{\nu }{u}^{″}(x+\nu h),$$

(3)

where u is any two times continuously differentiable function.

Without loss of generality we may take ${\alpha }_k=1$, hence Equation (2) can be written as

$$y_{n+k}+\sum _{i=0}^{k-1}{\alpha }_i{y}_{n+i}=h^2\sum _{i=0}^k{\beta }_{i}f(x_{n+i},y_{n+i}).$$

(4)

Definition 1.

[23] The multistep method Equation (2) is said to have an algebraic order q when the associated linear operator $\mathcal{L}$ vanishes for any algebraic polynomial in the linear span of $\{1,x,x^2,\dots ,x^{q+1}\}$.

Applying a symmetric $2m-$step method to the scalar test equation

$$y^{″}=-{\varphi }^{2}y,$$

(5)

leads to the following difference equation:

$${\mathcal{A}}_m\left(v\right)y_{n+m}+\dots +{\mathcal{A}}_1\left(v\right)y_{n+1}+{\mathcal{A}}_0\left(v\right)y_n+{\mathcal{A}}_1\left(v\right)y_{n-1}+\dots +{\mathcal{A}}_m\left(v\right)y_{n-m}=0,$$

(6)

where ${\mathcal{A}}_i\left(v\right)={\alpha }_i+{\beta }_i{v}^2,i=0,\dots ,m$, $v=\varphi h$, and h is the step size.

Equation (6) is associated with a characteristic equation given by:

$${\mathcal{A}}_m\left(v\right){\lambda }^m+\dots +{\mathcal{A}}_1\left(v\right)\lambda +{\mathcal{A}}_0\left(v\right)+{\mathcal{A}}_1\left(v\right){\lambda }^{-1}+\dots +{\mathcal{A}}_m\left(v\right){\lambda }^{-m}=0.$$

(7)

The following definition was formulated by Lambert and Watson in [9]:

Definition 2.

A symmetric $2m-$step method, which has a characteristic Equation (7), is said to have an interval of periodicity $(0,v_0^2)$ whenever the roots ${\lambda }_j,j=1,2,\dots ,,2m$ of Equation (7) satisfy:

$${\lambda }_1=e^{I\mu \left(v\right)},{\lambda }_2=e^{-I\mu \left(v\right)},and\left|{\lambda }_j\right|\le 1,for j=3,\dots ,2m,$$

(8)

for all $v\in (0,v_0^2),$ where μ is a real valued function.

Definition 3.

[24,25] The phase-lag, of any symmetric multistep method corresponding to the characteristic Equation (7), is defined as the first non vanishing term in the power series expansion of

$$\xi =v-\mu \left(v\right).$$

(9)

The phase-lag is said to have an order p whenever $\xi$ is $O\left(v^{p+1}\right)$ as $v\to \infty$. The phase-lag of any linear symmetric $2m-$step method can be computed using the following theorem proved in [24]:

Theorem 1.

Any linear symmetric $2m$-step method with characteristic equation given by Equation (7) has a phase-lag order q and a phase-lag constant c that satisfies

$$-c{v}^{q+2}+O\left(v^{q+4}\right)=\frac{{\sum }_{j=1}^{m}2{\mathcal{A}}_j\left(v\right)\cos \left(jv\right)+{\mathcal{A}}_0\left(v\right)}{{\sum }_{j=1}^{m}2j^2{\mathcal{A}}_j\left(v\right)}.$$

(10)

3. Method Development

In this article we consider an explicit linear symmetric four-step method as:

$$y_{n+2}+y_{n-2}+{\alpha }_1\left(y_{n+1}+y_{n-1}\right)+{\alpha }_0{y}_n=h^2\left[{\beta }_1\left(f_{n+1}+f_{n-1}\right)+{\beta }_0{f}_n\right],$$

(11)

where $f_{n+i}=y^{″}(x_{n+i},y_{n+i}),i=0,\pm 1$.

The following proposition follows radially from Theorem 1 with $m=2$.

Proposition 1.

The linear symmetric four-step method Equation (11) has a phase-lag order q and a phase-lag constant c given by:

$$-c{v}^{q+2}+O\left(v^{q+4}\right)=\frac{{\tau }_0\left(v\right)}{{\tau }_1\left(v\right)},$$

(12)

where

$${\tau }_0\left(v\right)=2{\mathcal{A}}_2\left(v\right)\cos \left(2v\right)+2{\mathcal{A}}_1\left(v\right)\cos \left(v\right)+{\mathcal{A}}_0\left(v\right),$$

$${\tau }_1\left(v\right)=8{\mathcal{A}}_2\left(v\right)+2{\mathcal{A}}_1\left(v\right),$$

and ${\mathcal{A}}_i\left(v\right)={\alpha }_i+{\beta }_i{v}^2,i=0,1,2$.

The parameters are chosen so that the method is consistent. Thus, setting ${\alpha }_1=-\frac{3}{40}$ and ${\alpha }_0=-\frac{37}{20}$ in Equation (11) we get the following method:

$$y_{n+2}+y_{n-2}-\frac{3}{40}\left(y_{n+1}+y_{n-1}\right)-\frac{37}{20}{y}_0=h^2\left[{\beta }_1\left(f_{n+1}+f_{n-1}\right)+{\beta }_0{f}_0\right].$$

(13)

Hence, in view of Equation (12), the method Equation (13) has a phase-lag given by:

$$PL\left(v\right)=\frac{{\mathsf{\Psi }}_1\left(v\right)}{{\mathsf{\Psi }}_2\left(v\right)},$$

(14)

where

$${\mathsf{\Psi }}_1\left(v\right):=2\cos \left(2v\right)+2\cos \left(v\right)\left(-\frac{3}{40}+{\beta }_1{v}^2\right)+{\beta }_0{v}^2-\frac{37}{20},$$

and

$${\mathsf{\Psi }}_2\left(v\right):=8+2\left(-\frac{3}{40}+{\beta }_1{v}^2\right).$$

The method in Equation (13) is subject to the following two conditions:

$$\left.\begin{array}{c}PL\left(v\right)=0,\hfill \\ P{L}^{\prime }\left(v\right)=0.\hfill \end{array}\right\}$$

(15)

Solving the linear system Equation (15) for the coefficients ${\beta }_0$ and ${\beta }_1$ we obtain

$$\left.\begin{array}{ccc}{\beta }_0\hfill & =\hfill & {\displaystyle \frac{-1}{20v^3}}\left(-160v{\cos }^2\left(v\right)+(-37+40\cos \left(2v\right))(v-2\cot \left(v\right))+6\cos \left(v\right)\cot \left(v\right)\right),\hfill \\ {\beta }_1\hfill & =\hfill & {\displaystyle \frac{1}{40v^3}}\left((3v+6\cot (v)+74\csc (v)-80\cos (2v\left)\csc \right(v)-80v\csc (v\left)\sin \right(2v)\right).\hfill \end{array}\right\}$$

(16)

As it is seen in Equation (16), the coefficients ${\beta }_i$ for $i=0,1,$ depend on $v=\varphi h$. Moreover, these coefficients have indeterminate forms as v approaches zero, as their numerators and denominators vanish as $v\to 0$. Hence, the division of such expressions will cause huge round-off errors. To overcome this situation, these coefficients are replaced by their Taylor series expansions, which are given by:

$$\left.\begin{array}{ccc}{\beta }_0\hfill & =\hfill & \frac{61}{48}+\frac{643}{2400}{v}^2-\frac{2053}{26880}{v}^4+\frac{51113}{7257600}{v}^6-\frac{23021}{63866880}{v}^8+\frac{486061}{41513472000}{v}^{10}-\frac{1727441}{5977939968000}{v}^{12}\hfill \\ \hfill & \hfill \\ \hfill & +\hfill & \frac{3495979}{1016249794560000}{v}^{14}-\frac{21785879}{108128978141184000}{v}^{16}+O\left(v^{17}\right),\hfill \\ \hfill & \hfill \\ {\beta }_1\hfill & =\hfill & \frac{637}{480}-\frac{643}{4800}{v}^2+\frac{421}{89600}{v}^4-\frac{1373}{14515200}{v}^6+\frac{161}{456192000}{v}^8-\frac{6101}{83026944000}{v}^{10}-\frac{2770787}{418455797760000}{v}^{12}\hfill \\ \hfill & \hfill \\ \hfill & -\hfill & \frac{9609433}{14227497123840000}{v}^{14}-\frac{221928661}{3243869344235520000}{v}^{16}+O\left(v^{17}\right).\hfill \end{array}\right\}$$

(17)

The behavior of the coefficients ${\beta }_0$ and ${\beta }_1$, given by Equation (16), is shown in Figure 1 and Figure 2.

4. Error and Stability Analysis

4.1. Comparative Error Analysis

The local truncation error, denoted by $LTE$, for the new developed method, named as $NL4SM$, is derived through substituting the values of the coefficients given by Equation (17) in Equation (13), then expanding both sides in Taylor series, to give:

$$LT{E}_{NL4SM}=\frac{643h^6}{9600}\left(y^{\left(6\right)}+2{\varphi }^2{y}^{\left(4\right)}+{\varphi }^4{y}^{\left(2\right)}\right).$$

(18)

To investigate the behavior of the LTE, given by Equation (18), we will use radial time independent Schrödinger’s equation

$$y^{″}\left(x\right)=\eta \left(x\right)y\left(x\right),$$

(19)

in which the function $\eta$ can be expressed as (see [17])

$$\eta \left(x\right)=\zeta \left(x\right)+G,$$

(20)

with $\zeta \left(x\right)=V\left(x\right)-V_c$, where $V\left(x\right)$ is a potential function, $V_c$ is a constant approximating the potential at a certain point x, and $G=V_c-E$, where E represents the energy.

Next, we will compare the behavior of the LTE in Equation (18) when it is applied to the above test equation with other known methods of the same type in the literature. Namely, we consider the following methods:

• The classical method, that is the method corresponding to Equation (13) when the coefficients are constants

  $$LT{E}_{CLM}=\frac{643h^6}{9600}{y}_n^{\left(6\right)}+O\left(h^8\right).$$
• The methods derived in [26]:

  $$LT{E}_{L4SMSI}=\frac{161h^6}{2400}\left(y_n^{\left(6\right)}+2w^2{y}_n^{\left(4\right)}+w^4{y}_n^{\left(2\right)}\right)+O\left(h^8\right),$$
• and

  $$LT{E}_{L4SMSII}=\frac{329h^6}{4800}\left(y_n^{\left(6\right)}+2w^2{y}_n^{\left(4\right)}+w^4{y}_n^{\left(2\right)}\right)+O\left(h^8\right).$$
• The newly developed method given by Equation (13)

  $$LT{E}_{NL4SM}=\frac{643h^6}{9600}\left(y_n^{\left(6\right)}+2{\varphi }^2{y}_n^{\left(4\right)}+{\varphi }^4{y}_n^{\left(2\right)}\right)+O\left(h^8\right).$$

To determine the asymptotic expressions of the LTEs corresponding to these methods, we have computed the derivatives occurred in the expressions of these errors, using the equation $y_n^{″}=(\zeta \left(x\right)+G)y\left(x\right)$, where $y_n$ represents an approximation of $y\left(x\right)$ at $x_n$, and we obtain

$$\begin{array}{ccc}\hfill y_n^{\left(2\right)}& =& \left(\zeta \right(x)+G)y\left(x\right),\hfill \\ \hfill y_n^{\left(4\right)}& =& {\zeta }^{″}\left(x\right)y\left(x\right)+2{\zeta }^{\prime }\left(x\right)y^{\prime }\left(x\right)+{(\zeta \left(x\right)+G)}^{2}y\left(x\right),\hfill \\ \hfill y_n^{\left(6\right)}& =& {\zeta }^{\left(4\right)}\left(x\right)y+4{\zeta }^{\left(3\right)}\left(x\right)y^{\prime }\left(x\right)+7(\zeta \left(x\right)+G)y\left(x\right){\zeta }^{″}\left(x\right)+4{\left({\zeta }^{\prime }\left(x\right)\right)}^{2}y\left(x\right)\hfill \\ & & +6(\zeta \left(x\right)+G)y^{\prime }\left(x\right){\zeta }^{\prime }\left(x\right)+{(\zeta \left(x\right)+G)}^{3}y\left(x\right).\hfill \end{array}$$

Thus, substituting these derivatives in the above error formulas, lead to the following asymptotic expansions of the LTEs for the above mentioned methods:

• Classical method

  $$LT{E}_{CLM}=h^6\left(\frac{643}{9600}y\left(x\right)G^3+\dots \right)+O\left(h^8\right).$$
• The four-step method developed in [26]

  $$LT{E}_{L4SMSI}=h^6\left[\left(\frac{161}{2400}{\left(\zeta \left(x\right)\right)}^{2}y\left(x\right)+\frac{161}{1200}{\zeta }^{{}^{\prime }}\left(x\right)y^{\prime }\left(x\right)+\frac{161}{480}{\zeta }^{{}^{″}}\left(x\right)y\left(x\right)\right)G+\dots \right]+O\left(h^8\right).$$
• The four-step method developed in [26]

  $$LT{E}_{L4SMSII}=h^6\left[\left(\frac{329}{4800}{\left(\zeta \left(x\right)\right)}^{2}y\left(x\right)+\frac{329}{2400}{\zeta }^{{}^{\prime }}\left(x\right)y^{\prime }\left(x\right)+\frac{329}{960}{\zeta }^{{}^{″}}\left(x\right)y\left(x\right)\right)G+\dots \right]+O\left(h^8\right).$$
• The new method given by Equation (11)

  $$LT{E}_{NL4SM}=h^6\left[\left(\frac{643}{9600}{\left(\zeta \left(x\right)\right)}^{2}y\left(x\right)+\frac{643}{4800}{\zeta }^{{}^{\prime }}\left(x\right)y^{\prime }\left(x\right)+\frac{643}{1920}{\zeta }^{{}^{″}}\left(x\right)y\left(x\right)\right)G+\dots \right]+O\left(h^8\right).$$

As it is seen from the above equations the LTEs involve powers of $G=V_c-E$. Thus, we will consider the following two cases:

• Case 1. $G\approx 0$, i.e., the energy, E, is close to the potential $V_c$. Therefore, in comparing the LTEs in these methods only the terms free of G are considered. Hence, for these terms all the methods are of comparable accuracy, as all their terms free of G are identical.
• Case 2. $G\ll 0\mathrm{or}G\gg 0$, that is $\left|G\right|$ is considerably large.

As it is seen from the above equations, the asymptotic error in the case of the classical method increases as the third power of G, while in the last three methods, it increases as the first power of G, but in the new constructed method it has the lowest coefficients. Thus, for the numerical integration of the radial time independent Schrödinger’s equation, it seems that the newly developed method is the most efficient, especially when $\left|G\right|$ is large.

4.2. Stability Analysis

Substituting the values of the coefficients given by Equation (16) in the method Equation (13), then applying this method to the test problem $y^{″}=-{\theta }^{2}y$, leads to a difference equation associated with the following characteristic equation

$${\lambda }^4+{\mathcal{A}}_1{\lambda }^3+{\mathcal{A}}_0{\lambda }^2+{\mathcal{A}}_1\lambda +1=0,$$

where $A_1=-\frac{3}{40}+{\beta }_1{s}^2,$ $A_0=-\frac{37}{20}+{\beta }_0{s}^2,$ and $s=\theta h$. Let us mention that the frequency $\varphi$ in the developed method Equation (13) is not equal to $\theta$, which indicates the frequency of the above test problem. The stability region of the new method is presented in Figure 3. The shaded region represents the region where the method is stable, from which it follows that the method has a periodicity interval equal to $(0,32.1489)$.

5. Numerical Results and Discussion

To test the efficiency of the new method it is applied to solve numerically the radial time independent Schrödinger’s equation:

$$y^{″}\left(\xi \right)=\left[\frac{l(l+1)}{{\xi }^2}+V\left(\xi \right)-{\gamma }^2\right]y\left(\xi \right),$$

(21)

subject to two conditions:

$$y\left(0\right)=0,$$

(22)

and one more boundary condition, imposed for large values of $\xi$, under certain physical considerations, where ${\gamma }^2$ is a real number denotes the energy, l is an integer represents the angular momentum, and V is the potential function which satisfies $V\left(\xi \right)\to 0$ as $\xi \to \infty$.

The coefficients in the new developed method are variables depending on the frequency of the problem under consideration. Therefore, we need to specify the value of the frequency $\varphi$. We will apply the new method to solve Equation (21) with $l=0$, thus we take $\varphi$ as:

$$\varphi =\sqrt{\left|V\right(\xi )-E|},$$

(23)

where $E:={\gamma }^2$ is the energy, and we will use the Woods–Saxon potential:

$$V\left(\xi \right)=\frac{50q}{a{(1+q)}^2}-\frac{50}{1+q}$$

(24)

where $q=e^{\frac{\xi -7}{a}}$ and $a=0.6$. The behavior of $V\left(\xi \right)$ is shown in Figure 4.

In Woods–Saxon potential, the parameter $\varphi$ is not presented as a function of $\xi$, instead it is estimated through approximating the potential, $V\left(\xi \right)$, at certain critical points of this potential [17]. For our numerical tests we choose $\varphi$ as follows (see [17])

$$\varphi =\left\{\begin{array}{ccc}\sqrt{50+E},\hfill & \hfill & \xi \in [0,6.5],\hfill \\ \sqrt{E},\hfill & \hfill & \xi \in [6.5,15]\hfill \end{array}\right.$$

(25)

For large values of $\xi$ the potential $V\left(\xi \right)$, in Equation (21), decays faster than $\frac{l(l+1)}{{\xi }^2}$, thus the radial Schrödinger equation is reduced to

$$y^{″}+\left({\gamma }^2-\frac{l(l+1)}{{\xi }^2}\right)=0,$$

(26)

which has two linearly independent solutions, namely, $\gamma \xi j_l\left(\gamma x\right)$ and $\gamma \xi N_l\left(\gamma \xi \right)$, where $j_l\left(\gamma \xi \right)$ and $N_l\left(\gamma \xi \right)$ are the well know Bessel and Neumann spherical functions, respectively. Hence, as $\xi \to \infty$, Equation (21) has a solution in asymptotic form which is given by:

$$\begin{array}{ccc}\hfill y& \approx & A\gamma \xi j_l\left(\gamma x\right)-B\gamma \xi N_l\left(\gamma \xi \right),\hfill \\ & \approx & AC\left[\sin \left(\gamma \xi -\frac{l\pi }{2}\right)+\tan \left({\delta }_l\right)\cos \left(\gamma \xi -\frac{l\pi }{2}\right)\right],\hfill \end{array}$$

where ${\delta }_l$ indicates the phase shift, and it can be computed using the equation

$${\delta }_l=\arctan \left(\frac{y\left({\xi }_i\right)S\left({\xi }_{i+1}\right)-y\left({\xi }_{i+1}\right)S\left({\xi }_i\right)}{y\left({\xi }_{i+1}\right)C\left({\xi }_{i+1}\right)-y\left({\xi }_i\right)C\left({\xi }_i\right)}\right),$$

(27)

where ${\xi }_{i+1}$ and ${\xi }_i$ are two different points in the asymptotic region, with $S\left(\xi \right)=\gamma \xi j_l\left(\gamma \xi \right)$ and $C\left(\xi \right)=\gamma \xi N_l\left(\gamma \xi \right)$.

Given the energy, E, we will use the developed method to solve the resonance problem, that is to approximate the phase shift, ${\delta }_l$, which has the exact value $\frac{\pi }{2}$. In order to start the method Equation (13) we need to determined the values of $y_j$ for $j=0,1,\dots ,3$. From the condition Equation (22) we have $y_0=0$. The remaining values are computed using Runge–Kutta–Nyström methods [27,28]. Based on these values, we will compute the value of ${\delta }_l$ at $x_{i+1}$ in the asymptotic region by using the following methods:

• The method of order eight produced in [22], which is denoted by $QT\left(8\right)$.
• The method of order ten developed in [22], which is denoted by $QT\left(10\right)$.
• The method of order eight developed in [22], which is denoted by $QT\left(12\right)$.
• The four-step method produced in [26], which is denoted by $SIMI$.
• The four-step method developed in [29], which is denoted by $SIM$.
• The method presented in [30], which is denoted by $SIMN$.
• The method developed in [31], which is denoted by $S2St$.
• The method produced in [20], which is denoted by $S3St$.
• The method developed in [32], which is denoted by $Y2D$.
• The method produced in [21], which is denoted by $Y3D$.
• The four-step classical method described in Section 2, which is denoted by $CLM$.
• The new four-step method presented in Section 2, as given in Equation (13), which is denoted by $NewM$.

The efficiency of the above methods is tested through computing the CPU time (in seconds) required for achieving different accuracy digits of the computed value of the phase shift, ${\delta }_l$, (which has the exact value $\frac{\pi }{2}$), for the three well known energy values; $E_1=53.588872$, $E_2=341.495874$ and $E_3=989.701916$.

Figure 5, Figure 6 and Figure 7 depict the absolute error $MER=\left|log\right(ER\left)\right|$ corresponding to the three energies $E_1,E_2$ and $E_3$ respectively, versus the values of the consumed CPU time (in seconds), where

$$ER=|\frac{\pi }{2}-{\delta }_{l,comp}|.$$

6. Conclusions

As it is seen from the results in Figure 5, Figure 6 and Figure 7, the newly developed explicit four-step method ($NewM$) is more efficient than the other methods in solving the resonance problem for the different energy values; $E_1=53.588872$, $E_2=341.495874$, and $E_3=989.701916$. It is found that among all the compared methods the new method ($NewM$) provides the highest accuracy with a lower cost of consumed CPU time.