Journal of Applied Mathematics

Volume 2012, Article ID 236281, 27 pages

http://dx.doi.org/10.1155/2012/236281

## A New Family of Phase-Fitted and Amplification-Fitted Runge-Kutta Type Methods for Oscillators

^{1}Department of Applied Mathematics, Nanjing Agricultural University, Nanjing 210095, China^{2}State Key Laboratory for Novel Software Technology at Nanjing University, Nanjing 210093, China^{3}Department of Computer Science and Engineering, Shanghai Jiaotong University, Shanghai 200240, China

Received 19 April 2012; Revised 17 July 2012; Accepted 23 July 2012

Academic Editor: Jesus Vigo-Aguiar

Copyright © 2012 Zhaoxia Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

In order to solve initial value problems of differential equations with oscillatory solutions, this paper improves traditional Runge-Kutta (RK) methods by introducing frequency-depending weights in the update. New practical RK integrators are obtained with the phase-fitting and amplification-fitting conditions and algebraic order conditions. Two of the new methods have updates that are also phase-fitted and amplification-fitted. The linear stability and phase properties of the new methods are examined. The results of numerical experiments on physical and biological problems show the robustness and competence of the new methods compared to some highly efficient integrators in the literature.

#### 1. Introduction

A vast of problems in applied fields, such as elastics, mechanics, astrophysics, electronics, molecular dynamics, ecology, and biochemistry, can be cast into the form of an initial value problem of the system of first-order ordinary differential equations as follows: where is a sufficiently smooth function. Traditionally, the problem (1.1) has been solved numerically by Runge-Kutta methods (RK) or linear multistep methods (LMMs) (see [1–3]). In many applications, the solution to the problem (1.1) is oscillatory or periodic. But the general-purpose methods do not take into account the oscillatory feature of the problem and the numerical results they produce are generally not satisfactory.

Recently, some authors have proposed to adapt traditional integrators to the oscillatory character of the solution to the problem (1.1) (see [4–7]). Bettis [8] constructs a three-stage method and a four-stage method which solve the equation without truncation error. Paternoster [9] develops a class of implicit methods of Runge-Kutta (RK) and Runge-Kutta-Nyström (RKN) types by the trigonometric fitting technique. For oscillatory problems which can be put in the form of a second-order equation , Franco [10] improves the update of the classical RKN methods and proposes a family of explicit RKN methods adapted to perturbed oscillators (ARKN) and a class of explicit adapted RK methods in [11]. Anastassi and Simos [12] construct a phase-fitted and amplification-fitted RK method of “almost” order five. Van de Vyver [13] investigates phase-fitted and amplification-fitted two step hybrid methods (FTSH). For other important work on frequency dependent integrators for general second-order oscillatory equations, the reader is referred to [14–18]. Many authors focus on effective numerical integration of specific categories of oscillatory problems. For example, Vigo-Aguiar and Simos [18] construct an exponentially fitted and trigonometrically fitted method for orbital problems. The papers [19–21] have designed highly efficient integrators for the Schrödinger equation.

Based on the previous work, we consider, in this paper, phase-fitted and amplification-fitted RK type integrators whose coefficients in the update depend on the product of the fitting frequency and the step size. In Section 2, we present a result on order conditions for frequency-depending modified Runge-Kutta type methods. Section 3 introduces the notion of phase-fitted and amplification-fitted RK type methods (FRK) and derives the phase-fitting and amplification-fitting conditions. In Section 4 two FRK methods of order four and two FRK methods of order five are constructed. Their error coefficients and error constants are also calculated. The linear stability and phase properties of the new FRK methods are analyzed in Section 5. In Section 6, numerical experiments are carried out to illustrate the effectiveness and superiority of our new methods compared to two well-known highly efficient integrators we have chosen from the recent literature. Section 7 is devoted to conclusive comments.

#### 2. Order Conditions for RK Type Methods with Frequency-Dependent Weights

Assume that the principal frequency of the problem (1.1) is known or can be accurately estimated in advance. This estimated frequency is denoted by , which is also called the fitting frequency. We consider the following -stage modified Runge-Kutta method: where is the step size, are real numbers, and are even functions of . The scheme (2.1) can be represented by the Butcher tableau as follows: or simply by . Conventionally, we assume that

As explained in [3], the order conditions for the modified RK method (2.1) can be derived by just considering the autonomous equation . Using Butcher's rooted tree theory in [1], the exact solution to the problem (1.1) and the numerical solution produced by the modified RK type method (2.1) have the -series expressions as follows:
where the trees , the functions (order), (number of equivalent trees), (density), (the vector of elementary weights), and the elementary differentials are defined in [3]. Then the local truncation error has the following series expansion:
If for any th differentiable function , when the scheme (2.1) is applied to the problem (1.1), the *local truncation error *, then the method (2.1) is said to have (algebraic) *order *.

Theorem 2.1 (see Franco [11]). *The modified RK type method (2.1) has order if and only if the following conditions are satisfied:
*

If the method (2.1) is of order , then we have
where is called the *error coefficient* associated with the tree of order which is involved in the leading term of the local truncation error. Denote
The positive number
is called the* error constant* of the method.

#### 3. Phase-Fitted and Amplification-Fitted Conditions

