Application of Higher Order Haar Wavelet Method for Solving Nonlinear Evolution Equations

The recently introduced higher order Haar wavelet method is treated for solving evolution equations. The wave equation, the Burgers’ equations and the Korteweg-de Vries equation are considered as model problems. The detailed analysis of the accuracy of the Haar wavelet method and the higher order Haar wavelet method is performed. The obtained results are validated against the exact solutions.


Introduction
From the mathematical viewpoint evolution equations under consideration are partial differential equations (PDEs). In case of numerical solution of PDEs one has to approximate the partial derivatives with respect to time and space coordinates. Local or global methods can be used for these approximations. The finite difference method is the most typical example of local methods. In the case of global methods, usually the function itself is approximated as a sum of basis functions (for example trigonometric functions). However, trigonometric functions are not the only choice for basis functions. For example, in [64] mixed Laguerre-Legendre interpolation approximation is applied. One of the more popular global methods is the fast Fourier transform (FFT) related pseudospectral method [20,22,55,56,58,59]. In the present paper another global approach is applied: Haar wavelets are used for approximation of partial derivatives.
In the current study the HOHWM is treated for solving evolution equations. In addition, the wave equation was chosen as the first model equation for its simplicity. The Burgers' equation [10,11] has previosuly been studied using the HWM [29,36,41,49] and also has an anlytical solution, thus the possibility of comparison made it a decent choice for a model equation. It finds applications in modeling of turbulence [10] and traffic [46], as well as in weak non-stationary shock wave in real fluids [31] and in nonlinear acoustics [24]. The HWM has also previously been applied to the Korteweg-de Vries (KdV) equation in [50,52]. The KdV equation was derived in order to describe the movement of long unidirectional shallow water waves in a rectangular channel but has been found to model different nonlinear phenomena nowadays [18,19,54]. A detailed analysis of the improved accuracy of HOHWM over HWM is performed for all three model equations.
The article is structured as follows. Firstly, the Haar wavelet family is introduced in Section 2. Secondly, the model equations used are described in Section 3. Section 4 outlines the HWM and HOHWM. The analysis of the results is presented in Section 5 and conclusions are drawn in the final section.

Haar wavelet family
In the following the Haar wavelet family is defined utilizing the notation introduced by Lepik in [38]. 2M subintervals of equal length ∆x = (B − A)/(2M ) form the integration domain x ∈ [A, B]. The maximum level of resolution J is defined as J = log 2 (M ). The Haar wavelet family for a fixed M can be described as

M. Ratas
Thus, any square integrable function f (x) can be expanded into Haar wavelets as where a i denote the Haar coefficients. According to [38] the integrals of the Haar functions (2.1) of order n can be calcualted analytically as follows , for x ∈ ξ 3 (i), B . (2.4) Within the present paper the matrix form of the above formulation is used. Therefore, the elements of (2M ) × (2M ) matrix H, are given as values of the Haar functions denotes the nth integral of the Haar wavelet matrix for a given resolution J. Using (2.5) and (2.6) and considering the coefficient vector a one obtains f (x) = a · H instead of (2.3) and · · · It must be noted, that the matrices H and P n depend on vector x of collocation points. The statements in (2.4) imply that in boundary points A and B hold P n (A) i = 0, ∀n > 0, ∀i, (2.7) P n (0) and P n (B) form column vectors which will be used when forming a system of equations that satsify the boundary conditions. These expressions can often be simplified due to equation (2.7).

Model equations
The wave equation subject to the initial and boundary conditions was considered as the first model equation. The particular exact solution used was a travelling wave solution of the form where c 1 is a parameter that varies the steepness of the shockwave, c is its travelling speed and x 0 its initial phase. As the second model equation, the Burgers' equation was chosen. Here ν is the diffusion coefficient. The Burgers' equation subjected to the initial and boundary conditions was considered. Its analytical solution has been found in the form of where R 0 = u 0 L/ν is the Reynolds number, E n = νn 2 π 2 /L 2 , I n represents the modified Bessel functions of first kind and L = B − A is the x domain range [9,16].
The third model equation was chosen to be the KdV equation of the form subjected to the initial and boundary conditions An exact one soliton solution for the KdV equation (3.5) is known in the form of where c is the speed of the travelling soliton, α is the nonlinear coefficient in the equation and x 0 denotes the initial phase. There exists also a two-soliton solution [1,32,57,70] in the form where , and c B and c S are the speeds of the bigger and smaller soliton, respectively, x 0B and x 0S are the initial phases for the bigger and smaller soliton, respectively and α is the nonlinear parameter. Using (3.7) soliton interactions can be observed.

