Implicit Exponentially Fitted RKNd Methods for Solving Oscillatory ODEs
Wenjuan Zhai^{1,}^{ *}, Bingzhen Chen^{2}
^{1}Department of Mathematics, Beijing Jiaotong University Haibin College, Cangzhou, P. R. China
^{2}Department of Applied Mathematics, Beijing Jiaotong University, Beijing, P. R. China
Email address:
To cite this article:
Wenjuan Zhai, Bingzhen Chen. Implicit Exponentially Fitted RKNd Methods for Solving Oscillatory ODEs. Science Journal of Applied Mathematics and Statistics. Vol. 4, No. 2, 2016, pp. 74-80. doi: 10.11648/j.sjams.20160402.19
Received: February 25, 2016; Accepted: March 9, 2016; Published: April 13, 2016
Abstract: In this paper, we derive the implicit exponentially fitted RKNd methods for solving oscillatory ODEs. The new methods integrate exactly differential systems whose solutions can be expressed as linear combinations of functions from the set {exp(λt), exp(−λt)}, λ ∈ C, or equivalently when λ = iω, ω ∈ R. Numerical experiments are accompanied to show the efficiency and competence of the implicit exponentially fitted RKNd methods compared with implicit RKNd methods.
Keywords: RKNd Method, Exponentially Fitted, Implicit, Stability, Efficiency, Oscillatory
1. Introduction
A large number of problems arising in physical fields such as elastics, celestial mechanics, quantum differential equations
(1)
In the begining, we use Euler method to solve these problems. But since Runge (1895) and Kutta (1905) derived the Runge-Kutta(RK) methods, we have constructed many RK methods of different orders (see Ref. [6,7]) to solve these problems. Bing-Zhen Chen and Xiong You gave a new class of method based on the analysis of the internal stage’s order of the RK method in [4]. The new method is denoted as RKNd and when it comes to the initial value problem (1), it is a scheme of the form
(2)
where Compared with the RK method, the RKNd method has two merits: (i) when they have the same stage number, RKNd method is able to achieve higher algebraic order than RK method; (ii) the internal stages k_{i} of RKNd method have the same coefficients of h^{2} with . In addition, this scheme requires the function of (1) has derivative of x. This restraint can’t be a problem in practice, because in most of the application problems the function has derivative or even higher derivative.
The RKNd method (2) can also solve the second order differential equations
(3)
If a new variable z is introduced to represent the first derivative y′, then it is turned into the partitioned system of first order equations
So, we can use the RKNd method (2) to solve second order differential equations.
In the last decade, a great interest in the research of new methods for the numerical integration of initial value problems having oscillatory or periodic solutions has arisen. A problem is called oscillatory if the solution of the problem is oscillatory. When it comes to the second order initial value problems of the form
(4)
where K is a symmetric positive semi-definite matrix (stiffness matrix) that contains implicitly the main frequencies of the problem, if the magnitude of the matrix K satisfies ≤, then the solutions of the initial value problem (4) are in general oscillatory or highly oscillatory. These initial value problems are of great interest in many problems of molecular dynamics, orbital mechanics, electronics, and so on, and they need efficient numerical methods because of the high accuracy demands.
Since Gautschi (see Ref. [8]) and Lyche (see Ref. [9]) gave the theoretical analysis of exponential fitting technic, many researchers have been working on the exponential fitting methods. A method is called exponentially fitted, that is to say, it integrates exactly differential systems whose solutions can be expressed as linear combinations of functions from the set{exp(λt), exp(−λt), λ ∈ C}, or equivalently {sin(ωt), cos(ωt)} when λ=iω, ω∈R. The construction of exponentially fitted RK(N) methods is originally due to B. Paternoster (see Ref. [10]), and a detailed exposition of exponentially fitted methods with an extensive bibliography on this subject can be found in Ixaru and Vanden Berghe [12]. And a lot of exponential fitting RK and RKN methods can be found in [3, 5, 11-14, 16]. When we solve oscillatory problems, if we can give a good estimate of oscillatory frequency, the exponential fitting method has an advantage over other methods. Therefore, it has become an indispensable tool for solving oscillatory problems.
Recently, Zhai Wenjuan and Chen Bingzhen has considered explicit EFRKNd methods in [21]. In this paper, we construct the implicit exponentially fitted RKNd (IEFRKNd) methods based on the implicit RKNd methods derived in [4]. The paper is organized as follows: In Section 1, we present the modified form of the RKNd method. Section 2 is concerned with the exponential fitting conditions of the modified RKNd method presented in Section 1. Section 3 is devoted to the algebraic order conditions of the modified EFRKNd method. In Section 4, we derive three two-stage implicit exponentially fitted RKNd methods of orders two, three and four respectively. In Section 5, we analyze the stability properties of the new methods, obtaining the stability regions. In Section 6, we present some numerical experiments to show the efficiency of the new methods compared with implicit RKNd methods given in [4].
2. Implicit Modified RKNd Methods
In this paper, we present a class of implicit modified RKNd methods which integrate exactly differential systems whose solutions can be expressed as linear combinations of the set of functions {exp(λt), exp(−λt)} or equivalently {sin(ωt), cos(ωt)} when λ = iω, ω ∈ R. This means that the internal stage equations and the final step equation have to integrate exactly these sets of functions. So, we introduce a modification of the algorithm of implicit RKNd method (s ≥ 2), given by
(5)
which can be expressed in the Butcher tableau form as
.
We note that the algorithm (5) coincides with the algorithm of RKNd method (2) when , and the remaining coefficients are constant. So, the coefficients and are introduced in the algorithm definition so that the family of functions exp{±iωt}, can be integrated exactly by the method.
3. Exponentially Fitted Conditions
Following Albrecht’s approach (see Refs. [1,2]), the each stage of algorithm can associate a linear functional as follows.
• for the internal stage
(6)
• for the final stage
(7)
Requiring that (6) and (7) will vanish for the functions from the set exp{±iωt}, we obtain the following equations:
(8)
Having in mind the relations cosh(z) = (e^{z} + e^{−z})/2 and sinh(z) = (e^{z} − e^{−z})/2, Eqs. (8) can be expressed in the form
(9)
(10)
The conditions defined by Eqs. (9) and (10) characterize when an implicit RKNd method (5) is exponentially fitted, and therefore they will be denominated as exponential fitting conditions (EF conditions).
4. Algebraic Order Conditions
Now, we make a study of the local truncation error for the IEFRKNd methods in order to obtain the order conditions for the algorithm (5).
As discussed in [21], we consider the following assumptions
(11)
Using the assumptions mentioned above and following the way given in Hairer [6, pp.143-148] for obtaining the terms of the local truncation error, the order conditions for IEFRKNd methods (up to fifth order) are the following ones:
Order 2 requires:
(12)
Order 3 requires in addition:
(13)
Order 4 requires in addition:
(14)
Order 5 requires in addition:
(15)
5. Construction of Implicit EFRKNd Methods
In this section, we analyze the construction of implicit two-stage EFRKNd methods (up to order 4) with the help of the order conditions obtained in the previous section.
In order to derive IEFRKNd methods corresponding to IRKNd methods obtained in [4], we assume the coefficients satisfy So the EF conditions (9) and (10) reduce to
(16)
(17)
From (16) and (17), we have
(18)
The taylor expansions of these coefficients are
We can verify that the coefficients are satisfied the second order condition (12). So the IEFRKNd method has algebraic order at least two.
As we can see, when z → 0, the method (18) reduces to IRKNd method in [4]. We can choose different values of to obtain different order methods. As discussed in [4], we can obtain the following different order methods.
First, we choosethen (18) becomes
We denote this method as IEFRKNd2s2. When z → 0, it reduces to the method IRKNd2s2. We can verify that the method IEFRKNd2s2 is of order 2.
Second, by choosing (18) becomes
We denote this method as IEFRKNd2s3. When z → 0, it reduces to the method IRKNd2s3. We can verify that the method IEFRKNd2s3 is of order 3.
Now we derive the fourth order method. From the order conditions (13) and (14), we obtain
So, and we obtain the fourth order method also obtain a fourth order IRKNd method
We denote this method as IEFRKNd2s4. Here, we also obtain a fourth order IRKNd method
.
6. Stability Properties and Phase Analysis
Before we put a method into application, we must analysis its stability or it will be dangerous. Roughly speaking, the stability of a numerical method is that the numerical solution has a small change when the initial value has a small perturbation. Since the coefficients of exponentially fitted method depend on z = λh, we consider its imaginary stability region (see Ref. [18]).
Applying the method (5) to the test equation
we obtain , where
θ = µh, i^{2}= −1, and R(iθ; z) is called the stability function. The region of imaginary stability is a region in the θ − z plane (θ > 0, z > 0), throughout which | R(iθ; z) |≤ 1. In Fig. 1, (a)- (e) show the stability regions for IEFRKNd2s2, IRKNd2s2, IEFRKNd2s3, IRKNd2s3 and IEFRKNd2s4 respectively.
(a)
(b)
(c)
(d)
(e)
From Fig. 1, we can see our IEFRKNd methods have larger stability regions than IRKNd methods.
7. Numerical Experiments
To test the numerical performance of the implicit EFRKNd methods, we carry out experiments on three oscillatory or periodic problems to illustrate the effectiveness and efficiency. The codes used in the comparisons are:
• IEFRKNd2s2: two-stage IEFRKNd method of order 2 derived in this paper;
• IEFRKNd2s3: two-stage IEFRKNd method of order 3 derived in this paper;
• IEFRKNd2s4: two-stage IEFRKNd method of order 4 derived in this paper;
• IRKNd2s2: two-stage IRKNd method of order 2 given in [3];
• IRKNd2s3: two-stage IRKNd method of order 3 given in [3];
The criterion used in the numerical comparisons is the decimal logarithm of the maximum global error (log10(GE)) versus the computational effort measured in the number of function evaluations required by each method.
Problem 1. Consider the second order ODE
whose analytic solution is given by y(x) = sin(30x)/30.
In our test we chose the parameter values ω = 30, λ = 30i, xend = 100, and the numerical results presented in Figure. 2(a) have been computed with the integration steps , m =1, 2, 3, 4, 5.
Problem 2. Consider the inhomogeneous equation (see Ref. [19])
y′′ + 100y = 99 sin(x), y(0) = 1, y′= 11(11),
whose analytic solution is given by
The equation has been solved in the interval [0, 10] with The numerical result is shown in Fig. 2 (b).
Problem 3. Consider the nonhomogeneous Duffing equation used in [17]
The ‘exact’ solution was first given by Van Dooren (see Ref. [20]), but a more accurate solution was given by Zhao et al. (see Ref. [22]), i.e.
with
The equation is integrated over the interval [0, 100] with , m = 1, 2, 3, 4, ω = 1.01. The numerical result is shown in Fig. 2(c).
(a)
(b)
(c)
From Fig. 2(a), we can see that when the solution of differential system can be expressed as linear combinations of functions from the set {exp(λt), exp(−λt)}, λ ∈ C, our EF IRKNd methods are more efficient than IRKNd methods which are not exponentially fitted. From Fig. 2(b) and Fig. 2(c), we can find that the efficient curves of the methods IEFRKNd2s2, IEFRKNd2s3 and IRKNd2s23 are almost same. But the method IEFRKNd2s4 is the most efficient method among them.
8. Conclusion
In this paper, we derive the implicit exponentially fitted RKNd methods of different order for solving oscillatory ODEs. These new methods integrate exactly differential systems whose solutions can be expressed as linear combinations of functions from the set {exp(λt),exp(−λt)}, λ ∈ C, or equivalently when λ = iω, ω ∈ R. We also analysis the stability of these methods. Numerical experiments are accompanied to show the efficiency and competence of the implicit exponentially fitted RKNd methods compared with implicit RKNd methods.
References