\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small Ninth MSU-UAB Conference on Differential Equations and Computational Simulations. \emph{Electronic Journal of Differential Equations}, Conference 20 (2013), pp. 119--132.\newline ISSN: 1072-6691. URL: http://ejde.math.txstate.edu or http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2013 Texas State University - San Marcos.} \vspace{9mm}} \begin{document}\setcounter{page}{119} \title[\hfilneg EJDE-2013/Conf/20/ \hfil Solving oscillatory problems] {Solving oscillatory problems using a block hybrid trigonometrically fitted method with two off-step points} \author[F. F. Ngwane, S. N. Jator \hfil EJDE-2013/conf/20 \hfilneg] {Fidele F. Ngwane, Samuel N. Jator} % in alphabetical order \address{Fidele F. Ngwane \newline Department of Mathematics, University of south Caroline, Salkehatchie, SC 29488, USA} \email{fifonge@yahoo.com} \address{Samuel N. Jator \newline Department of Mathematics and Statistics, Austin Peay State University, Clarksville, TN 37044, USA} \email{Jators@apsu.edu} \thanks{Published October 31, 2013.} \subjclass[2000]{65L05, 65L06} \keywords{First order system; hybrid method; trigonometrically fitted method; \hfill\break\indent off-step point; block method} \begin{abstract} A continuous hybrid method using trigonometric basis (CHMTB) with two `off-step' points is developed and used to produce three discrete hybrid methods which are simultaneously applied as numerical integrators by assembling them into a block hybrid trigonometrically fitted method (BHTFM) for solving oscillatory initial value problems (IVPs). The stability property of the BHTFM is discussed and the performance of the method is demonstrated on some numerical examples to show accuracy and efficiency advantages. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \newtheorem{definition}[theorem]{Definition} \newtheorem{example}[theorem]{Example} \newtheorem{remark}[theorem]{Remark} \allowdisplaybreaks \section{Introduction} In this article, we consider the first-order differential equation $$\label{eq-1} y' = f( x,y), \quad y(a) = y_{0}, \quad x \in [a,b],$$ with oscillating solutions where $f:\mathbb{R} \times \mathbb{R}^m\to \mathbb{R}^m$, $y, y_{0} \in \mathbb{R}^m$. Oscillatory IVPs frequently arise in areas such as classical mechanics, celestial mechanics, quantum mechanics, and biological sciences. Several numerical methods based on the use of polynomial basis functions have been developed for solving this class of important problems (see Lambert \cite{Lam1,Lam2}, Hairer et al in \cite{HWS}, Hairer \cite{HWS10}, and Sommeijer \cite{SJ}). Other methods based on exponential fitting techniques which take advantage of the special properties of the solution that may be known in advance have been proposed (see Simos \cite{SIMOS}, Vanden et al \cite{VAN}, Vigo-Aguiar et al \cite{RAMOS}, Franco \cite{FRANCO}, Fang et al \cite{FANG}, Nguyen et al \cite{CONG}, and Jator et al \cite{JSF}). The motivation governing the exponentially-fitted methods is inherent in the fact that if the frequency or a reasonable estimate of it is known in advance, these methods will be more advantageous than the polynomial based methods. The goal of this article is to construct a CHMTB which provides three discrete methods that are combined and applied as a BHTFM which takes the frequency of the solution as a priori knowledge. The coefficients of the BHTFM are functions of the frequency and the stepsize, hence the solutions provided by the proposed method are exact if \eqref{eq-1} has periodic solutions with known frequencies. We are motivated to use the hybrid method in order to increase the order of the method, while preserving good stability properties. Hybrid methods were first proposed to overcome the Dahlquist \cite{DQ} barrier theorem whereby the conventional linear multistep method was modified by incorporating off-step points in the derivation process (see Gear \cite{GE1}, Gragg et al in \cite{GS}, Butcher \cite{BU}, Gupta \cite{GUP}, Lambert \cite{Lam2}, and Kohfeld et al in \cite{KT1}). These methods were shown to enjoy both higher order and good stability properties, but included additional off-grid functions. Gupta \cite{GUP} noted that the design of algorithms for hybrid methods is more tedious due to the occurrence of off-step functions which increase the number of predictors needed to implement the methods. In order to avoid this deficiency, we adopt a different approach based on a block-by-block implementation instead of the traditional step-by-step implementation, generally performed in a predictor-corrector mode. Hence, we adopted the approach given in Jator et al in \cite{JSF}, where the CHMTB is used to generate a main and two additional methods which are combined and used as a BHTFM to simultaneously produce approximations $\{y_{n+\mu}, y_{n+v}, y_{n+1}\}\text{ at a block of points } \{x_{n+\mu}, x_{n+v}, x_{n+1}\},$ $h=x_{n+1}-x_n$, $n=0, \dots, N-1$, on a partition $[a,b]$, where $\mu, v \in (0,1)$, $a, b \in {\mathbb R}$, $h$ is the constant stepsize, $n$ is a grid index and $N>0$ is the number of steps. Block methods have also been considered by Shampine and Watts \cite{SHW}. We emphasize that the BHTFM simultaneously generates approximations $\{y_{n+\mu}, y_{n+v}, y_{n+1}\}$ to the exact solutions $\{y(x_{n+\mu}), y(x_{n+v}), y(x_{n+1})\}$. To apply the block method at the next block to obtain $y_{n+2}$, the only necessary starting value is $y_{n+1}$, and the loss of accuracy in $y_{n+1}$ does not affect subsequent points, thus the order of the algorithm is maintained. It is unnecessary to make a function evaluation at the initial part of the new block. Thus, at all blocks except the first, the first function evaluation is already available from the previous block. The organization of this article is as follows. In Section 2, we obtain a trigonometric basis representation $U(x)$ for the exact solution $y(x)$ which is used to generate two discrete methods which are combined into a BHTFM for solving \eqref{eq-1}. The analysis and implementation of the BHTFM are discussed in Section 3. Numerical examples are given in Section 4 to show the accuracy and efficiency of the BHTFM. Finally, we give some concluding remarks in Section 5. \section{Development of method} In this section, our objective is to construct a CHMTB which produces three discrete methods as by-products. The main method has the form $$\label{eq-2} y_{n+1}=y_n + h (\beta_{0}(u) f_n + \beta_{1}(u) f_{n+1} + \beta_{v}(u) f_{n+v} + \beta_{\mu}(u) f_{n+\mu}),$$ and the additional methods are given by $$\label{eq-3} \begin{gathered} y_{n+v}=y_n + h (\hat{\beta}_{0}(u) f_n + \hat{\beta}_{1}(u)f_{n+1} + \hat{\beta}_{v}(u) f_{n+v} + \hat{\beta}_{\mu}(u) f_{n+\mu}), \\ y_{n+\mu}=y_n + h (\check{\beta}_{0}(u) f_n + \check{\beta}_{1}(u)f_{n+1} + \check{\beta}_{v}(u) f_{n+v} + \check{\beta}_{\mu}(u) f_{n+\mu}), \end{gathered}$$ where $u = w h$, $\beta_{j}(u)$, $\hat{\beta}_{j}(u)$, $\check{\beta}_{j}(u)$, $\beta_{v}(u)$, $\hat{\beta}_{v}(u)$, $\check{\beta}_{v}(u)$, $\beta_{\mu}(u)$, $\hat{\beta}_{\mu}(u)$ and $\check{\beta}_{\mu}(u)$, $j=0,1$, are coefficients that depend on the stepsize and frequency. We note that $y_{n+v}$ and $y_{n+\mu}$ are the numerical approximation to the analytical solutions $y(x_{n+v})$, and $y(x_{n+\mu})$ respectively and $f_{n+v}=f(x_{n+v},y_{n+v})$, $f_{n+\mu}=f(x_{n+\mu},y_{n+\mu})$, $f_{n+j}=f(x_{n+j},y_{n+j})$ with $j = 0, 1$. To obtain equations \eqref{eq-2} and \eqref{eq-3}, we proceed by seeking to approximate the exact solution $y(x)$ on the interval $[x_n, x_n+h]$ by the interpolating function $U(x)$ of the form $$\label{eq-4} U(x)= a_0 +a_1 x + a_2 x^2 + a_3 \sin(w x) + a_4 \cos(w x),$$ where $a_0, a_1, a_2, a_3$ and $a_4$ are coefficients that must be uniquely determined. We then impose that the interpolating function equation \eqref{eq-4} coincides with the analytical solution at the end point $x_n$ to obtain the equation $$\label{eq-5} U(x_n)=y_n.$$ We also require that the function \eqref{eq-4} satisfy the differential equation \eqref{eq-1} at the points $x_{n+\mu}, x_{n+v}, x_{n+j}$, $j=0,1$ to obtain the following set of three equations: % $$\label{eq-6} U'(x_{n+\mu})=f_{n+\mu}, \quad U'(x_{n+v})=f_{n+v}, \quad U'(x_{n+j})=f_{n+j}, \quad j=0,1.$$ Equations \eqref{eq-5} and \eqref{eq-6} lead to a system of five equations which is solved by Cramer's rule to obtain $a_j$, $j = 0, 1, 2, 3, 4$. Our continuous CHMTB is constructed by substituting the values of $a_j$ into equation \eqref{eq-4}. After some algebraic manipulation, the CHMTB is expressed in the form $$\label{eq-7} U(x)=y_n + h(\beta_{0}(w, x) f_n + \beta_{1}(w, x) f_{n+1} + \beta_{v}(w, x) f_{n+v} + \beta_{\mu}(w, x) f_{n+\mu}),$$ where $w$ is the frequency, $\beta_0(w,x)$, $\beta_1(w,x)$, $\beta_v(w,x)$, and $\beta_\mu(w,x)$, are continuous coefficients. The continuous method \eqref{eq-7} is used to generate the main method of the form \eqref{eq-2} and two additional methods of the form \eqref{eq-3} by choosing, $v=\frac{1}{2}$ and $\mu=\frac{1}{4}$. Thus, evaluating \eqref{eq-7} at $x=\{x_{n+1}, x_{n+\frac{1}{2}}, x_{n+\frac{1}{4}}\}$ and letting $u = w h$, we obtain the coefficients of \eqref{eq-2} and \eqref{eq-3} as follows: $$\label{eq-8} \begin{gathered} \beta_0 = \frac{ \cos(\frac{u}{8}) (\csc(\frac{u}{4})^3) \sin(\frac{u}{8}) (u - 2 \sin(\frac{u}{2})) }{2u} \\ \beta_1 = -( \frac{ \cos(\frac{u}{8}) \csc(\frac{u}{4})^3 \sin(\frac{u}{8}) (-u + 2 \sin(\frac{u}{2}))}{2u})\\ \beta_v = -( \frac{\cos(\frac{u}{8}) \csc(\frac{u}{4})^3 \sin(\frac{u}{8}) (u \cos(\frac{u}{2}) - 2 \sin(\frac{u}{2}))}{u}), \end{gathered}$$ and $$\label{eq-9} \begin{gathered} \hat{\beta}_0 = \frac{ (\csc(\frac{u}{8})^2) (u - 4 \sin(\frac{u}{4})) }{8 u} \\ \hat{\beta}_v = \frac{(\csc(\frac{u}{8})^2) (u - 4 \sin(\frac{u}{4}))} {8 u} \\ \hat{\beta}_\mu = -(\frac{ (\csc(\frac{u}{8})^2) (u \cos(\frac{u}{4}) - 4 \sin(\frac{u}{4}))}{4 u}) \\ \check{\beta}_0 = \frac{(\csc(\frac{u}{4})^3) \sin(\frac{u}{8}) (8 u \cos(\frac{u}{8}) + 3 u \cos(\frac{(3 u)}{8}) - 8 (2 \sin(\frac{(3 u)}{8}) + \sin(\frac{(5 u)}{8})))} {16 u} \\ \check{\beta}_1 = -(\frac{(\csc(\frac{u}{4})^3) (u \cos(\frac{u}{8}) - 8 \sin(\frac{u}{8})) \sin(\frac{u}{8})} {16 u}) \\ \check{\beta}_v = \frac{(3 + 3 \cos(\frac{u}{4}) + \cos(\frac{u}{2})) ( \csc(\frac{u}{4})^3) (u \cos(\frac{u}{8}) - 8 \sin(\frac{u}{8})) \sin(\frac{u}{8})} {8 u} \\ \check{\beta}_\mu = -( \frac{(\cos(\frac{u}{8})^2) (\csc(\frac{u}{4})^3) \sin(\frac{u}{8}) (3 u \cos(\frac{u}{8}) + 3 u \cos(\frac{(3 u)}{8}) - 16 \sin(\frac{(3 u)}{8}))} { 4 u}), \end{gathered}$$ with $\hat{\beta}_1 = 0$ and $\beta_\mu = 0$. \section{Error analysis and stability} \subsection{Local truncation error} The Taylor series is used for small values of $u$ (see Simos \cite{SIMOS}). Thus the coefficients in equation \eqref{eq-8} can be expressed as \begin{gather*} %10 \beta_0 = \frac{1}{6} + \frac{u^2}{720} + \frac{u^4}{80640} + \frac{u^6}{9676800} + \frac{u^8}{1226244096} + \frac{691 u^{10}}{111588212736000}+ \dots \\ \beta_1 = \frac{1}{6} + \frac{u^2}{720} + \frac{u^4}{80640} + \frac{u^6}{9676800} + \frac{u^8}{1226244096} + \frac{691u^{10}}{111588212736000} + \dots \\ \beta_v = \frac{2}{3} - \frac{u^2}{360} - \frac{u^4}{40320} - \frac{u^6}{4838400} - \frac{u^8}{613122048} - \frac{691u^{10}}{55794106368000} + \dots \end{gather*} and the coefficients in equation \eqref{eq-9} can be expressed as \begin{gather*} %11 \begin{aligned} \hat{\beta}_0 &= \frac{1}{12} + \frac{u^2}{5760} + \frac{u^4}{2580480} + \frac{u^6}{1238630400} + \frac{u^8}{627836977152} \\ &\quad + \frac{691 u^{10}}{228532659683328000} + \dots \end{aligned}\\ \begin{aligned} \hat{\beta}_v &= \frac{1}{12} + \frac{u^2}{5760} + \frac{u^4}{2580480} + \frac{u^6}{1238630400} + \frac{u^8}{627836977152} \\ &\quad + \frac{691 u^{10}}{228532659683328000} + \dots \end{aligned}\\ \begin{aligned} \hat{\beta}_\mu &= \frac{1}{3} - \frac{u^2}{2880} - \frac{u^4}{1290240} - \frac{u^6}{619315200} - \frac{u^8}{313918488576} \\ &\quad - \frac{691 u^{10}}{114266329841664000} + \dots \end{aligned}\\ \begin{aligned} \check{\beta}_0 &= \frac{37}{384} + \frac{67u^2}{184320} + \frac{401u^4}{165150720} + \frac{1649u^6}{79272345600} + \frac{1711u^8}{9132174213120} \\ &\quad + \frac{12094087 u^{10}}{7313045109866496000} + \dots \end{aligned}\\ \begin{aligned} \check{\beta}_1 &= \frac{1}{384} + \frac{13 u^2}{184320} + \frac{37 u^4}{33030144} + \frac{1091 u^6}{79272345600} + \frac{2087 u^8}{14350559477760} \\ &\quad + \frac{10191073 u^{10}}{7313045109866496000} + \dots \end{aligned}\\ \begin{aligned} \check{\beta}_v & = -\frac{7}{192} + \frac{7 u^2}{46080} - \frac{11 u^4}{11796480} - \frac{29 u^6}{1415577600} - \frac{12503 u^8}{50226958172160}\\ &\quad - \frac{659969 u^{10}}{261180182495232000} + \dots \end{aligned}\\ \begin{aligned} \check{\beta}_\mu & = \frac{3}{16} - \frac{3 u^2}{5120} - \frac{3 u^4}{1146880} - \frac{31 u^6}{2202009600} - \frac{13 u^8}{155021475840} \\ &\quad - \frac{11747 u^{10}}{22571126882304000} + \dots. \end{aligned} \end{gather*} Thus, the Local Truncation Errors (LTEs) for methods \eqref{eq-2} and \eqref{eq-3} are given by $$%12 \begin{gathered} \text{LTE}\eqref{eq-2} = -\frac{h^5}{2880}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)), \\ \text{LTE}\eqref{eq-3} = -\frac{h^5}{92160}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)), \\ \text{LTE}\eqref{eq-3} = -\frac{53h^5}{1474560}(w^2y^{(3)}(x_n) + y^{(5)}(x_n)). \end{gathered}$$ \begin{remark} \label{rmk3.1} \rm We noticed that as $u \to 0$, the method (2) reduces to the fourth-order method of Gragg and Stetter \cite{GS}, which is based on polynomial basis functions. \end{remark} \subsection{Stability} The BHTFM is constructed by combining equations \eqref{eq-2} and \eqref{eq-3}, where the coefficients of the method are explicitly given by equations \eqref{eq-8} and \eqref{eq-9}. We then define the block-by-block method for computing vectors $Y_0, Y_1, \cdots$ in sequence (see \cite{FT}). Let the $\eta-$vector ($\eta = 3$ is the number of points within the block) $Y_\gamma, Y_{\gamma -1}, F_\gamma, F_{\gamma - 1}$ be given as \begin{gather*} Y_\gamma = (y_{n+ \frac{1}{4}}, y_{n+ \frac{1}{2}}, y_{n+1})^T, \quad Y_{\gamma-1} = (y_{n - \frac{1}{4}}, y_{n - \frac{1}{2}}, y_n)^T, \\ F_\gamma = ( f_{n+ \frac{1}{4}}, f_{n+ \frac{1}{2}}, f_{n+1})^T, \quad F_{\gamma-1} = (f_{n - \frac{1}{4}}, f_{n - \frac{1}{2}}, f_n)^T, \end{gather*} then the $1$-block $3$-point method for equation \eqref{eq-1} is given by $$\label{eq-13} Y_\gamma = \sum_{i=1}^1 A^{(i)}Y_{\gamma-i} + \sum_{i=0}^1 B^{(i)}F_{\gamma-i},$$ where $A^{(i)}, B^{(i)}$, $i = 0, 1$ are $3 \times 3$ matrices whose entries are given by the coefficients of \eqref{eq-2} and \eqref{eq-3}. \subsection*{Zero-stability} \begin{definition}\label{df-1} \rm The block method \eqref{eq-13} is zero stable provided the roots $R_j$, $j=1, 2, 3$ of the first characteristic polynomial $\rho(R)$ specified by $$\label{eq-14} \rho(R) =\det \Big[ \sum_{i=0}^1 A^{(i)}R^{1-i} \Big] = 0, \quad A^{(0)} = -I,$$ satisfies $|R_j| \leq 1$, $j= 1, 2, 3$ and for those roots with $|R_j| =1$, the multiplicity does not exceed $1$ (see \cite{FT}). \end{definition} \subsection*{Consistency of BHTFM} We note that the block method \eqref{eq-13} is consistent as it has order $p>1$. We see from equation \eqref{eq-14} and definition \eqref{df-1} that the block method \eqref{eq-13} is zero-stable since $\rho(R)=R^2(R-1)=0$ satisfies $|R_j| \leq 1$, $j=1, 2, 3$, and for those roots with $|R_j| = 1$, the multiplicity does not exceed $1$. The block method \eqref{eq-13} is thus convergent, as zero-stability plus consistency equals convergence. \subsection*{Linear stability of the BHTFM} To analyze the linear stability of the BHTFM, we apply the method to the test equation $y' =\lambda y$, where $\lambda$ is expected to run through the (negative) eigenvalues of the Jacobian matrix. Then an application of \eqref{eq-13} to the test equation yields $$\label{eq-15} Y_\gamma = M(q;u)Y_{\gamma-i}, \quad q = h \lambda, \quad u = w h,$$ where $M(q;u) = (A^{(0)}- qB^{(0)})^{-1} (A^{(1)}+ qB^{(1)}),$ where the matrix $M(q;u)$ is the amplification matrix which determines the stability of the method. \begin{definition}\label{df-2} \rm A region of stability is a region in the $q-u$ plane, throughout which $|\rho(q;u)| \leq 1$, where $\rho(q;u)$ is the spectral radius of $M(q;u)$. \end{definition} It is easily seen that the eigenvalues of $M(q;u)$ are $\lambda_1 =0$, $\lambda_2 =0$, and \begin{align*} \lambda_3 &= \Big(-16 (2 + q) (q^2 + u^2) \cos(\frac{u}{8}) + (6q^3 + 16u^2 + 6qu^2 + q^2(16 + u^2)) \cos(\frac{3u}{8}) \\ &\quad + 16 q^2 \cos(\frac{5u}{8}) + 10q^3 \cos(\frac{5u}{8}) + 16u^2 \cos(\frac{5u}{8}) + 10qu^2 \cos(\frac{5u}{8}) \\ &\quad + 3q^2u^2\cos(\frac{5u}{8}) + q^3u \sin(\frac{3u}{8}) + 3q^3u \sin(\frac{5 u}{8})\Big)\\ &\quad\div \Big(16 (-2 + q) (q^2 + u^2) \cos(\frac{u}{8}) + (-10 q^3 + 16 u^2 - 10 q u^2 \\ &\quad + q^2 (16 + 3 u^2)) \cos(\frac{3 u}{8}) + 16 q^2 \cos(\frac{5 u}{8}) - 6 q^3 \cos(\frac{5 u}{8}) + 16 u^2 \cos(\frac{5 u}{8}) \\ &\quad - 6 q u^2 \cos(\frac{5 u}{8}) + q^2 u^2 \cos(\frac{5 u}{8}) - 3 q^3 u \sin(\frac{3 u}{8}) - q^3 u \sin(\frac{5 u}{8})\Big). \end{align*} \begin{figure}[ht] \begin{center} \includegraphics[width=0.4\textwidth]{fig1} % StaRegP3 \end{center} \caption{The shaded region represents the truncated region of absolute stability} \label{StaReg} \end{figure} % \begin{figure}[ht] \begin{center} \includegraphics[width=0.4\textwidth]{fig2} % ZeroPolesP3 \end{center} \caption{$\lambda_{3}$ has three zeros$(\square)$ and no poles$(+)$ in $\mathbb{C}^-$, hence the BHTFM is $A_0$-stable.} \label{ZerosPoles} \end{figure} We observed that in the $q-u$ plane the BHTFM is stable for $q \in [-100, 100]$, and $u \in [-\pi, \pi ]$. Figure \ref{StaReg} is a plot of the stability region. Figure \ref{ZerosPoles} shows the zeros and poles of $\lambda_3$. \subsection{Implementation} We emphasize that methods \eqref{eq-2} and \eqref{eq-3} are combined to form the block method \eqref{eq-13}, which is implemented to solve \eqref{eq-1} without requiring starting values and predictors. For instance, if we let $n=0$ and $\gamma = 1$ in \eqref{eq-13}, then $(y_{\frac{1}{4}}, y_{\frac{1}{2}}, y_{1})^{T}$ are simultaneously obtained over the sub-interval $[x_{0},x_{1}]$, as $y_{0}$ is known from the IVP. Similarly, if $n = 1$ and $\gamma = 2$, then $(y_{\frac{5}{4}}, y_{\frac{3}{2}}, y_{2} )^{T}$ are simultaneously obtained over the sub-interval $[x_{1},x_{2}]$, as $y_{1}$ is known from the previous block, and so on; until we reach the final sub-interval $[x_{N-1},x_{N}]$. We note that for linear problems, we solve \eqref{eq-1} directly using the feature solve[~~] in Matlab, while nonlinear problems were solved by implementing the Newton's method in Matlab enhanced by the feature fsolve[~~]. \section{Numerical examples} In this section, we give numerical examples to illustrate the accuracy (small errors) and efficiency (fewer number of function evaluations (FNCs)) of the BHTFM. We find the approximate solution on the partition $\pi_{N}$, where \label{pi} \pi_{N}:a=x_{0}