Higher Order Haar Wavelet Method for Solving Differential Equations

1. Introduction

Wavelets are most commonly used in signal processing applications to denoise the real signal, to cut a signal into different frequency components, to analyze the components with a resolution matched to its scale, also in image compression, earthquake prediction and other algorithms.

However, the current study is focused on the area where the use of wavelet methods shows a growth trend, i.e., in the solution of differential equations. Many different wavelets based methods have been introduced for solving differential and integro-differential equations. The Legendre wavelets are utilized to solve fractional differential equations in [1, 2, 3, 4] and integro-differential equations in [5, 6]. In [7, 8], the Daubechies wavelet based approximation algorithms are derived to solve ordinary and partial differential equations. In [9], the Lucas wavelets are combined with Legendre–Gauss quadrature for solving fractional Fredholm–Volterra integro-differential equations. The series solution of partial differential equations through separation of variables is developed by using the Fourier wavelets in [10]. The Riesz wavelets- based method for solving singular fractional integro-differential equations was developed in [11]. In the studies in [12], the Galerkin method was combined with the quadratic spline wavelets for solving Fredholm linear integral equations and second-order integro-differential equations. The Chebyshev wavelets method for partial differential equations with boundary conditions of the telegraph type is examined in [13].

The simplest of all wavelet-based approaches was introduced by Alfred Haar already in 1910 [14]. The Haar wavelet-based approach for solving differential and integro-differential equations was introduced in 1997 [15, 16]. Based on the Haar wavelet method (HWM), Chen and Hsiao in [15, 16] proposed an approach where the higher order derivative involved in the differential or integro-differential equations is expanded into the series of Haar wavelets. This approach is based on the nature of Haar functions. Due to the piece-wise constant nature of the Haar functions they are not differentiable but are integrable. In [15, 16], the problems of the lumped and distributed parameter system and those of linear time delayed systems were solved. The Chen and Hsiao approach-based HWM was adapted successfully for solving a wide class of differential, integro-differential and integral equations [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45]. Pioneering work in the development of Haar wavelet-based techniques was conducted by Lepik [17, 18, 19, 20, 21, 22, 23], covering ordinary and partial differential equations [17, 19, 21], integro-differential equations [18], integral equations [20], and fractional integral equations [22]. The HWM approaches and their applications are summarized in a monograph [23]. The HWM is adapted for the analysis of nonlinear integral and integro-differential equations in [24, 25, 26, 27], covering one- and multi-dimensional problems. Solid mechanics, particularly composite structures, are examined using the HWM in [28, 29, 30, 31, 32, 33]. These studies cover free vibration analysis of orthotropic plates [28], functionally graded composite structures [30, 31, 32], delamination detection in composite beams [29], and other structures.

Some recent trends in the development and application of the HWM can be outlined as solutions of fractional differential and integro-diffrential equations [34, 35, 36, 37, 38] as well as the development of a non-uniform and adaptive grid. In the case of fractional differential or integro-differential equations, two principally different HWM approaches regarding to wavelet expansion are available in the literature. The aim of the first approach is to expand the highest order fractional derivative included in the differential equation directly into Haar wavelets, i.e., direct conversion of the Chen and Hsiao approach for fractional differential equations. In [34, 35, 36, 37, 38], the Haar wavelet operational matrix of fractional order integration is introduced and implemented for solving differential and integro- differential equations. The aim of the second approach is to utilize the definitions of fractional derivatives (Caputo derivative, etc.) and convert fractional differential terms into integrals, which contain integer derivatives only. Such an approach has been introduced by Lepik in [22] and utilized in a number of papers [39, 40, 41]. The two approaches considered are implemented and compared in [42]. It is pointed out in [42] that the two approaches have the same rate of convergence if the order of the fractional derivative exceeds one (α > 1). However, if the order of the fractional derivative is less than one (α < 1), the second approach has the rate of convergence equal to two, but the rate of convergence of the first approach is 1 + α, i.e., less than two. Thus, in the case of α < 1, the second approach has a higher convergence rate and can be preferred.

HWM with a nonuniform grid was introduced in [43] using a proportionally changing grid size. The same approach was utilized for the free vibration analysis of non-uniform axially graded beams in [44] and for solving singularly perturbed differential difference equations of neuronal variability in [45].

In most of the studies [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45], it was concluded that HWM is simple to implement. In the review paper [46], it was pointed out that the HWM is efficient and powerful in solving a wide class of linear and nonlinear reaction–diffusion equations. However, the convergence theorem and accuracy estimates derived for the HWM in [47, 48] state that the order of convergence of the Chen and Hsiao approach-based HWM is equal to two. The latter result is rather modest in the context of engineering. Comparison of the HWM with widely used numerical methods in engineering reveals that HWM needs principal improvement in order to compete with the differential quadrature method (DQM) [49].

