**6. Numerical Examples**

The methods developed originally by Gautschi [1] and those that follow by Neta and Ford [2] fit low order monomials and the sine and cosine functions of multiples of the frequency. On the other hand, the methods developed later by Vanden Berghe and Van Daele [44] use monomials and product of monomials and sine and cosine functions of the frequency. We will demonstrate via the first three examples the difference between the two strategies. Vanden Berghe and Van Daele [50] compared the two approaches in some cases but not used schemes developed by Neta and Ford [2]. In the latter examples we also compare the results to Runge-Kutta-Nyström based method (46)–(49), see [49].

First, we list a method of trigonometric order 2 based on the idea of Vanden Berghe and Van Daele [44], which we obtained using Maple software, see Chun and Neta [51].

$$a\_1y\_{n+1} + a\_1y\_n + a\_2y\_{n-1} = h^2(b\_1f\_n + b\_2f\_{n-1} + b\_3f\_{n-2}),\tag{51}$$

where

$$\begin{aligned} a\_1 &= -\sin(v)v - 2\cos(v), \\ a\_2 &= -1 - a\_1, \\ b\_1 &= \frac{v(v\sin(v) - 1)(\cos(v) + 1) + 2\sin(v)}{v^3(1 + \cos(v))}, \\ b\_2 &= \frac{v(2 - v\sin(v))(\cos(v) + 1) - 4\sin(v)\cos(v)}{v^3(1 + \cos(v))}, \\ b\_3 &= \frac{2 - v\sin(v) - 2\cos(v)}{v^3\sin(v)}. \end{aligned} \tag{52}$$

The Taylor series expansion of the coefficients are

$$\begin{aligned} a\_1 &= -2 + \frac{1}{12}v^4 - \frac{1}{180}v^6, \\ a\_2 &= -1 - a\_1, \\ b\_1 &= \frac{13}{12} - \frac{19}{120}v^2 + \frac{37}{4032}v^4 - \frac{41}{362880}v^6, \\\\ b\_2 &= -\frac{1}{6} + \frac{3}{20}v^2 - \frac{59}{10080}v^4 + \frac{13}{36288}v^6, \\\\ b\_3 &= \frac{1}{12} + \frac{1}{120}v^2 + \frac{17}{20160}v^4 + \frac{31}{362880}v^6. \end{aligned} \tag{53}$$

**Example 1.** *The first example is chosen so that the exact solution in a combination of sine and cosine of multiples of the frequency, i.e.,*

$$y''(x) + 9y(x) = 3\sin(6x), \ 0 \le x \le 40\pi \tag{54}$$

*subject to the initial conditions*

$$y(0) = 1, \ y'(0) = 3.\tag{55}$$

*The exact solution is*

$$y\_{exact}(\mathbf{x}) = \frac{11}{9}\sin(3\mathbf{x}) + \cos(3\mathbf{x}) - \frac{1}{9}\sin(6\mathbf{x}).\tag{56}$$

The results using *h* = *π*/500 are given in Table 2. We expect the methods that fit sine and cosine of multiples of the frequency will do better.

**Table 2.** The *L*2 norm of the error for the first example using three methods for three different values around the exact frequency.


Based on the results we see that (25) is best when the frequency is known exactly. If it is not known exactly, the method prefers underestimation of the frequency. The second best is (30). This method will have no preference to underestimation.

**Example 2.** *The second example is very similar*

$$y''(x) + 9y(x) = 3\sin(3x), \ 0 \le x \le 40\pi \tag{57}$$

*subject to the initial conditions*

$$y(0) = 1, \ y'(0) = 3.\tag{58}$$

*The exact solution is*

$$y\_{\text{exact}}(\mathbf{x}) = \frac{7}{6}\sin(3\mathbf{x}) + \cos(3\mathbf{x}) - \frac{1}{2}\mathbf{x}\cos(3\mathbf{x}).\tag{59}$$

*The results are given in Table 3. Now we expect that the method (51) due to Chun and Neta [51] will perform better, since the exact solution has a product of monomial and cosine. In fact this is the case followed by (25). The method (30) requires smaller step size to converge and the results are not as good as those of the other two methods. Note that for this example all methods have no preference to underestimation of the frequency.*

**Table 3.** The *L*2 norm of the error for the second example using three methods for three different values around the exact frequency.


**Example 3.** *What if the frequency of the forcing term is not an integer multiple of the frequency of the homogeneous solution? We now consider the following example*

$$y''(x) + 9y(x) = 3\sin(4x), \ 0 \le x \le 40\pi \tag{60}$$

*subject to the initial conditions*

$$y(0) = 1, \ y'(0) = 3.\tag{61}$$

*The exact solution is*