Haar wavelet methods
In many numerical studies [2,12,25,41], spatial derivatives are expanded into Haar wavelet series while finite difference type schemes are used for integration with respect to time. In the present study, MATLAB's ode45 [61] solver based on Runge-Kutta (4,5) formula [17] is used for integration with respect to time. This allows one to calculate the value of the function at each time moment from its values at previous time moments without the need to expand both axes into the Haar series. The latter often results in extremely large matrices which can get computationally expensive. The well known HWM involves expanding in the Haar wavelet series the highest order derivative present within the equation. One then obtains for the nth spatial derivative where the subscript notes the axis along which the HWM (as well as the HO-HWM) is deployed and a is the Haar wavelet coefficient vector. After integrating n times one arrives at the function itself as where the unkown coefficients c i can be calculated by using the boundary conditions and x i denotes the collocation vector with its elements raised to the power of i. The superscript of the vector of collocation points will refer to the element wise multiplication throughout the rest of the present study.
When it comes to the HOHWM, one needs to start with a derivative of a higher order than that which is present in the equation. In general, an additional 2s derivatives are used. However, in the current study, only s = 1 is considered. Thus, one starts with and after integration arrives at where two extra coefficients, c n and c n+1 are introduced. In order to calculate those extra coefficients, some extra information is needed. In the present study the equation is evaluated at the boundary. When the unkown coefficients c i have been obtained, the function u can be described as where R xm is a matrix and S xm a vector obtained only from the Haar matrices and boundary conditions. They both depend on x. S xm can also depend on t (depending on the boundary conditions). Thus for the kth derivative (k < m) of u one can write It must be noted that R x0 = H x and S x0 = 0 due to (4.1) and (4.2) in case of HWM. Once R xm and S xm have been found, one can arrive at where m = n (HWM) or m = n + 2 (HOHWM). Equation (4.4) can be used to find the wavelet coefficient vector a from the previous iteration. Finally, the equation in question will need to be arranged in the form u t = f (u, u x , u xx , ..., u nx ) or u tt = f (u, u x , u xx , ..., u nx ). Then, (4.4) can be used to calculate a, after which (4.3) can be used to substitute u and its derivatvies into the right hand side of the equation in question. This routine has to be followed for each time step.
In the following subsections, the numerical approach is described for each model equation. Firstly, the HWM is described. Then the extra conditions to find the extra coefficients added by HOHWM are introduced. Finally, HOHWM is described. The domain is fixed as A = 0 and B = 1.

M. Ratas
one again obtains the same results (4.6) for HWM. However, for HOHWM the equation is again evaluated at the boundary points as Since the boundary conditions (4.8) are homogeneous u t (0, t) = u t (1, t) = 0 and (4.9) simplifies to The addition of (4.10) gives

KdV equation
Given homogeneous boundary conditions for HWM one obtains For HOHWM the equation is yet again evaluated at the boundary points as Since the boundary conditions (4.11) are homogeneous u t (0, t) = u t (1, t) = 0 and the above simplifies to The addition of (4.12) gives