The HOHWM as an improvement of HWM was recently introduced by Majak et al. in [50]. The convergence rate of the method was improved from 2 to 2 + 2 s, where s stands for the method parameter. This new method is currently underused, but the first results obtained have shown that a principal growth of the accuracy can be achieved with a minimum growth of complexity [51, 52, 53]. In [52], the free vibrations analysis of the Euler-Bernoulli nanobeam was performed. Figure 1 shows the numerical complexity estimates of the HWM and HOHWM solutions yielding a similar absolute error. Here the numerical complexity is determined by number of main operations of the most complex subtask - solution of discrete algebraic system of equations [52].

The logarithmic scale is used in Figure 1 since the complexity of the HWM appears several orders higher (10ˆ8) than that of the HOHWM (10ˆ3…10ˆ5). These results were obtained using the method parameter s = 1 (i.e., fourth order convergence). In practice, one of most important factors is the computational cost. In the case of the considered problem, the computational cost of the HOHWM solution is 10ˆ3…10ˆ5 times lower than that of the HWM. The obtained results hold good in the case of all four boundary conditions considered: pinned-pinned (P–P), clamped-pinned (C-P), clamped-clamped (C-C), and clamped-free (C-F).

2. Theoretical basis of the HWM and the HOHWM

This section introduces the Haar functions and presents the theoretical basis of both, the HWM and the HOHWM, covering basic principles, algorithms, convergence and accuracy issues.

3. Numerical convergence analysis and Richardson extrapolation

The derivations of the numerical estimates of the order of convergence, as well as extrapolation formulas can be found in [54] and are omitted herein for the sake of conciseness. Let us denote the numerical solutions on a sequence of nested grids by Fi−2,Fi−1, Fi, corresponding to grid sizes hi−2/hi−1=hi−1/hi=2. Then the order of convergence of the numerical method can be estimated by the formula.

if the exact solution Fexactis known. If the exact solution is unknown, the following formula can be employed [54].

The accuracy of the results can be improved by employing the Richardson extrapolation formula as [54].

The accuracy of the extrapolated results Rican be estimated by applying formulas (28) or (29) to these results. The numerical rates of convergence computed for case study, described in Section 5.1, are depicted in Figure 2.

The numerical rates of convergence determined are in agreement with convergence theorem for the HWM (Section 2.2.2) and relation given for the HOHWM in Section 2.3 (the rate of convergence of the HOHWM is equal to 2 + 2 s).

4. Complexity analysis

As pointed out above, the accuracy and numerical/time complexity are two key characteristics for any numerical method, algorithm. The computing time is often used as a measure of complexity of algorithms using particular software. More general approach applied commonly in algorithm theory is to estimate the number of basic operations required by each algorithm. The latter approach is independent of software used and does not even require execution of the algorithms. For this reason, in the current study the numerical complexity of the algorithms is estimated based on the number of basic operations.

According to the HWM and the HOHWM algorithms the solution of the differential equation is obtained from the solution of the discrete algebraic system of equations and certain additional operations for composing the linear system and evaluation of the solution in given points. The mentioned additional operations are similar for both methods and have lower asymptotic complexity than the solution of the algebraic system of equations. Thus, the numerical complexity of the HWM and the HOHWM can be compared based on number of basic operations needed for solving algebraic system of equations determined by the rank of the algebraic system of equations (systems are similar by structure).

In the case of the same number of collocations points N, the ranks of the algebraic systems corresponding to the HWM and the HOHWM are equal to N and N + 2 s (here s = 1, 2 or 3), respectively. Furthermore, in the cases where the 2 s complementary integrations constants are determined analytically, the rank of the algebraic system of equations of the HOHWM reduces to N. Thus, in the case of the same mesh used the numerical complexity of the HWM and HOHWM is similar (or equal depending on implementation). However, these solutions have principally different accuracy (see Tables 1-5) and such comparison is rather theoretical.

In practice, it is important to compare methods, providing the same accuracy. In the following the given accuracy is fixed by absolute error less than 2.0e-10 and the complexities of the HWM and the HOHWM are compared in Figure 3. The logarithmic scale is used in Figure 3, since the complexities of the HWM and the HOHWM differ by several magnitudes.

In the case of the HWM the absolute error 2.0e-10 was reached by use of 16,384 collocation points (corresponding algebraic system has 16,384 equations). In the case of the HOHWM the same accuracy was achieved by use just 64, 16 or 4 collocations points corresponding to the s = 1, s = 2 or s = 3, respectively (i.e. the algebraic system needed to solve is reduced to 64, 16 or 4 equations). Thus, it can be concluded, that making use of the HOHWM instead of the HWM will lead to principal reduction of numerical complexity of the solution.