$$y\_{exact}(\mathbf{x}) = \frac{11}{7}\sin(3\mathbf{x}) + \cos(3\mathbf{x}) - \frac{3}{7}\sin(4\mathbf{x}).\tag{62}$$

The results are given in Table 4. Now (51) is best followed by (30).


**Table 4.** The *L*2 norm of the error for the third example using three methods for three different values around the exact frequency.

Based on the three examples, we find that (51) is best in the last two examples, but not in the first case where the frequency of the forcing term is a multiple of the frequency of the homogeneous solution.

Before moving to the rest of the experiments, we have decided to rerun the first example on a much longer interval. This will test the quality of those methods in long-term integration. We have taken the same step size *h* = *π*/500 and integrated for 0 ≤ *x* ≤ 4000*π*. The results are given in Table 5. It is clear that the method due to Neta and Ford is no longer viable. The method (51) gave same errors when *ω* = 3 but all other cases show lower accuracy at the end of the longer interval

**Table 5.** The *L*2 norm of the error for the first example using three methods for three different values around the exact frequency.


**Example 4.** *The fourth example is a system of two second order initial value problems*

$$\begin{aligned} u''(\mathbf{x}) &= -\frac{u(\mathbf{x})}{(u^2(\mathbf{x}) + v^2(\mathbf{x}))^{3/2}}, \ 0 \le \mathbf{x} \le 12\pi\\ v''(\mathbf{x}) &= -\frac{v(\mathbf{x})}{(u^2(\mathbf{x}) + v^2(\mathbf{x}))^{3/2}}, \ 0 \le \mathbf{x} \le 12\pi \end{aligned} \tag{63}$$

*subject to the initial conditions*

$$\begin{aligned} u(0) &= 0, \ u'(0) = 1, \\ v(0) &= 1, \ v'(0) = 0. \end{aligned} \tag{64}$$

*The exact solution is given by*

$$u\_{\text{exact}}(\mathbf{x}) = \sin(\mathbf{x}), \ v\_{\text{exact}}(\mathbf{x}) = \cos(\mathbf{x}).\tag{65}$$

*We have converted this to a system of first order equations and solved it numerically using implicit Adams methods of trigonometric orders 2 and 3 (denoted AI2 and AI3, respectively) and generalized Milne-Simpson* *methods (GMS) of the same order, which are (32) and (66), respectively. We also included results from [51] and Runge-Kutta-Nyström method [49]. The results are given in Table 6. For Adams implicit, we have used the Taylor series coefficients as given in [1]. For GMS with q* = 2*, we used the coefficients as given in [2]. They did not give the coefficients for q* = 3 *but suggested to numerically solve the system for the coefficients. We were able now to get the coefficients*

$$y\_{n+5} = y\_{n+3} + h \left( b\_0 f\_n + b\_1 f\_{n+1} + b\_2 f\_{n+2} + b\_3 f\_{n+3} + b\_4 f\_{n+4} + b\_5 f\_{n+5} \right),\tag{66}$$

*where*

$$\begin{aligned} b\_{0} &= \frac{1}{6} \frac{\sin(v)}{d\_{1}},\\ b\_{1} &= -\frac{1}{3} \frac{\sin(v)(2\cos^{2}(v)-1)}{d\_{2}},\\ b\_{2} &= \frac{1}{3} \frac{\sin(v)\left(16\cos^{5}(v)+8\cos^{4}(v)-16\cos^{3}(v)-6\cos^{2}(v)+4\cos(v)+1\right)}{d\_{1}},\\ b\_{3} &= -\frac{1}{3} \frac{\sin(v)\left(8\cos^{3}(v)-2\cos(v)+1\right)\left(4\cos^{3}(v)-4\cos(v)-1\right)}{d\_{1}},\\ b\_{4} &= \frac{1}{6} \frac{\sin(v)\left(16\cos^{4}(v)+24\cos^{3}(v)+4\cos^{2}(v)-2\cos(v)+1\right)}{d\_{2}},\\ b\_{5} &= \frac{2}{3} \frac{\sin(v)\cos^{2}(v)(4\cos(v)+3)}{d\_{1}}.\end{aligned} \tag{67}$$

*where*

*and*

$$\begin{aligned} d\_1 &= v \cos(v) \left( 8 \cos^3(v) + 8 \cos^2(v) - 1 \right), \\\\ d\_2 &= v \cos(v) \left( 4 \cos^2(v) + 2 \cos(v) - 1 \right). \end{aligned}$$

**Table 6.** The *L*2 norm of the error for the fourth example using two implicit methods of trigonometric orders 2 and 3 and one explicit from [51] and one based on Runge-Kutta-Nyström.

