Solving Nonlinear Boundary Value Problems Using the Higher Order Haar Wavelet Method

Abstract

The current study is focused on development and adaption of the higher order Haar wavelet method for solving nonlinear ordinary differential equations. The proposed approach is implemented on two sample problems—the Riccati and the Liénard equations. The convergence and accuracy of the proposed higher order Haar wavelet method are compared with the widely used Haar wavelet method. The comparison of numerical results with exact solutions is performed. The complexity issues of the higher order Haar wavelet method are discussed.

1. Introduction

Recently, the higher order Haar wavelet method (HOHWM) has been developed by the workgroup of the current study [1]. The HOHWM can be considered a refinement of the widely used HWM introduced in [2]. The HOHWM has been found to increase accuracy as well as convergence over the regular HWM [1,3,4,5,6].

The HWM was proposed originally for solving differential equations [2,7,8]. It has been shown in [9,10,11,12,13,14,15,16,17] that the HWM can also be successfully utilized for solving integral and integro-differential and integral equations. In pioneering works of Lepik [14,15,18,19,20,21], the integration techniques for the HWM were developed. The monograph [22] gives a thorough overview of the HWM and its application in different research areas. The weak formulation based HWM was introduced by Majak et al. in [23].

Engineering seems to be one of the most common areas of application for the HWM. This includes solid mechanics [14,20,21], composites [24,25,26,27,28,29,30,31,32,33], etc. In [28], free vibration analysis of the multilayer composite plate is performed and in [27,30,31,32,33] the HWM is successfully utilized for vibration analysis of plate and shell structures. The propagation of a solitary wave in a composite material is studied in [29]. Besides engineering, wavelet methods are used in a wide range of applications in areas of natural sciences, medicine, social and live sciences [2,34,35,36].

Pioneering study in the area of fractional differential and integro-differential equations was performed by Lepik in [15], where the Caputo derivatives are used to convert the fractional differential equations into integral equations, which only include integer order derivatives. Thus, the wavelet expansion is applied for integer order derivatives. An alternate approach was employed in [37,38,39]. According to the latter approach, the Caputo derivatives are expanded into Haar wavelets by employing a fractional operational matrix.

Recent developments of the HWM cover extension to the multidimensional case [40,41,42], treatment of adaptive methods [43], nonuniform Haar wavelet method [44], inverse problems [45], identification of software piracy [46], solutions to variational problems [47] and boundary layer flow problems [48], etc.

The accuracy and convergence are considered the most critical characteristics of the HWM and any numerical method. These aspects are studied in [49,50,51]. In [50], the convergence theorem is proven and the order of the convergence of the Chen and Hsiao approach based HWM has been found to be equal to two. Obviously, the HWM needs principal refinement in order to compete with similar strong formulation based method like differential quadrature method, used widely in engineering design. These results were confirmed in comparative study, performed by the workgroup, where HWM was compared with DQM and finite difference method [52]. Motivated by the latter conclusions Majak et al. developed the HOHWM in [1].

In the current study, the HOHWM is implemented for solving nonlinear differential equations. The Riccati equations are often selected for evaluation of numerical methods [53,54,55,56,57], especially for new ones [2,58,59,60], since the analytical solutions are known for a number of particular Riccati equations. The Riccati equations cover various actual research problems including Kalman filtering, model reduction optimal Control, etc. [59,61,62]. Liénard Equation [63,64,65] is chosen as an example of a nonlinear second order equation. Liénard equations arise when studying nonlinear oscillations [66,67]. A detailed analysis of the HOHWM is performed covering accuracy, convergence and complexity.

The structure of the article is following. In Section 2, the Haar functions and operational matrix of integration are given. In Section 3, the two model equations are introduced. In Section 4, the both Haar wavelet method approaches (HWM and HOHWM) for solving Riccati and Liénard problems are described. In Section 5, the numerical results are presented and analysed. Finally, the conclusions are given in Section 6.

2. Haar Wavelet Family

In the current study, the notation introduced in [21] is utilized. The integration domain $[A,B]$ can be divided into $2M$ equal subinteravals of length $\Delta x=(B-A)/\left(2M\right)$. The maximal level of resolution J is defined as $M=2^J$. The Haar functions are defined as