For oscillatory problems, as suggested by Paternoster [9] and Van der Houwen and Sommeijer [22], apart from the algebraic order, the analysis of phase lag and dissipation is important. To this end, we consider the following linear scalar equation: The exact solution of this equation with the initial value satisfies where . This means that after a period of time , the exact solution experiences a phase advance and the amplification remains constant.

An application of the modified RK method (2.1) to (3.1) yields
where
Thus, after a time step , the numerical solution obtains a phase advance and the amplification factor . is called the *stability function* of the method (2.1). Denote the real and imaginary part of by and , respectively. Then, for small , we have
For small , and . The above analysis leads to the following definition.

*Definition 3.1 (see [22]). *The quantities
are called the phase lag (or dispersion) and the error of amplification factor (or dissipation) of the method, respectively. If
then the method is called *dispersive of order * and *dissipative of order *, respectively. If
the method is called *phase-fitted* (or *zero-dispersive*) and *amplification-fitted* (or *zero-dissipative*), respectively. A modified RK method which is phase-fitted and amplification-fitted is called in short an *FRK method*.

It is interesting to consider the phase properties of the update of the scheme (2.1). Suppose that the internal stages have been exact for the linear equation (3.1), that is, , then the update gives where Denote the real and imaginary part of by and , respectively. Then, for small ,

*Definition 3.2. *The quantities
are called the *phase lag* (or *dispersion*) and the *error of amplification factor* (or *dissipation*) of the update of the method, respectively. If
then the update of the method is called *dispersive of order * and *dissipative of order *, respectively. If
the update of the method is called phase-fitted (or zero-dispersive) and amplification-fitted (or zero-dissipative) in the update, respectively.

Generally speaking, a traditional RK method with constant coefficients inevitably carries a nonzero phase lag and a nonzero error of amplification factor when applied to the linear oscillatory equation (3.1). Therefore, they are neither phase-fitted nor amplification-fitted. So does the update. For example, the classical RK method of order four with constant coefficients (denoted as RK4, see [3]) given by has a phase lag and an error of amplification factor as follows: For the update, Therefore, both the method and its update are dispersive of order four and dissipative of order five.

The following theorem gives the necessary and sufficient conditions for a modified RK method and its update to be phase-fitted and amplification-fitted, respectively.

Theorem 3.3. *(i) The method (2.1) is phase-fitted and amplification-fitted if and only if
** (ii) The update of the method (2.1) is phase-fitted and amplification-fitted if and only if
*

The proof of this theorem is immediate.

#### 4. Construction of New Methods

Now we proceed to construct modified RK type methods that are both phase-fitted and amplification-fitted based on the internal coefficients of two classical RK methods. For convenience we restrict ourselves to explicit methods.

##### 4.1. Fourth-Order Methods

Consider a four-stage modified RK method with the following Butcher tableau: For , the phase-fitting and amplification-fitting conditions (3.18) have the following form: On the other hand, the first-order and second-order conditions become Here, for the purpose of deriving the weights of the method, the higher order terms in the order conditions are omitted. We will follow this convention in the sequel. Solving (4.2) and (4.3), we obtain As , have the following Taylor expansions: By direct calculation, we obtain the following expansions: where the dot “” between two vectors indicates componentwise product, and and indicate the componentwise square of the vector. We will follow the convention in the sequel. Thus, the coefficients given by (4.1) and (4.4) satisfy all the conditions of order four in Theorem 2.1. But they cannot satisfy all the conditions of order five. For instance, Therefore, this method is of order four. The method was originally obtained by Simos [23] and we denote it by Simos4.

Corresponding to the nine fifth-order rooted trees , the error coefficients of Simos4 are given by Then for Simos4, where

Now we require that the update of the method (4.1) is phase-fitted and amplification-fitted. By Theorem 3.3 (ii), we have Solving (4.2) and (4.11), we obtain where It can be verified that the method given by (4.1) and (4.12) is of order four. We denote the method by FRK4.

The error coefficients of FRK4 are given by where Then for FRK4, we have

It can be seen that as , both the methods Simos4 and FRK4 reduce to a classical RK method of order four (on page 138 of [3]), which we denote by RK4. The method RK4 is called the *prototype method* or the *limit method* of the methods Simos4 and FRK4. Moreover, as , both Simos4 and FRK4 have the same error constant as that of RK4: .

##### 4.2. Fifth-Order Methods

Consider the following modified RK type method with FSAL property, the prototype of which can be found on page 167 of [2] or on page 178 of [3],

For the method (4.17), the phase-fitting and amplification-fitting conditions (3.18) become Combining these two equations and the following sufficient conditions in Theorem 2.1 corresponding the trees of orders one, two, and three, we have Solving (4.18) and (4.19), we obtain As , we obtain the following expansions:

By simple computation, we have

By Theorem 2.1, the coefficients given by (4.17) and (4.20) satisfy all the conditions of order five. But it cannot satisfy all the conditions of order six. For example, Therefore, this method is of order five and we denote it as FRK5a.

Corresponding to the twenty sixth-order rooted trees , the error coefficients of FRK5a are given by Then for FRK5a, we have where

Now we require that the update of the method (4.17) is phase-fitted and amplification-fitted. By Theorem 3.3 (ii), we have We can solve the algebraic systems (4.18), (4.27), and the last two equations in (4.19) for . The closed expressions for is too complicated. As , they have the series expressions as follows:

It can be verified that the method given by (4.17) and (4.28) is of order five. We denote the method by FRK5b.

The error coefficients of FRK5b are given by