## Abstract

In this paper we present an improved P-stable trigonometrically fitted Obrechkoff method with phase-lag (frequency distortion) infinity. Compared with the previous P-stable trigonometrically fitted Obrechkoff method developed by Simos, our new method is simpler in structure and more stable in computation. We have also improved the accuracy of the first-order derivative formula. From the numerical illustration presented, we can show that the new method is much more accurate than the previous methods.

## 1. Introduction

Numerical integration of periodic initial-value problems,(1.1)has attracted great attention for many decades. Many important equations such as the Schrödinger equation (Landau & Lifshitz 1974; Simos & Raptis 1990; Simos 1991*a*, 1992), the orbital problem (Simos 2003, 2004) and the periodic driven Duffing equation without damping (Mickens 1981; van Dooren 1974) can be placed into this category. Seeking a new method to integrate (1.1) numerically with a wide interval of periodicity, high accuracy and efficiency has been an important research activity (Allison 1970; Wang *et al*. 2004; Cash 1981; Chawla 1981; Chawla & Neta 1986; Chawla & Rao 1984, 1986, 1987; Chawla *et al.* 1986; Gautschi 1961; Gladwell & Thomas 1983; Jain 1988; Simos 1991*b*; Stiefel & Bettis 1969; Thomas 1984).

Since Lambert *et al*. introduced the idea of P-stability (Lambert & Watson 1976), it has been the criterion by which to test whether a new linear multistep method can be applied to periodic initial-value problems. Using the high-derivative envisaged by Obrechkoff (Obrechkoff 1942) opens up prospects for improving the accuracy and stability of the difference method. Ananthakrishnaiah, the first to apply Obrechkoff's idea to the numerical method for integration of (1.1), developed an adaptive method (Ananthakrishnaiah 1982), and two P-stable Obrechkoff methods with the local truncation error, *O*(*h*^{8}) and *O*(*h*^{10}), respectively (Ananthakrishnaiah 1987). Rai & Ananthakrishnaiah (1997) developed Obrechkoff methods with the additional parameters *O*(*h*^{4}), *O*(*h*^{6}) and *O*(*h*^{8}). Recently, we took advantage of the accuracy of Obrechkoff methods to solve the Schrödinger equation (Wang & Dai 2003*a*,*b*, Wang *et al*. 2004). Nevertheless, in order to satisfy the requirement of P-stable in the traditional sense, as shown in Ananthakrishnaiah (1987), we have to sacrifice accuracy. For example, for a two-step Obrechkoff method, which has from second- to sixth-order derivatives, the accuracy can be as high as *O*(*h*^{14}) for a non-P-stable method (Wang *et al*. 2004). However, for a P-stable method with the same structure as the method (2.12) in Ananthakrishnaiah (1987), the accuracy is only as low as of *O*(*h*^{8}). Simos intelligently put forward a trigonometrically fitting method to derive a P-stable Obrechkoff method using the same structure (Simos 1993). The accuracy of his method is of *O*(*h*^{12}), a little lower than the accuracy of a non-P-stable method (Wang *et al*. 2004: see (2.10) and (2.11) in §2). It seems unnecessary for an Obrechkoff method to lose accuracy to become a P-stable method. However, this trigonometrically fitted method has some small shortcomings that could be improved.

The purpose of this paper is to present a new trigonometrically fitted P-stable two-step Obrechkoff method to improve the previous one. This paper is arranged as follows: in §2 a brief review of the previous trigonometrically fitted P-stable Obrechkoff method is made; in §3 the derivation of the new method is given; in §4, numerical results from the application of the new method to some well-known problems are compared with the results obtained from the other methods; conclusions are given in §5.

## 2. Preliminaries

The Simos method is briefly reviewed in order to appreciate the advantages of the new, compared with the previous, method (Simos 1993). The following structure has been considered by Simos:(2.1)

By using the Taylor series expansion to the algebraic order 10, the following set of equations is obtained:(2.2)

By solving this system of equations, we obtain(2.3)

After substituting (2.3) into (2.1), the following test equation to the resultant expression is applied,(2.4)and the following characteristic equation is obtained,(2.5)

The characteristic roots of (2.5) are(2.6)where(2.7)

If *B*(*H*)=*A*(*H*) cos(*H*) is valid for *H*∈(0, ∞), then the characteristic roots |*λ*_{±}|=|*e*^{±iH}|≤1 can always satisfy. Therefore, *γ*_{2} can be written as(2.8)

As shown in figure 1*b*, the denominator of *γ*_{2} in (2.8) contains so many zeros that it cannot work properly for a small step size. In the previous Simos paper, the denominator of (2.8) is expanded into a Taylor series as −59/2*H*^{12}, and the numerator is modified, so that *γ*_{2} can be approximately expressed as(2.9)where

Although the zeros have been successfully eliminated by using expression (2.9), the resultant factor of *H*^{−12} will increase the error which arises from the neighbourhood around a zero numerator as shown in figure 1*a*. Hence the instability has not been solved for the small step size (see figure 6).