$$h_i\left(x\right)=\left\{\begin{array}{cc}\hfill 1& \mathrm{for}x\in [{\xi }_1\left(i\right),{\xi }_2\left(i\right)),\hfill \\ \hfill -1& \mathrm{for}x\in [{\xi }_2\left(i\right),{\xi }_3\left(i\right)),\hfill \\ \hfill 0& \mathrm{elsewhere},\hfill \end{array}\right.$$

(1)

where

$$\begin{array}{cc}\hfill {\xi }_1\left(i\right)& =A+2k\mu \Delta x,{\xi }_2\left(i\right)=A+(2k+1)\mu \Delta x,\hfill \\ \hfill {\xi }_3\left(i\right)& =A+2(k+1)\mu \Delta x,\mu =M/m,\Delta x=(B-A)/\left(2M\right).\hfill \\ \hfill j& =0,1,\dots ,J,k=0,1,\dots ,m-1,i=m+k1\hfill \end{array}$$

(2)

Index i is calculated from $i=m+k+1$. Parameter $m=2^j$ corresponds to the maximum number of square waves that can be sequentially deployed in interval $[A,B]$ for the given j and the parameter k indicates the location of the particular square wave. The scaling function $h_1\left(x\right)$ forms a special case in which $m=0,{\xi }_1=A,{\xi }_2=B,{\xi }_3=B$. The Haar series is orthogonal and therefore forms a good transform basis

$${\int }_A^B{h}_i\left(x\right)h_l\left(x\right)dx=\left\{\begin{array}{cc}\hfill 2^{-j}& i=l=2^j+k\hfill \\ \hfill 0& i\ne l.\hfill \end{array}\right.$$

(3)

Thus, any square integrable function $f\left(x\right)$ can be expanded into Haar wavelets as

$$f\left(x\right)=\sum _{i=1}^{\infty }{a}_i{h}_i\left(x\right),$$

(4)

where $a_i$ are the Haar coefficients.

Integrating Haar functions (1) one obtains

$$p_{n,i}\left(x\right)=\left\{\begin{array}{cc}\hfill 0& \mathrm{for}x\in [A,{\xi }_1\left(i\right)),\hfill \\ \hfill \frac{{(x-{\xi }_1\left(i\right))}^n}{n!}& \mathrm{for}x\in [{\xi }_1\left(i\right),{\xi }_2\left(i\right)),\hfill \\ \hfill \frac{{(x-{\xi }_1\left(i\right))}^n-2{(x-{\xi }_2\left(i\right))}^n}{n!}& \mathrm{for}x\in [{\xi }_2\left(i\right),{\xi }_3\left(i\right)),\hfill \\ \hfill \frac{{(x-{\xi }_1\left(i\right))}^n-2{(x-{\xi }_2\left(i\right))}^n+{(x-{\xi }_3\left(i\right))}^n}{n!}& \mathrm{for}x\in [{\xi }_3\left(i\right),B).\hfill \end{array}\right.$$

(5)

The Haar matrix can be written in terms of Haar functions as

$$H_{il}=h_i\left(x_l\right),$$

(6)

where $x_l=(l-1/2)\Delta x$ stand for collocation points. The operational matrix of integration $P_n$ is given as

$${\left(P_n\right)}_{il}=p_{n,i}\left(x_l\right).$$

(7)

Combining (6) and (7), the function approximation can be expressed in matrix notation as

$$f\left(x\right)=aH.$$

(8)

Based on (6) and (7), the n-th order operational matrix of integration can be written as

$$\underset{n}{\underbrace{{\int \dots \int }_A^x}}H\underset{n}{\underbrace{d\xi \cdots d\xi }}=P_n.$$

(9)

Obviously, the matrices H and $P_n$ depend on x.

In boundary points the relation (5) simplifies to

$$P_n\left(A\right)=0\forall n>0$$

(10)

and

$${\left(P_n\left(B\right)\right)}_i=p_{n,i}\left(B\right)=\frac{{(B-{\xi }_1\left(i\right))}^n-2{(B-{\xi }_2\left(i\right))}^n+{(B-{\xi }_3\left(i\right))}^n}{n!}.$$

(11)

In general, the relation (10) leads to simplified boundary conditions.

3. Model Equations

3.1. First Order Problems

Any nonlinear ODE of the form

$$u^{\prime }\left(x\right)+f\left(x\right)u^2\left(x\right)+g\left(x\right)u\left(x\right)+h\left(x\right)=0$$

(12)

can be called a Riccati equation. The boundary condition for Equation (12) can be given as

$$u\left(x_c\right)=u_c,$$

(13)

where $x_c$ is the boundary (A or B) and $u_c$ is a given constants. In the following the three particular sets of coefficients ($f\left(x\right)$, $g\left(x\right)$ and $h\left(x\right)$) of the Riccati equation are considered for which the exact solution is known. The domain $x\in [0,5]$ was chosen since it reflects the qualitative behaviour of the exact solution.

Problem 1 (constant coefficients):

$$\begin{array}{cc}\hfill u^{\prime }\left(x\right)& =1-u^2\left(x\right),u\left(0\right)=0,x\in [0,5],\hfill \\ \hfill u_e\left(x\right)& =\frac{e^{2x}-1}{e^{2x}+1},\hfill \end{array}$$

(14)

where $u_e\left(x\right)$ is the known exact solution.

Problem 2 (constant coefficients):

$$\begin{array}{cc}\hfill u^{\prime }\left(x\right)& =1+2u\left(x\right)-u^2\left(x\right),u\left(0\right)=0,x\in [0,5],\hfill \\ \hfill u_e\left(x\right)& =\frac{\left(-1+\sqrt{2}\right)\left(-1+e^{2\sqrt{2}x}\right)}{1+\left(3-2\sqrt{2}\right)e^{2\sqrt{2}x}}.\hfill \end{array}$$

(15)

Problem 3 (variable coefficients):

$$\begin{array}{cc}\hfill u^{\prime }\left(x\right)& =2x-u\left(x\right)/x+u^2\left(x\right)/x^3,u\left(5\right)=40,x\in [0,5],\hfill \\ \hfill u_e\left(x\right)& =\frac{2x^3+30x^2}{2x+15}.\hfill \end{array}$$

(16)

In order to use the HWM and HOHWM, the Riccati Equations (14)–(16) can be quasilinearized. The approach shown in [68] is used. By applying quasilinearization one obtains:

Problem 1

$$u_{r+1}^{\prime }\left(x\right)+2u_{r+1}\left(x\right)u_r\left(x\right)=u_r^2\left(x\right)+1,u_{r+1}\left(0\right)=0$$

(17)

Problem 2

$$u_{r+1}^{\prime }\left(x\right)+2u_{r+1}\left(x\right)\left(u_r\left(x\right)-1\right)=u_r^2\left(x\right)+1,u_{r+1}\left(0\right)=0,$$

(18)

Problem 3

$$u_{r+1}^{\prime }\left(x\right)+u_{r+1}\left(x\right)/x-2u_{r+1}\left(x\right)u_r\left(x\right)/x^3=2x-u_r^2\left(x\right)/x^3,u_{r+1}\left(5\right)=40.$$

(19)

In (17)–(19) r denotes the iteration step. The vector $u_{r+1}\left(x\right)$ has to be computed using the HWM or the HOHWM for each iteration until the difference between two consecutive solutions becomes sufficiently small (${10}^{-8}$ was used here). Clearly, an initial value for $r=0$ must be given to start the iteration. In the present study, the linear relation $u_0\left(x\right)=x$ was used. The same applies to the upcoming subsection as well.

3.2. Second Order Problems

A nonlinear differential equation of the form

$$u^{″}\left(x\right)+f\left(u\right)u^{\prime }\left(x\right)+g\left(u\right)=0$$

(20)

is known as a Liénard equation. In the present study, $f\left(u\right)=u$ and $g\left(u\right)=0$ are considered. In such a case, the exact solution is known. The two boundary conditions can be defined as

$$u\left(A\right)=c_A,u\left(B\right)=c_B,$$

(21)

where A and B are on the boundary of the domain, $c_A$ and $c_B$ are arbitrary constants.

Problem 4

$$\begin{array}{cc}\hfill u^{″}\left(x\right)& +u^{\prime }\left(x\right)u\left(x\right)=0,u\left(0\right)=0,u\left(1\right)=2\tan \mathrm{h}\left(1\right),x\in [0,1],\hfill \\ \hfill u_e\left(x\right)& =2\tan \mathrm{h}\left(x\right).\hfill \end{array}$$

(22)

Problem 5

$$\begin{array}{cc}\hfill u^{″}\left(x\right)& +u^{\prime }\left(x\right)u\left(x\right)=0,u\left(0\right)=0,u\left(1\right)=4\tan \mathrm{h}\left(2\right),x\in [0,1],\hfill \\ \hfill u_e\left(x\right)& =4\tan \mathrm{h}\left(2x\right).\hfill \end{array}$$

(23)

The model equations can again be quasiliearised. This results in the following problems:

Problem 4

$$\begin{array}{c}\hfill u_{r+1}^{″}\left(x\right)+u_{r+1}^{\prime }\left(x\right)u_r\left(x\right)+u_r^{\prime }{u}_{r+1}\left(x\right)=u_r^{\prime }\left(x\right)u_r\left(x\right),\\ \hfill u_{r+1}\left(0\right)=0,u_{r+1}\left(1\right)=2\tan \mathrm{h}\left(1\right).\end{array}$$

(24)

Problem 5

$$\begin{array}{c}\hfill u_{r+1}^{″}\left(x\right)+u_{r+1}^{\prime }\left(x\right)u_r\left(x\right)+u_r^{\prime }{u}_{r+1}\left(x\right)=u_r^{\prime }\left(x\right)u_r\left(x\right),\\ \hfill u_{r+1}\left(0\right)=0,u_{r+1}\left(1\right)=4\tan \mathrm{h}\left(2\right).\end{array}$$

(25)

4. Haar Wavelet Methods

The well known HWM involves expanding the highest order derivative present in the equation into the Haar wavelet series as

$$u_{r+1}^{\prime }\left(x\right)=a_{r+1}H,$$

(26)

where H is the Haar wavelet matrix described in (6) and $a_{r+1}$ is the row vector of Haar wavelet coefficients for the given iteration. By integrating the above, one can arrive at

$$u_{r+1}\left(x\right)=a_{r+1}{P}_1+c_1,$$

(27)

where $P_1$ is the first integration matrix introduced in (7) and $c_1$ is an unknown coefficient. The boundary condition (for other boundary value problems the process is analogous) $u\left(0\right)=0$ can be expressed as

$$u_{r+1}\left(0\right)=a_{r+1}{P}_1\left(0\right)+c_1=0.$$

(28)

Obviously, (10) and (28) imply that $c_1=0$. Substituting the obtained result into (27) yields

$$u_{r+1}\left(x\right)=a_{r+1}{P}_1.$$

(29)

The process is similar in case of a second order equation, in which case, two coefficients are obtained and the two boundary conditions are used to obtain said coefficients.

In the case of the HOHWM, one expands a higher order derivative into the Haar wavelet series. In general, according to [1], $2s$ (where $s>0$) extra derivatives are used.

$$u_{r+1}^{(1+2s)}\left(x\right)=a_{r+1}H.$$

(30)

In the simplest case where $s=1$ one obtains

$$u_{r+1}^{‴}\left(x\right)=a_{r+1}H.$$

(31)

Integration of Equation (31) results in expression

$$u_{r+1}\left(x\right)=a_{r+1}{P}_3+c_1{x}^2/2+c_{2}x+c_3,$$

(32)

which includes two extra unknown coefficients $c_2$ and $c_3$. The algorithms for determining integration constants for the HOHWM are described in details in [1] and are omitted herein. In the simplest case, where $s=1$, according to these algorithms, the differential equation should be fulfilled at boundary points.

Thus, utilizing the boundary condition $u\left(0\right)=0$ as well as the conditions $DE\left(0\right)=0$ and $DE\left(5\right)=0$ for (17), where $DE$ refers to the quasilinearized differential equation, one can arrive at

$$u_{r+1}\left(x\right)=a_{r+1}{P}_3+x+\frac{x^2\left(-a_{r+1}{P}_2\left(5\right)-2u_r\left(5\right)(a_{r+1}{P}_3\left(5\right)+5)+u_r{\left(5\right)}^2\right)}{50u_r\left(5\right)+10}.$$

(33)

The process is similar to the second order equation and is thus omitted. It must be noted, however, that a total of 4 coefficients are present in this case.

Either (29) or (33) can be substituted into differential Equation (17) depending on which wavelet method is implemented (HWM or HOHWM). The obtained linear algebraic system of equations should be solved with respect to the coefficient vector $a_{r+1}$. Finally, for known vector $a_{r+1}$ the solution of the differential equation can be computed by use of (29) for the HWM and (33) for the HOHWM. These steps must be carried out at each iterations step.

The expressions of the solutions for the HWM and the HOHWM, corresponding to problems 2–5 can be derived similarly to expressions (29) and (33). It can be shown that for problem 2, one obtains (29) in case of the HWM and

$$\begin{array}{cc}\hfill u_{r+1}\left(x\right)=& a_{r+1}{P}_3+x\hfill \\ & +\frac{x^2\left(10+2a_{r+1}{P}_3\left(5\right)-a_{r+1}{P}_2\left(5\right)-2u_r\left(5\right)(a_{r+1}{P}_3\left(5\right)+5)+u_r{\left(5\right)}^2\right)}{50u_r\left(5\right)+10}\hfill \end{array}$$

(34)

in case of the HOHWM. In case of problem 3, one can obtain

$$u=a_{r+1}{P}_1-a_{r+1}{P}_1\left(5\right)+40$$

(35)

and

$$u=a_{r+1}{P}_3+a_{r+1}\left(\left(\frac{x^2}{25}-\frac{2x}{5}\right)P_3\left(5\right)+\left(x-\frac{x^2}{5}\right)P_2\left(5\right)\right)$$

(36)

for the HWM and the HOHWM, respectively. When dealing with problems 4 and 5, it can be shown that for the HWM

$$u=a_{r+1}{P}_2-a_{r+1}x{P}_2\left(1\right)+2x{b}_r$$

(37)

can be obtained and for the HOHWM

$$\begin{array}{cc}\hfill u=& a_{r+1}{P}_4+a_{r+1}\left(\left(b_r\frac{x^3-3x}{6+2b_r}-\frac{3x}{3+b_r}\right)P_4\left(1\right)-b_r\frac{x^3-x}{6+2b_r}{P}_3\left(1\right)-\frac{x^3-x}{6+2b_r}{P}_2\left(1\right)\right)\hfill \\ \hfill & +\frac{3b_{r}x}{3+b_r}-b_r^2\frac{x^3-3x}{6+2b_r}\hfill \end{array}$$

(38)

can be obtained. In (37) and (38), $b_r$ corresponds to the right hand side boundary condition. In case of problem 4, $b_r=2\tan \mathrm{h}\left(1\right)$ and in case of problem 5, $b_r=4\tan \mathrm{h}\left(2\right)$.

Comparing the expressions (29), (35), (37) to (33) and (34), one can observe some principal differences between the HWM and the HOHWM. In case of the HOHWM, the solution depends strictly on the boundary values of the previous iteration, which is not the case with the HWM. However, in case of (36) and (38), such a dependence is not realized for the HOHWM.

5. Analysis of Numerical Results

The numerical experiments were carried out for three particular Riccati Equations, (14)–(16), as well as two particular Liénard Equations, (22) and (23). The numerical solutions for problem 1 are computed according to (29) or (33) and validated against exact solutions at the midpoint ($x=5/2$ for problems 1–3; $x=1/2$ for probelms 4 and 5). The order of convergence of the solution can be computed with ($k_J^e$) and without ($k_J$) use of exact solution as [1,49,50]

$$k_J^e={log}_2\left(\frac{F_{J-1}-F_e}{F_J-F_e}\right),k_J={log}_2\left(\frac{F_{J-2}-F_{J-1}}{F_{J-1}-F_J}\right),$$

(39)

where $F_J$ and $F_e$ stand for the numerical solution at resolution J and exact solution, respectively. In the present study the numerical solution at the midpoint $u(5/2)$ is used. Since $k_J^e$ depends on the solution at the current resolution and the previous resolution (as well as the exact solution), it cannot be calculated for the lowest resolution. Similarly, $k_J$ depends on the solution at the current resolution as well as two lower resolutions so it cannot be calculated for the two lowest resolutions. In the following $r_f$ refers to the number of iterations that were carried out in order to get a satisfactory result.

The numerical solution of Riccati equation for problem 1 as well as the absolute error $\Delta u(5/2)=|u_e(5/2)-u(5/2)|$, the rates of convergence and the number of iterations used to obtain the solution by applying HWM and HOHWM are given in

Table 1. Solutions to (17) where $u_e(5/2)=0.986614298151$ using HWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	1.0041392106	$1.75\times {10}^{-2}$	-	-	6
2	8	0.9924682634	$5.85\times {10}^{-3}$	1.5819	-	7
3	16	0.9881249553	$1.51\times {10}^{-3}$	1.9542	1.4261	7
4	32	0.9869944074	$3.80\times {10}^{-4}$	1.9907	1.9418	7
5	64	0.9867094717	$9.52\times {10}^{-5}$	1.9978	1.9883	7
6	128	0.9866381006	$2.38\times {10}^{-5}$	1.9995	1.9972	7
7	256	0.9866202493	$5.95\times {10}^{-6}$	1.9999	1.9993	7
8	512	0.9866157860	$1.49\times {10}^{-6}$	2.0000	1.9998	7
9	1024	0.9866146701	$3.72\times {10}^{-7}$	2.0000	2.0000	7

and

Table 2. Solutions to (17) where $u_e(5/2)=0.986614298151$ using HOHWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	0.9759511782	$1.07\times {10}^{-2}$	-	-	7
2	8	0.9865266209	$8.77\times {10}^{-5}$	6.9262	-	7
3	16	0.9866085182	$5.78\times {10}^{-6}$	3.9231	7.0127	7
4	32	0.9866139091	$3.89\times {10}^{-7}$	3.8930	3.9252	8
5	64	0.9866142729	$2.53\times {10}^{-8}$	3.9432	3.8894	8
6	128	0.9866142965	$1.61\times {10}^{-9}$	3.9705	3.9413	8
7	256	0.9866142981	$1.02\times {10}^{-10}$	3.9849	3.9695	8
8	512	0.9866142981	$6.40\times {10}^{-12}$	3.9922	3.9844	8
9	1024	0.9866142982	$4.01\times {10}^{-13}$	3.9973	3.9918	8

, respectively. Similar results are obtained for problem 2 (

Table 3. Solutions to (18) where $u_e(5/2)=2.400281157801$ using HWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	2.4071314442	$6.85\times {10}^{-3}$	-	-	6
2	8	2.4126017209	$1.23\times {10}^{-2}$	−0.8468	-	6
3	16	2.4041410973	$3.86\times {10}^{-3}$	1.6744	−0.6292	6
4	32	2.4012885555	$1.01\times {10}^{-3}$	1.9379	1.5685	6
5	64	2.4005355806	$2.54\times {10}^{-4}$	1.9853	1.9216	6
6	128	2.4003449232	$6.38\times {10}^{-5}$	1.9964	1.9816	6
7	256	2.4002971091	$1.60\times {10}^{-5}$	1.9991	1.9955	6
8	512	2.4002851463	$3.99\times {10}^{-6}$	1.9998	1.9989	6
9	1024	2.4002821545	$9.97\times {10}^{-7}$	1.9999	1.9997	6

and

Table 4. Solutions to (18) where $u_e(5/2)=2.400281157801$ using HOHWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	2.3223263560	$7.80\times {10}^{-2}$	-	-	6
2	8	2.3999392337	$3.42\times {10}^{-4}$	7.8328	-	6
3	16	2.4002409325	$4.02\times {10}^{-5}$	3.0875	8.0070	6
4	32	2.4002789788	$2.18\times {10}^{-6}$	4.2063	2.9873	6
5	64	2.4002810290	$1.29\times {10}^{-7}$	4.0807	4.2139	6
6	128	2.4002811500	$7.83\times {10}^{-9}$	4.0392	4.0833	6
7	256	2.4002811573	$4.83\times {10}^{-10}$	4.0193	4.0405	6
8	512	2.4002811578	$3.00\times {10}^{-11}$	4.0094	4.0200	6
9	1024	2.4002811578	$1.87\times {10}^{-12}$	4.0028	4.0098	6

) and problem 3 (

Table 5. Solutions to (19) where $u_e(5/2)=10.9375$ using HWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	10.4949337507	$4.42\times {10}^{-1}$	-	-	9
2	8	10.8281017790	$1.09\times {10}^{-1}$	2.01631	-	10
3	16	10.9102237718	$2.72\times {10}^{-2}$	2.0039	2.0204	11
4	32	10.9306854555	$6.81\times {10}^{-3}$	2.00104	2.0048	10
5	64	10.9357966448	$1.70\times {10}^{-3}$	2.00024	2.0012	10
6	128	10.9370741787	$4.26\times {10}^{-4}$	2.00006	2.0003	10
7	256	10.9373935457	$1.06\times {10}^{-4}$	2.00001	2.0001	10
8	512	10.9374733865	$2.66\times {10}^{-5}$	2.00000	2.0000	11

and

Table 6. Solutions to (19) where $u_e(5/2)=10.9375$ using HOHWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	10.9376310324	$1.31\times {10}^{-4}$	-	-	9
2	8	10.9375106883	$1.07\times {10}^{-5}$	3.6158	-	10
3	16	10.9375006943	$6.94\times {10}^{-7}$	3.9444	3.5900	10
4	32	10.9375000437	$4.37\times {10}^{-8}$	3.9891	3.9413	10
5	64	10.9375000027	$2.74\times {10}^{-9}$	3.9984	3.9884	10
6	128	10.9375000002	$1.71\times {10}^{-10}$	4.0002	3.9983	10
7	256	10.9375000000	$1.07\times {10}^{-11}$	4.0012	4.0002	10
8	512	10.9375000000	$6.63\times {10}^{-13}$	4.0101	4.0006	11

). Finally, the results for problems 4 and 5 along with their corresponding absolute error rates $\Delta u(1/2)=|u_e(1/2)-u(1/2)|$ are given in

Table 7. Solutions to (24) where $u_e(1/2)=0.9242343145$ using HWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	0.92171227335	$2.52\times {10}^{-3}$	0	0	4
2	8	0.92361442525	$6.20\times {10}^{-4}$	2.0245	0	4
3	16	0.92407999821	$1.54\times {10}^{-4}$	2.0061	2.0306	4
4	32	0.92419577632	$3.85\times {10}^{-5}$	2.0015	2.0076	4
5	64	0.92422468252	$9.63\times {10}^{-6}$	2.0004	2.0019	4
6	128	0.92423190668	$2.41\times {10}^{-6}$	2.0001	2.0005	4
7	256	0.92423371257	$6.02\times {10}^{-7}$	2.0000	2.0001	4
8	512	0.92423416403	$1.50\times {10}^{-7}$	2.0000	2.0000	4

–

Table 8. Solutions to (24) where $u_e(1/2)=0.9242343145$ using HOHWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	0.92412675792	$1.08\times {10}^{-4}$	0	0	4
2	8	0.92422791226	$6.40\times {10}^{-6}$	4.0704	0	4
3	16	0.92423391886	$3.96\times {10}^{-7}$	4.0162	4.0739	4
4	32	0.92423428986	$2.47\times {10}^{-8}$	4.0040	4.0171	4
5	64	0.92423431298	$1.54\times {10}^{-9}$	4.0010	4.0042	4
6	128	0.92423431442	$9.63\times {10}^{-11}$	4.0002	4.0010	4
7	256	0.92423431451	$6.02\times {10}^{-12}$	4.0001	4.0003	4
8	512	0.92423431452	$3.76\times {10}^{-13}$	4.0001	4.0001	4

and

Table 9. Solutions to (25) where $u_e(1/2)=3.046376624$ using HWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	3.0425991607	$3.78\times {10}^{-3}$	0	0	5
2	8	3.0455841415	$7.92\times {10}^{-4}$	2.2530	0	5
3	16	3.0461865972	$1.90\times {10}^{-4}$	2.0602	2.3088	5
4	32	3.0463296049	$4.70\times {10}^{-5}$	2.0149	2.0748	5
5	64	3.0463648993	$1.17\times {10}^{-5}$	2.0037	2.0186	5
6	128	3.0463736946	$2.93\times {10}^{-6}$	2.0009	2.0046	5
7	256	3.0463758916	$7.32\times {10}^{-7}$	2.0002	2.0012	5
8	512	3.0463764408	$1.83\times {10}^{-7}$	2.0001	2.0003	5

–

Table 10. Solutions to (25) where $u_e(1/2)=3.046376624$ using HOHWM.

J	$2\mathit{M}$	$\mathit{u}(5/2)$	$|{\mathit{u}}_{\mathit{e}}(5/2)-\mathit{u}(5/2)|$	${\mathit{k}}_{\mathit{J}}^{\mathit{e}}$	${\mathit{k}}_{\mathit{J}}$	${\mathit{r}}_{\mathit{f}}$
1	4	3.0485346088	$2.16\times {10}^{-3}$	0	0	5
2	8	3.0464549533	$7.83\times {10}^{-5}$	4.7840	0	5
3	16	3.0463809046	$4.28\times {10}^{-6}$	4.1936	4.8117	5
4	32	3.0463768823	$2.58\times {10}^{-7}$	4.0499	4.2024	5
5	64	3.0463766398	$1.60\times {10}^{-8}$	4.0126	4.0524	5
6	128	3.0463766248	$9.99\times {10}^{-10}$	4.0032	4.0132	5
7	256	3.0463766239	$6.24\times {10}^{-11}$	4.0008	4.0033	5
8	512	3.0463766238	$3.90\times {10}^{-12}$	4.0006	4.0008	5

, respectively.

All model problems showed that the HOHWM outperforms HWM considerably. The rate of convergence of the HOHWM tends to four, i.e., doubles that of the HWM. In the case of sample problem 1, the accuracy of $10\times {10}^{-7}$ was not achieved by the HWM even with 512 collocation points, but by employing the HOHWM, such accuracy was reached with 32 collocation points. Similar results are also obtained in the case of the other sample problems. The number of iterations necessary to arrive at the solution corresponding to the HWM and the HOHWM is same in the case of the second and the third sample problems. In the case of the first sample problem, the number of iterations of the HOHWM is higher for some resolutions, but just by one iteration.

When considering the computation times needed to reach the same accuracy using the HWM and HOHWM (

Table 11. CPU times for solutions of Equations (17) and (18) at the same level of accuracy.

Equation (17)		Equation (18)
HWM, $J=9$	HOHWM, $J=4$	HWM, $J=8$	HOHWM, $J=4$
6.59283	1.8728	3.42749	1.87265

), one can see that the same level of accuracy is achieved using less CPU time by employing the HOWHM instead of the HWM. This is because the size of the matrices involved is considerably smaller in case of the HOHWM.

6. Conclusions

The HOHWM has been implemented for solving nonlinear differential equations. Three particular Riccati equations and two Liénard equations with known exact solutions are examined.

The iterated numerical solution in case of the HOHWM can depend strictly on the value of the previous iterated solution at the boundary (Equations (33) and (34)). That is not the case when using the HWM. These dependencies are caused by satisfying differential equations in points, determined by the HOHWM algorithms (boundary points for $s=1$, etc.).

Both, the HWM and the HOHWM produce numerical solutions which were shown to be in good agreement with the exact solution. However, the HOHWM has shown principally higher accuracy in the case of the same mesh.

Computationally, the most expensive operation performed is the calculation of the solution of the system of algebraic equations. However, the algebraic systems of equations corresponding to the HWM and the HOHWM, have the same dimensions and symmetry properties. From a practical point of view, it is interesting to estimate computational complexity of the solutions providing the same accuracy. Thus, for example in the case of model problem 1, to obtain the same accuracy, one can use the HWM with 512 collocation points and the HOHWM with just 16–32 collocation points. Thus, by applying the HOHWM, it is possible to obtain the results with the same accuracy as by applying the HWM with a substantially lower computational cost. It is correct to note that the conclusions made are based on solution of five sample problems related to the Riccati and Liénard equations. Therefore, the study of nonlinear differential equations by employing the HOHWM needs further attention.

Future study should be related to development and application of the HOHWM for analysis of the nonlinear solid mechanics problems [69] and simulation of solitons and shock waves, which are known for their complex behaviour and providing challenges for any new numerical method [70,71].