The practical value of the developed HOHWM approach is reduction of computational cost of the solution by several magnitudes (directly determined by numerical complexity). It should be noted that making use of the HOHWM instead of the HWM, especially in the cases s > 1 will increase implementation complexity, but not substantially.

5. Case studies

In the following, the two case studies are performed in order to validate the accuracy and convergence of the recently introduced HOHWM and compare results with HWM.

5.1 Linear ordinary differential equations

As a rule, the new methods are validated on the samples where the exact solution is known. Herein, the linear ordinary differential equations are considered as the first sample problem. Let us consider a sample problem solved in [17] by applying the HWM

Gtutdutdtd2utdt2=d2utdt2+αdutdt+βut−γft=0,E31

where ftis a given function (ft=cos2t), p, qand rare constant parameters (α=0.05, β=0.15, γ=1). The initial conditions u0=0,dudt0=1are utilized.

In the case of the HWM, the second order derivative is expanded into Haar wavelets as

d2utdt2=aH.E32

In Eq. (32), aand Hstand for the coefficient vector (row vector) and the discrete Haar matrix given by formulas (6), respectively. The solution of the differential Eq. (31) is obtained by integrating relation (32) twice with respect to tand satisfying initial conditions

ut=aP2+t,E33

where the elements of the matrix P2are defined by (7) and the coefficient vector ais determined by substituting the solution (33) and its derivatives in Eq. (31) as

a=γft−α−βtH+αP1+βP2−1.E34

In the case of the HOHWM and s = 1, the fourth order derivative is expanded into Haar wavelets as

d4utdt4=aH.E35

The solution of the differential Eq. (31) is obtained by integrating relation (35) four times with respect to tand satisfying initial conditions

ut=aP4+c3t36+c2t22+t,E36

The remaining two integration constants in Eq. (36) can be determined by satisfying Eq. (31) at the boundary points t=0and t=1.The latter two algebraic equations can be added to the algebraic system obtained by substituting the solution (36) and its derivatives in Eq. (31). In the latter case, the algebraic system includes 2M+2equations. An alternate approach is to determine the remaining two integration constants analytically from the same conditions and replace to algebraic system.

The numerical results obtained by utilizing the HWM and the HOHWM are compared in Tables 1-3.

It can be observed from Tables 1-3 that in the case of the HWM, the order of convergence tends to two and in the case of the HOHWM, it tends to 2 + 2 s, i.e., to four if s = 1, to six if s = 2 and to eight if s = 3. Use of HOHWM provides a principal increase of accuracy. The maximum accuracy obtained by the use of the HWM at 256 collocation points (N = 256) has been achieved by using the HOHWM at 16 collocation points if s = 1, and at 4 collocation points if s = 2. In the case of the HOHWM and s = 3, the accuracy achieved at 4 collocation points was significantly higher than that of the HWM with 256 collocation points.

In Figure 4 are shown the error ratios for different mesh (N = 4,16,64 and 256). The absolute error of the HWM is divided by error of the HOHWM, where blue, green and gray colors correspond to the HOHWM parameter s values 1,2 and 3, respectively. Thus, in the case of mesh N = 4, making use of the HOHWM instead of the HWM reduced the absolute error 8.91E+01 (s = 1) to 9.83E+06 (s = 3) times. In the case of mesh N = 256, the use of the HOHWM reduced the absolute error 2.59E+05 (s = 1) to 4.53E+15 (s = 3) times. Since the error ratio depends strongly on the mesh used, the logarithmic scale was used in Figure 4.

The numerical analysis is performed using MATLAB software. Since the accuracy achieved by the use of the HOHWM in the case of s = 2 and s = 3 exceeds the limits of the double precision computing, the variable precision computing (VPA) was used.

Note that this is needed only in the case of particular problems and large mesh where the accuracy exceeds the limits of double precision computing.

6. Conclusions

The HOHWM introduced recently by authors as an improvement of the HWM in order to compete with the numerical methods widely used in engineering. It was shown that using the HOHWM instead of the HWM will improve principally the accuracy of the solution and increase the rate of convergence in the case of all problems studied. It was found that the rate of convergence of the HOHWM depends on the model parameter s and is equal to 2 + 2 s.

From a practical point of view, it is important that the HOHWM can achieve the same accuracy as the HWM with significantly lower mesh and reduced computational cost.

In the simplest case of the HOHWM where s = 1, the order of the convergence of the HOHWM is equal to four. The user can select suitable s value depending on the accuracy requirements of a particular problem considered.

In future study, the new method proposed can be extended/adapted for solving a wide class of differential and integro-differential equations, including fractional differential equations, multidimensional problems, nonlinear boundary value problems arising in engineering design.

Acknowledgments

The study was supported by Estonian Centre of Excellence in Zero Energy and Resource Efficient Smart Buildings and Districts, ZEBE, TK146 funded by the European Regional Development Fund (grant 2014–2020.4.01.15–0016).