By using (2.9), the final local truncation error is(2.10)

We remark that if we let *γ*_{2}=−2923/3925152, then (2.2) represents the non-P-stable method (Wang *et al*. 2004). The local truncation error is(2.11)

In Simos (1993), the following first-order derivative formula was used:(2.12)

Obviously, the trigonometrically fitted P-stable method has a great advantage in terms of accuracy over the traditional method developed by Ananthakrishnaiah (1982). However, the coefficient *γ*_{2} of (2.9) contains many zeros in its numerator and denominator, which makes this method more complex in practical application. We find that it is instable to integrate a periodic driven Duffing equation without damping when a small step size is used (see figure 6). On the other hand, the P-stable condition cannot hold whenever the zeros appear in the numerator or (2.9) is used for *γ*_{2}. It can lead to a large error for a high-frequency problem when it is not a P-stable. In order to overcome these faults, a new trigonometrically fitted Obrechkoff method was put forward.

## 3. Derivation

In this section, the new method is introduced. The derivation is divided into two parts, the main structure (§3*a*) and the first-order derivative formula (§3*b*).

### (a) Main structure

The same structure, (2.1), will be used to derive the new trigonometrically fitted Obrechkoff method. The derivation procedure used by Simos is reversed. Firstly, the P-stable requirement is considered to restrict the coefficients of (2.1), and then Taylor's series expansion is used to obtain all the coefficients.

Firstly, it is assumed that (2.1) is P-stable. Therefore, when applying the test equation, (2.4) to (2.1), , and , the following characteristic equation can be obtained:(3.1)where(3.2)

From (3.1) and (3.2), *α*_{2} is allowed as(3.3)

Next, (3.3) is placed into (2.1) and expanded into a Taylor series, such that(3.4)

From (3.4), the remaining coefficients can be determined and the truncation error can be found as(3.5)

So, the main structure of the new two-step P-stable Obrechkoff's method can be written as(3.6)with the truncation error(3.7)where *α*_{2} is given by (3.3) and (3.5).

Compared with the previous method in §2, it can be seen that the new method (3.6) is simpler and more accurate, and the local truncation error is the same as that of (2.11) for the non-P-stable method. Furthermore, it is reasonable to expect that the new method will be more stable in numerical computation, because only one coefficient, *α*_{2}, depends on the step size, and the *h*^{2} factor in the denominator of *α*_{2} will be cancelled by the *h*^{2} arising from the second derivative, *y*″(*x*).

### (b) First-order derivative formula

For a complex problem, we need to calculate *y*′(*x*+*h*), which results from *y*^{(4)}(*x*+*h*) and *y*^{(6)}(*x*+*h*). Therefore, an Obrechkoff method should usually include a main structure and a first-order derivative formula. The accuracy of the first-order derivative formula should be matched with the main structure, such that the Obrechkoff method can fully display its advantage in accuracy. The following structure is used to derive our first-order derivative formula:(3.8)

The coefficients can easily be obtained by a Taylor series expansion as shown in Wang *et al*. (2004). Finally, the following explicit first-order derivative formula is obtained:(3.9)

Its local truncation error is .

## 4. Numerical illustrations

To illustrate the new method, we will apply it to three problems. Since *y*_{n}=cos(*nhω*) is the exact solution of test equation (2.4), the test equation is replaced by a test-like problem as the first problem. The remaining two are the same as in Simos (1993), the ‘almost periodic’ problem of Stiefel & Bettis (1969) and the nonlinear undamped Duffin's equation.

*Problem 1*. A test-like problem is designed that is highly oscillatory with a more complex pattern compared with the test equation,(4.1)

The exact solution to this problem is . Its complex oscillatory pattern can be seen in figure 2.

Equation (4.1) is directly integrated numerically by three methods: the new method (3.6), Simos method in §2 and method (2.11) in Ananthakrishnaiah (1987). The numerical results are shown in table 1 and figure 3. Obviously, both trigonometrically fitted P-stable methods are much better than the traditional one if the vertical scale used in the figure is noted.

*Problem 2*. We consider the following ‘almost periodic’ problem studied by Stiefel & Bettis (1969):(4.2)whose theoretical solution represents a motion of the perturbed circular orbit in the complex plane and can be written as(4.3)

It is assumed that the new method (3.6) is able to deal with complex problems; let . According to (4.2), we obtain(4.4)which are valid for any node *n*.

After substituting (4.4) into (3.6) and letting *ω*=1, the following iterative formula for *z*_{n} can be obtained(4.5)where

