A Trigonometrically Fitted Block Method for Solving Oscillatory Second-Order Initial Value Problems and Hamiltonian Systems

Hindawi International Journal of Differential Equations Volume 2017, Article ID 9293530, 14 pages https://doi.org/10.1155/2017/9293530 Research Article A Trigonometrically Fitted Block Method for Solving Oscillatory Second-Order Initial Value Problems and Hamiltonian Systems F. F. Ngwane1 and S. N. Jator2 1 Department of Mathematics, University of South Carolina, Salkehatchie, Walterboro, SC 29488, USA Department of Mathematics and Statistics, Austin Peay State University, Clarksville, TN 37044, USA 2 Correspondence should be addressed to F. F. Ngwane; Received 25 July 2016; Revised 13 December 2016; Accepted 14 December 2016; Published 22 January 2017 Academic Editor: Julio D. Rossi Copyright © 2017 F. F. Ngwane and S. N. Jator. 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. In this paper, we present a block hybrid trigonometrically fitted Runge-Kutta-Nyström method (BHTRKNM), whose coefficients are functions of the frequency and the step-size for directly solving general second-order initial value problems (IVPs), including Hamiltonian systems such as the energy conserving equations and systems arising from the semidiscretization of partial differential equations (PDEs). Four discrete hybrid formulas used to formulate the BHTRKNM are provided by a continuous one-step hybrid trigonometrically fitted method with an off-grid point. We implement BHTRKNM in a block-by-block fashion; in this way, the method does not suffer from the disadvantages of requiring starting values and predictors which are inherent in predictor-corrector methods. The stability property of the BHTRKNM is discussed and the performance of the method is demonstrated on some numerical examples to show accuracy and efficiency advantages. 1. Introduction In what follows, we consider the numerical solution of the general second-order IVPs of the form 𝑦󸀠󸀠 = 𝑓 (𝑥, 𝑦, 𝑦󸀠 ) , 𝑦 (𝑥0 ) = 𝑦0 , (1) 𝑦󸀠 (𝑥0 ) = 𝑦0󸀠 , 𝑥 ∈ [𝑥0 , 𝑥𝑁] , where 𝑓 : R × R2𝑚 → R2𝑚 , 𝑁 > 0 is an integer, and 𝑚 is the dimension of the system. Problems of form (1) frequently arise in several areas of science and engineering such as classical mechanics, celestial mechanics, quantum mechanics, control theory, circuit theory, astrophysics, and biological sciences. Equation (1) is traditionally solved by reducing it into a system of first-order IVPs of double dimension and then solved using the various methods that are available for solving systems of first-order IVPs (see Lambert [1, 2], Hairer and Wanner in [3], Hairer [4], and Brugnano and Trigiante [5, 6]). Nevertheless, there are numerous methods for directly solving the special second-order IVPs in which the first derivative does not appear explicitly and it has been shown that these methods have the advantages of requiring less storage space and fewer number of function evaluations (see Hairer [4], Hairer et al. [7], Simos [8], Lambert and Watson, and [9], Twizell and Khaliq [10]). Fewer methods have been proposed for directly solving second-order IVPs in which the first derivative appears explicitly (see Vigo-Aguiar and Ramos [11], Awoyemi [12], Chawla and Sharma [13], Mahmoud and Osman [14], Franco [15], and Jator [16, 17]). It is also the case that some of these IVPs possess solutions with special properties that may be known in advance and take advantage of when designing numerical methods. In this light, several methods have been presented in the literature which take advantage of the special properties of the solution that may be known in advance (see Coleman and Duxbury [18], Coleman and Ixaru [19], Simos [20], Vanden Berghe et al. [21], VigoAguiar and Ramos [11], Fang et al. [22], Nguyen et al. [23], 2 International Journal of Differential Equations Ramos and Vigo-Aguiar [24], Franco and Gómez [25], and Ozawa [26]). However, most of these methods are restricted to solving special second-order IVPs in a predictor-corrector mode. Our objective is to present a BHTRKNM that is implemented in a block-by-block fashion; in this way, the method does not suffer from the disadvantages of requiring starting values and predictors which are inherent in predictorcorrector methods (see Jator et al. [27], Jator [16], and Ngwane and Jator [28]). We note that multiderivative trigonometrically fitted block methods for 𝑦󸀠󸀠 = 𝑓(𝑥, 𝑦, 𝑦󸀠 ) have been proposed in Jator [29] and Jator [16]. However, the BHTRKNM proposed in this paper avoids the computation of higher order derivatives which have the potential to increase computational cost, especially, when applied to nonlinear systems. In this paper, we propose a BHTRKNM which is of order 3 and its application is extended to solving oscillatory systems, PDEs, and Hamiltonian systems including the energy conserving equation. The organization of this article is as follows. In Section 2, we derive the BHTRKNM for solving (1). The analysis and implementation of the BHTRKNM are discussed in Section 3. Numerical examples are given in Section 4 to show the accuracy and efficiency of the BHTRKNM. Finally, the conclusion of the paper is given in Section 5. 2. Development of the BHTRKNM In order to numerical integrate (1) we define the BHTRKNM as consisting of the following four discrete formulas: 1 𝑦𝑛+1 = 𝑦𝑛 + ℎ𝑦𝑛󸀠 + ℎ2 ( ∑ 𝛽𝑗 𝑓𝑛+𝑗 + 𝛽𝑛+V 𝑓𝑛+V ) , 𝑗=0 1 𝑦𝑛+V = 𝑦𝑛 + ℎ𝑦𝑛󸀠 + ℎ2 ( ∑ 𝛾𝑗 𝑓𝑛+𝑗 + 𝛾𝑛+V 𝑓𝑛+V ) , 𝑗=0 (2) 1 󸀠 󸀠 = ℎ𝑦𝑛󸀠 + ℎ2 ( ∑ 𝛽𝑗󸀠 𝑓𝑛+𝑗 + 𝛽𝑛+V 𝑓𝑛+V ) , ℎ𝑦𝑛+1 𝑗=0 derive the main method and additional methods we initially seek a continuous local approximation Π(𝑥) on the interval [𝑥𝑛 , 𝑥𝑛+1 ] of the form Π (𝑥) = 𝛼0 (𝑥) 𝑦𝑛 + 𝛿0 (𝑥) ℎ𝑦𝑛󸀠 1 + ℎ2 ( ∑𝛽𝑗 (𝑥) 𝑓𝑛+𝑗 + 𝛽𝑛+V (𝑥) 𝑓𝑛+V ) , (4) 𝑗=0 where 𝛼0 (𝑥), 𝛿0 (𝑥), and 𝛽𝑗 (𝑥), 𝑗 = 0, V, 1, are continuous coefficients. The first derivative of (4) is given by Π󸀠 (𝑥) = 𝑑 Π (𝑥) . 𝑑𝑥 (5) We assume that 𝑦𝑛+𝑗 = Π(𝑥𝑛+𝑗 ) is the numerical approxima󸀠 tion to the analytical solution 𝑦(𝑥𝑛+𝑗 ), 𝑦𝑛+𝑗 = Π󸀠 (𝑥𝑛+𝑗 ) is the numerical approximation to 𝑦󸀠 (𝑥𝑛+𝑗 ), and 𝑓𝑛+𝑗 = Π󸀠󸀠 (𝑥𝑛+𝑗 ) is an approximation to 𝑦󸀠 (𝑥𝑛+𝑗 ), 𝑗 = 0, V, 1. The following theorem shows how the continuous method (4) is constructed. This is done by requiring that on the interval from 𝑥𝑛 to 𝑥𝑛+1 = 𝑥𝑛 + ℎ the exact solution is locally approximated by function (4) with (5) obtained as a consequence. Theorem 1. Let 𝐹𝑖 (𝑥) = 𝑥𝑖 , 𝑖 = 0, 1, 2, 𝐹3 (𝑥) = sin 𝑤𝑥, and 𝐹4 (𝑥) = cos 𝑤𝑥 be basis functions and let 𝑉 = (𝑦𝑛 , 𝑦𝑛󸀠 , 𝑓𝑛 , 𝑓𝑛+V , 𝑓𝑛+1 )𝑇 be a vector, where 𝑇 is the transpose. Define the matrix 𝐺 by 𝐹0 (𝑥𝑛 ) ⋅⋅⋅ 𝐹4 (𝑥𝑛 ) 𝐹0󸀠 (𝑥𝑛 ) ( 󸀠󸀠 𝐺=( ( 𝐹0 (𝑥𝑛 ) 𝐹0󸀠󸀠 (𝑥𝑛+V ) 󸀠󸀠 (𝐹0 (𝑥𝑛+1 ) ⋅⋅⋅ 𝐹4󸀠 (𝑥𝑛 ) ⋅⋅⋅ 𝐹4󸀠󸀠 (𝑥𝑛+V ) 𝐹4󸀠󸀠 (𝑥𝑛+1 )) ) ⋅ ⋅ ⋅ 𝐹4󸀠󸀠 (𝑥𝑛 ) ) ) ⋅⋅⋅ (6) and 𝐺𝑖 is obtained by replacing the 𝑖th column of 𝐺 by the vector 𝑉. Let the following conditions be satisfied: Π (𝑥𝑛 ) = 𝑦𝑛 , 1 Π󸀠 (𝑥𝑛 ) = 𝑦𝑛󸀠 , 󸀠 󸀠 = ℎ𝑦𝑛󸀠 + ℎ2 ( ∑ 𝛾𝑗󸀠 𝑓𝑛+𝑗 + 𝛾𝑛+V 𝑓𝑛+V ) , ℎ𝑦𝑛+V 𝑗=0 where 𝛽𝑗 , 𝛽𝑗 (...truncated)