New Seven-Step Numerical Method for Direct Solution of Fourth Order Ordinary Differential Equations

1 Introduction

The general fourth order ordinary differential equations (ODEs) of the form

\[y^{(4)} = f(x, y, y', y'', y''') \quad y^{(s)}(x_0) = y_s, s = 0 (1)3, x \in [x_0, b]\] (1)

are considered in this paper. Eq. (1) can be solved by reducing it to its equivalent first order system as mentioned in [1-9]. However, this approach suffers some setbacks, such as evaluation of too many functions and heavy computation (see [10-12]).

Direct methods of solving Eq. (1) have been examined by several researchers [10,12-14]. They developed linear multistep methods using interpolation and collocation whereby the use of the approximated power series as a basis function was considered. Kayode [12] developed an efficient zero-stable numerical method with step number k=4 and 5 for fourth order initial value problems (IVPs), which was implemented in predictor-corrector mode. Furthermore, a five-step block method for solving fourth order ODEs directly is presented in [10]. In addition, in [13] and [14] six-step block methods are developed for solving fourth order ODEs using a multistep collocation approach whereby collocation points are selected at some grid points. These methods, however, have low accuracy.

In order to improve the accuracy of the existing methods, this article proposes a new block method for directly solving general fourth order IVPs of ODEs by increasing step number k.

2 Methodology

Let the approximate solution to Eq. (1) be a power series of the form

\[y(x) = \sum_{j=0}^{k+4} a_j x^j\] (2)

where k is the step number. The fourth derivative of Eq. (2) gives

\[y^{(4)}(x) = \sum_{j=4}^{k+4} j(j-1)(j-2)(j-3)a_j x^{j-4} = f(x, y, y', y'', y''')\] (3)

Eq. (2) is interpolated at \(x = x_{n+i}\), i = (k-5)(1)(k-2) and Eq. (3) is collocated at \(x = x_{n+1}, i = 0(1)k\). This gives a system of equations in the form

\[AX = B \tag{4}\] where