Numerical results
Numerical experiments were carried out for equations  Table 1.  Table 2.
Results for the Burgers' equation: maximal time t f at which the numerical solutions deviates from the exact solution less than 10 −3 against the resolution parameter J with initial conditions (5.2) and boundary conditions (4.8) (x ∈ [0, 1], t ∈ [0, 1], ν = and boundary conditions (4.8) were used. In case of HOHWM the additional conditions (4.10) were added. Due to the near-singular nature of the solution, the numerical solution can get unstable near the steep slope. Thus the time at which sufficient accuracy was lost t n f is shown instead of the maximum error. The numerical result was taken to be sufficiently accurate if the maximum deviation with respect to the exact solution was smaller than 10 −3 . Corresponding results can be seen in Table 2. Figure 1 shows how the maximum error behaves in time for different resolutions at different values of ν.
and boundary conditions (4.11) were used. In case of HOHWM the additional conditions (4.12) were added. The maximum deviation from the exact solution can be seen in Table 3.
In case of the KdV equation, the two-soliton exact solution (3.7) was also used. In this case the initial condition was taken from the exact solution as Boundary conditions (4.11) as well as additional conditions (4.12) were used in this case as well. The results of these calculations can be found in Tables 5  and 6.  Table 4. Results for the KdV equation: maximal deviation of the numerical solution from the exact solution max∆u and maximal relative deviation max rel ∆u against the resolution parameter J with initial conditions (5.3) and boundary conditions (4.11) (α = 6, c = 1000, Typical numerical results for the wave equation, the Burgers' equation and the KdV equation can be seen in Figure 2.
In case of both the Burgers' equation as well as the KdV equation the steep slope can be observed. The soliton interaction and subsequent phase shift can also be observed in case of the KdV interaction. Table 5.
Results for the KdV equation soliton interaction: maximal deviation of the numerical solution from the exact solution max∆u and maximal relative deviation max rel ∆u against the resolution parameter J with initial conditions (5.4) and boundary conditions (4.11) (α = 1, c B = 10000, x 0B = 1/5, c S = 10000/3, Tables 1-6 show that HOHWM generally gives a more accurate result than HWM. In case of the wave equation (Table 1), the increase in c 1 (which translates into an increase in the slope) also increases the error of both methods while HOHWM remains more accurate. The Burgers' equation at the values of ν that were used introduces such a steep slope that at lower resolutions,  Table 6.
Results for the KdV equation soliton interaction: maximal deviation of the numerical solution from the exact solution max∆u and maximal relative deviation max rel ∆u against the resolution parameter J with initial conditions (5.4) and boundary conditions (4.11) (α = 6, c B = 10000, x 0B = 1/5, c S = 10000/3, neither method is able to successfully calculate the solution. This is due to the fact that the collocation points are used and at a low resolution, no collocation point falls within the slope. However, it is clear that HOHWM is able to perform calculations further in time before diminishing in accuracy, making it advantageous here as well. Figure 1 shows that HOHWM remains stable for a longer time than HWM since the former does not reach such high error throughout the integration. In case of the one soliton solutions of the KdV equation, HOHWM again shows a higher accuracy than HWM. It can be noted from Tables 3 and 4 that the HOHWM is more accurate at J = 5 (2M = 64) than HWM is at J = 7 (2M = 256).
While the results are not directly comperable with [50] due to the difference in domain as well as values of c and differences in presenting error analysis, some comparison can be made. The maximum relative deviation from the exact solution can be compared. Such a comparison between [50] and results in Table  4 is carried out in Table 7. While the result with HWM is not as accurate as that of [50] the results with HOHWM surpass the referenced results. The KdV soliton interactions show an interesting phenomenon, however. From Tables 5 and 6 one can see that at lower resolutions neither method can successfully solve the problem. This is caused by the fact that one now has two relatively localized solitary waves and at lower resolutions no collocation point falls within these localized solitary waves. However, at higher resolutions, one can clearly see the advantages of HOHWM over HWM as the former shows a significant boost in accuracy.

Conclusions
The HOHWM has been adapted for solving partial differential equations numerically. The results obtained by the widely used HWM were compared with those obtained by the HOHWM. The Burgers' equation and the KdV equation were considered as model equations.
Our numerical experiments demonstrated that both methods, the HWM and the HOHWM, were able to provide numerical solutions that are in good agreement with the exact analytical solutions. The detailed comparison of the accuracy of the HWM and the HOHWM was performed in Section 5. It is shown that HOHWM can be preferred where high accuracy is required. Furthermore, the HOHWM can also be preferred in cases where high accuracy is not required by applying the method at a lower resolution in comparison with HWM.
Solving the system of algebraic equations is the most computationally expensive task within Haar wavelet based methods. However, the algebraic systems of equations are of the same dimension and have the same symmetric properties for the HWM as well as for the HOHWM. Even though some expressions needed to evaluate are more complex in case of HOHWM, in general it can be concluded that the numerical complexity of the HOHWM is only slightly higher than that of the HWM at the same resolution. It would be more pragmatic to estimate computational complexity of the solution providing the same accuracy, however. For example, in the case of single KdV-soliton solution (Tables 3 and 4), one can use HWM with 256 collocation points or HOHWM with 64 collocation points in order to obtain the same degeree of accuracy. Solving a 64 × 64 algebraic system is substantially computationally cheaper than solving a 256 × 256 system. Therefore, it is possible to obtain the results with the same accuracy as HWM with lower computational cost by applying HOHWM.
The Fourier transform related pseodospectral method is known as a powerful tool for numerical solution of evolution equations because it is able to produce high accuracy at relatively low number of collocation points [22]. However, M. Ratas this method has a disadvantage: because of the nature of the Fourier transform one must apply periodic boundary conditions. Haar wavelet related methods do not have such a disadvantage and one can solve PDEs numerically applying arbitrary boundary conditions in case of the HWM and arbitrary time independent boundary conditions in case of the HOHWM. In order to use time dependent boundary conditions with HOHWM a different numerical scheme must be applied to approximate the temporal derivatives.