If letting *z*_{0}=1 and *z*_{1}=*z*(*h*), according to the exact solution (4.3), the problem can easily be solved in a numerical way by using (4.5). Following Simos (1993), *z*_{k} has been calculated by using the step size of *h*=*π*/2, *π*/3,…,*π*/6 in the range of 0<*x*<40*π* (which corresponds to 20 orbits of the point *z*(*x*)). The errors are , where *N*_{k}=Int(40*π*/*h*). For comparison purposes, in table 2, the errors are listed in the computed values of *z*(40*π*) obtained by method (4.5), method (2.11) of Ananthakrishnaiah (1987) and the Simos method in §2.

*Problem 3*. The nonlinear undamped Duffing's equation is considered:(4.6)with *B*=0.002 and *v*=1.01.

The following solution is used(4.7)(4.8)as the ‘exact’ solution of (4.7), which is derived by a numerical approach method. Its maximum error is approximately 10^{−14}, much higher than the one given by Simos (1993).

In order to show the error behaviour of (4.7), the error can be checked by the following expression:(4.9)

In figure 4, Δ*g* is plotted for the range of {0, 2000*π*}. From the figure, the maximum absolute error increases with *x* slowly. The maximum absolute error of Δ*g* is smaller than 10^{−14} within 0≤*x*≤2000*π*, which is sufficient for the purpose of the test.

The new method to integrate the Duffing equation is implemented in the following two-stage way. It is assumed that all the values of *y*_{k}, *y*_{k+1}, , ,… are known. We want to find *y*_{k+2}, , ,…. Since , and can be expressed as a function of *y*_{k+2} and , the purpose of the two-stage way is to find *y*_{k+2} and as accurately as possible. The first stage is to find the preliminary value of *y*_{k+2}. The second stage is to find and *y*_{k+2} iteratively, until the required accuracy is reached.

In the first stage, the following two formulae to calculate the preliminary value of *y*_{k+2} are used.(4.10)

The error of *y*_{k+2} can be estimated to be approximately

The following formula is used to recalculate *y*_{k+2}:(4.11)

This can further reduce the error of *y*_{k+2} because the local truncation error of the above formula is about .

In the next stage, the first-order derivative formula (3.9) and the main structure (3.6) are used to find and *y*_{k+2} iteratively. The following formula to calculate is used:(4.12)where(4.13)(4.14)with(4.15)

If *y*_{k+2} is ‘exact’, then the error in computed should be approximately . Once is obtained, all the high-order derivatives can be calculated. In order to increase efficiency, the following sequence should be used:(4.16)

In the iterative procedure, the following formula is used to calculate *y*_{k+2}. Comparing (4.9) and (4.10), the accuracy will now be greatly improved:(4.17)where(4.18)

(4.19)

(4.20)

In order to obtain the numerical result with the minimum error, the iterative procedure, from (4.12) to (4.20), should be repeated once or twice. On the other hand, it should be noted that (4.14), (4.15), (4.18) and (4.19) contain no *y*_{k+2} and , *Λ*_{0}, *Φ*_{0}, *Λ*_{1} and *Φ*_{1} only needs to be computed once.

For convenience, a simple Mathematica program based on the above algorithm for problem 3 is given in the appendix.

In table 3 the absolute errors in *y*(*x*) with *h*=*π*/5 are presented along with a list of the numerical results by the previous methods. In table 4 the absolute errors in *y*(*x*) with *h*=*π*/2, *π*/4, *π*/6, *π*/8, *π*/10 are presented.

In order to see how the error changes with the integral destination, Δ*y*_{n}=*y*_{n}*−g*(*nh*) with *h*=*π*/5 is calculated by the new method (3.6) and Simos method in §2 for problem (4.6). The plots are given in figure 5. Surprisingly, it is found that two of the plots are the same apart from the amplitude. The accuracy of the new method is higher than the Simos method with five orders. Both two error behaviours are stable and periodic and the period is about 155*π*.

According to the error periodicity shown in figure 5, if the average value of the error is defined asfor a period of 155*π*, the error behaviour and the dependence of the error on the step size can be conveniently compared. In figure 6, |Δ*y*|_{av} is plotted as a function of step size for these two methods. In order to show the two plots in the same figure, |Δ*y*|_{av}×5×10^{−6} for the Simos method is used. From the figure, it can be seen that the new method surpasses the Simos method in accuracy, by up to five to six orders overall, and the Simos method is shown to be unstable when *h*≤*π*/20, while the new method is stable even when *h* approaches *π*/200.

From these results it is easy to see that the new trigonometrically fitted P-stable method is much better than the previous one in stability and accuracy. Furthermore, the numerical data with high accuracy can be attributed to the highly accurate first-order derivative formula.

All the computations were carried out by the algorithm system Mathematica 5.0 on an IBM PC-AT with AMD Athlon XP 2600+, 512Mb memory.

## 5. Conclusion

In this paper, a new trigonometrically fitted P-stable Obrechkoff method for periodic initial-value problems is presented which improves the previous method developed by Simos in three aspects: it simplifies the structure eliminates the computation instability for small step size and increases the accuracy of the first-order derivative formula. From the numerical illustrations, it is shown that the new method is a powerful one for periodic initial-value problems.

## Acknowledgments

The authors would like to thank the Department of Physics, Shanghai University, for its continuous support during the progress of this work. This work is supported by National Natural Science Foundation of China (grant no. 60371033).

## Footnotes

- Received June 23, 2004.
- Accepted November 3, 2004.

- © 2005 The Royal Society