\[X = \begin{bmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{k+4} \end{bmatrix}, B = \begin{bmatrix} y_{n+k-5} \\ y_{n+k-4} \\ y_{n+k-3} \\ \vdots \\ y_{n+k} \end{bmatrix},\]

\[\text{[rumus tidak dapat ditampilkan dengan baik — lihat PDF asli]}\]

By using Gaussian elimination, the values of \(a_j\)'s, j = 0(1)k + 4 in Eq. (4) are obtained and then substituted into Eq. (2) to give a continuous linear multistep method in the form

\[y(t) = \sum_{i=k-5}^{k-2} \alpha_j(t) y_{n+j} + h^4 \sum_{i=0}^{k} \beta_j(t) f_{n+j}\] (5)

For k = 7, we have

\[\begin{pmatrix} \alpha_{2}(t) \\ \alpha_{3}(t) \\ \alpha_{4}(t) \\ \alpha_{5}(t) \end{pmatrix} = \begin{pmatrix} -1 & \frac{-11}{6} & -1 & \frac{-1}{6} \\ 4 & 7 & \frac{7}{2} & \frac{1}{2} \\ -6 & \frac{-19}{2} & -4 & \frac{-1}{2} \\ 4 & \frac{13}{3} & \frac{3}{2} & \frac{1}{6} \end{pmatrix} \begin{pmatrix} t^{0} \\ t^{1} \\ t^{2} \\ t^{3} \\ t^{4} \\ t^{5} \\ t^{6} \\ t^{7} \\ t^{8} \\ t^{9} \\ t^{10} \\ t^{11} \\ t^{2} \end{pmatrix}\] where

\[\text{[rumus tidak dapat ditampilkan dengan baik — lihat PDF asli]}\]

The values of w, u and v are w = 119750400, u = 13305600, v = 23950080, for \(t = \frac{x - x_{n+6}}{h} .\)

Eq. (5) is evaluated at the non-interpolating points, i.e. t = -6, -5, 0 and 1, to produce the discrete schemes. The first, second and third derivatives of Eq. (5) are evaluated at all the points within the interval, i.e. t = -6, -5, -4, -3, -2, -1, 0and 1, to give the derivatives of the discrete schemes. These schemes are combined in a matrix, whereby both the y and the f function are multiplied by the inverse of the coefficients of \(y_{n+j}\), j = 0(1)k. This yields a block of the form

\[A^{0}Y_{N} = A'Y_{N-1} + hA''Y_{N-1}' + h^{2}B'Y_{N-1}'' + h^{3}B''Y_{N-1}''' + h^{4}(E^{0}F_{N} + E^{1}F_{N-1})\] (6)

where

\[Y_{N} = \begin{pmatrix} y_{n+1} \\ y_{n+2} \\ \vdots \\ y_{n+k} \end{pmatrix}, Y_{N-1} = \begin{pmatrix} y_{n-k+1} \\ y_{n-k+2} \\ \vdots \\ y_{n} \end{pmatrix}, Y_{N-1}' = \begin{pmatrix} y_{n-k+1}' \\ y_{n-k+2}' \\ \vdots \\ y_{n}' \end{pmatrix}, Y_{N-1}'' = \begin{pmatrix} y_{n-k+1}'' \\ y_{n-k+2}'' \\ \vdots \\ y_{n}'' \end{pmatrix}, Y_{N-1}'' = \begin{pmatrix} y_{n-k+1}'' \\ y_{n-k+2}'' \\ \vdots \\ y_{n}'' \end{pmatrix}, Y_{N-1}'' = \begin{pmatrix} y_{n-k+1}'' \\ y_{n-k+2}' \\ \vdots \\ y_{n}'' \end{pmatrix}, Y_{N-1}'' = \begin{pmatrix} f_{n-k+1} \\ f_{n-k+2} \\ \vdots \\ f_{n} \end{pmatrix}.\]

If k = 7, we obtain

\[A^{0} = \begin{pmatrix} 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{pmatrix}, \qquad A' = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 \end{pmatrix},\]

\[A'' = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 2 \\ 0 & 0 & 0 & 0 & 0 & 0 & 3 \\ 0 & 0 & 0 & 0 & 0 & 0 & 4 \\ 0 & 0 & 0 & 0 & 0 & 0 & 6 \\ 0 & 0 & 0 & 0 & 0 & 0 & 7 \end{pmatrix}, \quad B' = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & \frac{1}{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 2 \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{9}{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{9}{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{25}{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{25}{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & \frac{49}{2} \end{pmatrix},\]

\[B'' = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & \frac{1}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{4}{3} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{9}{2} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{32}{3} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{32}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{125}{6} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{343}{6} \end{pmatrix}, \quad E^1 = \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & \frac{3102701}{119750400} \\ 0 & 0 & 0 & 0 & 0 & \frac{315461}{467775} \\ 0 & 0 & 0 & 0 & 0 & \frac{526077}{492800} \\ 0 & 0 & 0 & 0 & 0 & \frac{312608}{467775} \\ 0 & 0 & 0 & 0 & 0 & \frac{25842625}{4790016} \\ 0 & 0 & 0 & 0 & 0 & \frac{110052}{11550} \\ 0 & 0 & 0 & 0 & 0 & 0 & \frac{262892693}{17107200} \end{pmatrix}\]

The corresponding derivatives of Eq. (6) are given by

\[\begin{pmatrix} y'_{n+1} \\ y'_{n+2} \\ y'_{n+3} \\ y'_{n+4} \\ y'_{n+5} \\ y'_{n+6} \\ y'_{n+7} \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} y'_{n} + \begin{pmatrix} 1 \\ 2 \\ 3 \\ 4 \\ 5 \\ 6 \\ 7 \end{pmatrix} hy'' + \begin{pmatrix} \frac{1}{2} \\ \frac{9}{2} \\ 8 \\ \frac{25}{2} \\ 18 \\ \frac{49}{2} \end{pmatrix} h^{2} y''' + h^{3} T \begin{pmatrix} f_{n} \\ f_{n+1} \\ f_{n+2} \\ f_{n+3} \\ f_{n+4} \\ f_{n+5} \\ f_{n+6} \\ f_{n+7} \end{pmatrix}\] where

\[T = \begin{bmatrix} \frac{335799}{3628800} & \frac{562618}{362800} & \frac{-662757}{362800} & \frac{694230}{362800} & \frac{-506675}{362800} & \frac{239406}{362800} & \frac{-65823}{362800} & \frac{8002}{362800} \\ \frac{13376}{28350} & \frac{38762}{28350} & \frac{-32823}{28350} & \frac{34730}{28350} & \frac{-25300}{28350} & \frac{11934}{28350} & \frac{-3277}{28350} & \frac{398}{28350} \\ \frac{51327}{44800} & \frac{183654}{44800} & \frac{-105381}{44800} & \frac{135450}{44800} & \frac{-98955}{44800} & \frac{46818}{44800} & \frac{-12879}{44800} & \frac{1566}{44800} \\ \frac{29976}{14175} & \frac{118432}{14175} & \frac{-46608}{14175} & \frac{86880}{14175} & \frac{-7632}{14175} & \frac{928}{14175} \\ \frac{2451875}{145152} & \frac{2052250}{145152} & \frac{-577125}{145152} & \frac{1603750}{145152} & \frac{-891875}{145152} & \frac{456750}{145152} & \frac{-125375}{145152} & \frac{15250}{145152} \\ \frac{6912}{1400} & \frac{30024}{1400} & \frac{-6156}{1400} & \frac{24840}{1400} & \frac{-10800}{1400} & \frac{7128}{1400} & \frac{-1764}{1400} & \frac{216}{1400} \\ \frac{3520209}{518400} & \frac{15697738}{518400} & \frac{-2369787}{518400} & \frac{13613670}{518400} & \frac{-4621925}{518400} & \frac{4336206}{518400} & \frac{-655473}{518400} & \frac{114562}{518400} \\ \hline \end{cases}\]

\[\text{[rumus tidak dapat ditampilkan dengan baik — lihat PDF asli]}\] for

\[\begin{pmatrix} y_{n+1}^{"'} \\ y_{n+2}^{"'} \\ y_{n+3}^{"'} \\ y_{n+4}^{"'} \\ y_{n+5}^{"'} \\ y_{n+6}^{"'} \\ y_{n+7}^{"'} \end{pmatrix} = \begin{pmatrix} 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \\ 1 \end{pmatrix} y_{n}^{"'} + hN \begin{pmatrix} f_{n+1} \\ f_{n+2} \\ f_{n+3} \\ f_{n+4} \\ f_{n+5} \\ f_{n+6} \\ f_{n+7} \end{pmatrix}\] for

3 Analysis of the Properties of the Block Method

3.1 Order of the Method

The linear operator associated with Eq. (6) can be defined as

\[L\{y(x):h\} = A^{0}Y_{N} - A'Y_{N-1} - hA''Y_{N-1}' - h^{2}B'Y_{N-1}'' - h^{3}B''Y_{N-1}'''\]\[-h^{4}(E^{0}F_{N} + E'F_{N-1})\] (7)

Eq. (7) is expanded in Taylor series, which gives

\[L[y(x):h] = C_0 y(x) + C_1 h y'(x) + \dots + C_p h^{(p)} y^{(p)}(x) + C_{p+1} h^{(p+1)} y^{(p+1)}(x) + \dots\]

The Eq. (6) and the associated linear operator are said to have order p if \(C_0 = C_1 = C_2 = \dots = C_p = C_{p+1} = C_{p+2} = C_{p+3} \, 0, C_{p+4} \neq 0\). Therefore, our method's Eq. (6) has order (p) [8,8,8,8,8,8,8]^T with error constants

\[C_{p+4} = \left[\frac{-40}{90159}, \frac{-44}{6727}, \frac{-73}{2765}, \frac{-53}{781}, \frac{-454}{3269}, \frac{-332}{1343}, \frac{-317}{791}\right]^T \cdot C_{p+4}h^{p+4}y^{(p+4)}(x_n)\] is known as the principal local truncation error at point \(x\).

3.2 Zero Stability of the Method

Block Eq. (6) is said to be zero-stable if the roots \(z_s = 1, 2, ..., N\) of the first characteristic polynomial \(\rho(z) = \det(zA^0 - A')\) satisfies \(|z| \le 1\) and the root |z| = 1 has multiplicity not greater than the order of the differential equation, which is 4.

Now, \(\rho(z) = \det(zA^0 - A') = 0\) implies that \(\rho(z) = z^6(z-1)\). Hence z = 0,0,0,0,0,0,0,0,0,0,0,0. Therefore, the new method's Eq. (6) is convergent because it is zero-stable and has order greater than one.

4 Numerical Experiment

The accuracy of the new method is examined by solving the following differential problems.

Problem 1: \[y^{iv} = x\], \(y(0) = 0\), \(y'(0) = 1\), \(y''(0) = y'''(0) = 0\), \(h = 0.1\)

Exact solution: \[y(x) = \frac{x^5}{120} + x\]

Problem 2: \[y^{iv} + y'' = 0\], \(0 \le x \le \frac{\pi}{2}\), \(y(0) = 0\), \(y'(0) = \frac{-1.1}{72 - 50\pi}\),

\[y''(0) = \frac{1}{144 - 100\pi}, \ y'''(0) = \frac{1.2}{144 - 100\pi}, h = 0.01\]

Exact solution: \[y(x) = \frac{1 - x - \cos x - 1.2 \sin x}{144 - 100\pi}\]

Problem 3: \[y^{iv} = (x^4 + 14x^3 + 49x^2 + 32x - 12)e^x\], \(y(0) = y'(0) = 0\), \(y''(0) = 2\), \(y'''(0) = -6\), \(0 \le x \le 1\).

Exact Solution: \(y(x) = x^2 (1-x)^2 e^x\)

Problem 4: \[y^{iv} = y\] \(y(0) = y'(0) = y''(0) = y'''(0) = 1\), \(0 \le x \le 1\)

Exact Solution: \(y(x) = e^x\)

The results generated after solving the above problems are shown in Tables 1–4.

The following notations are used in Tables 3 and 4:

S2PEB: Sequential implementation of the 2-Point Explicit Block Method.

P2PEB: Parallel implementation of the 2-Point Explicit Block Method.

S3PEB: Sequential implementation of the 3-Point Explicit Block Method.

P3PEB: Parallel implementation of the 3-Point Explicit Block Method.

Table 1 Comparison of new method with [13] and [14] for solving problem 1.

Table 2 Comparison of new method with [10] and [12] for solving problem 2.

Table 3 Comparison of new method with [7] for solving problem 3.

Table 4 Comparison of new method with [7] for solving problem 4.

5 Conclusion

A seven-step block method for the solution of fourth order ODEs is proposed in this paper. The new method was used to solve fourth order IVPs of ODEs. The numerical results were compared with the existing methods. The new method performed better than the methods in [10,12-14], which employed 5 and 6 steps (refer to Tables 1 and 2). This implies that better accuracy can be achieved when step number k is increased. The accuracy of the new method was also found to be better than the method in [7], which was developed through numerical integration using 8 steps (refer to Tables 3 and 4). Thus, based on the numerical results, the new method outperformed the existing methods in terms of accuracy and should be considered as a viable alternative for directly solving fourth order initial value problems.