Solving Nonlinear PDEs Using the Higher Order Haar Wavelet Method on Nonuniform and Adaptive Grids

The higher order Haar wavelet method (HOHWM) is used with a nonuniform grid to solve nonlinear partial differential equations numerically. The Burgers’ equation, the Korteweg–de Vries equation, the modified Korteweg–de Vries equation and the sine–Gordon equation are used as model equations. Adaptive as well as nonadaptive nonuniform grids are developed and used to solve the model equations numerically. The numerical results are compared to the known analytical solutions as well as to the numerical solutions obtained by application of the HOHWM on a uniform grid. The proposed methods of using nonuniform grid are shown to significantly increase the accuracy of the HOHWM at the same number of grid points.

1 Introduction of piecewise constant functions and are therefore mathematically the simplest of all the wavelet families. They can be integrated analytically arbitrary times.
The Haar Wavelet Method (HWM) was originally proposed in [26] for solving differential equations. It was later extended for solving a wide class of integro-differential and integral equations [4,5,6,23,53,57,83]. According to this method, the highest order derivative included in the differential equations is expanded into the Haar series. Lepik developed the integration techniques for the HWM in [52,53,54,57,58]. A thorough overview of the HWM and its applications can be found in [59]. The weak formulation based HWM was introduced and complexity issues of the strong as well as weak formulation based HWM were discussed in [61].
Recently, the higher order Haar wavelet method (HOHWM) has been developed [63] as an improvement of the HWM. The HOHWM has been found to increase accuracy as well as convergence over the regular HWM [42,43,46,62,63]. In case of the HWM the highest derivative present within the differential equation is expanded into the Haar series. However, the HOHWM proposes that a derivative of 2s higher order is expanded into the Haar series. This procedure introduces 2s extra integration coefficients, although, in case of s = 1, the equation can be evaluated at its boundary in order to solve for those extra coefficients. Therefore, this paper focuses on the use of the HOHWM over the conventional HWM.
The nonuniform Haar wavelets were first introduced in [30]. They were first used to solve integral and differential equations by Lepik [55]. The nonuniform Haar wavelets have since been used in multiresolution analysis [3], boundary value problems [32], fractional order problems [64,78] as well as two dimensional problems [68]. An overview of uniform and non-uniform wavelet based methods can be found in [48].
In order to deal with PDEs with rapid solution variation, adaptive grids have been developed [10,15,31,75]. They have seen use in many areas, including astrophysics and turbulence problems in hydrodynamics [35,47,81,84], among others. The general idea of adaptive grids is to be able to reshape the grid in such a way that areas with large variations have more grid points. Such a system uses error estimates as weight functions to determine where the grid needs to be more concentrated. Adaptivity can be achieved by changing the number of grid points or the movement of grid points. Since the HWM and the HOHWM require a fixed number of grid points, movement of grid points is considered within this paper. With a fixed number of grid points resolution is raised locally at the expense of decreased resolution in other regions. The effect of the decreased resolution can be minor if the other regions previously have more grid points than are needed for the required accuracy. The basic premise of an adaptive grid is that when moving from coordinate grid ξ to s with a weight function of w, wds = cdξ holds for a proportionality constant c.
In the current study, the HOHWM is used to solve nonlinear partial differential equations which have abrupt changes in their solutions. The benefits of the nonuniform grid can best be observed in solutions with such abrupt changes. To the best knowledge of the authors, this is the first time the HOHWM is used on a nonuniform grid. The Burgers' equation [17,18] and the Kortewegde Vries (KdV) equation were chosen as model equations because they have previously been studied using the HOHWM on the uniform grid [76] and have analytical solutions [2,7,49,67] which make comparisons easy to draw. The Burgers' equation finds applications in modelling of turbulence [17] and traffic flow [66] as well as in nonlinear acoustics [37] and non-stationary shock waves in fluids [50]. The KdV equation was first used as a description of the propagation of long unidirectional shallow water waves in a rectangular channel, but it has been used to model different nonlinear phenomena since then [29]. The modified Koretweg-de Vries equation (mKdV) was chosen as a model equation because of its close relation to the KdV equation as well as the existence of an analytical solution [1,41]. The mKdV equation finds use in modelling of behaviour of anharmonic lattices [41]. The sine-Gordon equation was chosen as the last model equation because of the existance of an analytical traveling wave solution [51] and the fact that its analytical solution is not asymptotically homogeneous like those of the KdV and mKdV. The sine-Gordon equation was originally developed to describe surfaces with constant mean curvature [13] and has also seen use in one-dimensional crystal dislocation theory [33,34] as well as in many other fields [14,65,79,80]. In the current study, the adaptive grid with a constant number of grid points according to the constrained least-squares statement defined in [31] is used. This paper is structured as follows. Firstly, the spatial discretization and nonuniform grids are described in Section 2. Section 3 focuses on the Haar wavelet family and describes how it is used in the HOHWM. Model equations and their corresponding exact solutions are introduced in Section 4, and the nonuniform HOHWM is applied to them in Section 5. The numerical results are discussed in Section 6 and conclusions are drawn in the final section.

Discretization
The Haar wavelet method is used along the spatial axis only within this paper. The ODEPACK [40] provided by SciPy [85] is used for integration with respect to time. It automatically switches between Adams and BDF [8,28,36] solvers according to [74]. Since this involves calculating the value of the function at each moment in time, this allows for adaptive changes to the grid.
The collocation method is used alongside the HOHWM. As such, the collocation points lie in the middle of an adjacent pair of grid points. Given the domain x ∈ [A, B], the 2M + 1 grid points can be described as x g (l), l = 0, 1, . . . , 2M, x g (l + 1) > x g (l) ∀l, x g (0) = A, x g (2M ) = B.
(2.1) Within this paper, three different approaches are used for numerical differentiation with respect to the space coordinate. Two of those approaches are static (the space grid does not change during usage) and the third is adaptive (the space grid changes during operation).

Nonuniform grids
The first type of nonuniform grid for domain x ∈ [0, 1] follows the formula where q is an arbitrary constant. It is clear that for q < 1 the above leads to a grid which concentrates grid points around 1 and that in case of q > 1 the grid points are concentrated around 0. It is also clear that given q → 1 the grid (2.2) becomes uniform. This type of grid is most applicable for cases where the characteristic behaviour of a solution tends to concentrate near a boundary. The same type of nonuniform grid was used in [68]. The second type of nonuniform grid for domain x ∈ [0, 1] is perhaps a little simpler. It follows the formula where γ is the gap parameter. It describes how far from either boundary the coarse part of the grid ends and the finer part begins. Such a grid is most applicable for cases where the characteristic behaviour of a solution stays in the middle of the domain, yet its exact position cannot be determined.

Adaptive grid
In the current paper the adaptive grid based on the constrained least-squares statement [31] is used. According to it, when changing the grid, the new grid x will have grid points such that where x min and x max denote the minimal and maximal grid point, respectively and w is the weight function. It must be noted, that w i+1/2 is used. In the current paper, the weight function is calculated at collocation points, which lie in the middle of two subsequent grid points according to (2.1) and thus the values at the collocation points can directly be used. For the weight function, either the function itself or its first derivative is scaled and a constant e is added (to ensure the weights never reach 0). The weight function from the function u is thus obtained as where d is the scaling factor and e is the offset. Similarly, when using the derivative one obtains The algorithm (2.4) is carried out iteratively as many times as necessary for conversion of the grid points. During each iteration, the weights are estimated at the new collocation points by interpolation.
In case of the adaptive grid, a new grid is only calculated if the characteristic point in the numerical result (i.e the maximal point of the function or its derivative) has moved by δ. This helps to save on CPU time.

Haar wavelet family
In the following, the Haar wavelet family is defined utilizing the notation introduced by Lepik in [55] and Oruç in [68]. The integration domain [A, B] can be divided into 2M subintervals. The maximal level of resolution J is defined as M = 2 J . The base Haar wavelet family can be described as The coefficient c i is calculated from The parameter m = 2 j corresponds to the maximum number of rectangular waves that can be sequentially deployed in the interval [A, B] for the given dilation parameter j = 0, 1, . . . , J. The index i is calculated from i = m + k + 1. While the scaling function h 1 (x) = 1 is constant, the other Haar functions contain a single rectangular wave. Since the scaling function h 1 (x) does not include a wave, in its case m = 0, ξ 1 = A, ξ 2 = B, ξ 3 = B. A general example of the Haar wavelets for J = 2 on a nonunfiorm grid is shown in Figure 1. In this figure locations where the vertical lines cross the x axis along with the two boundary points form the grid points.
The Haar functions are orthogonal to each other and therefore form a good transform basis Figure 1. Haar wavelets on a nonuniform grid for J = 2 as described by (3.2)-(3.4).
Thus, any square integrable function f (x) can be expanded into Haar wavelets as where a i are the Haar coefficients. The integrals of order n of the Haar functions (3.1) can be calculated analytically as [58] Within this article 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 (3.7) and (3.8) and considering the vector a = (a 1 , a 2 , . . . , a 2M ) one obtains f (x) = a · H instead of (3.5) and It must be noted, that the matrices H and P n depend on space coordinate x and that the vector a is finite in a discrete setting. It implies from (3.6) that in boundary points A and B hold Due to equation (3.9) the boundary conditions for particular problems are often simplified.

Model equations
The four model equations along with their corresponding exact solutions are described here. The Burgers' equation where ν is the viscosity parameter, was used as the first model equation. Its analytical solution has been shown to be where R = R 0 /(2π) and 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 [2,7,9,27,49]. The KdV equation where α is the nonlinear parameter and β is the dispersion parameter, was used as the second model equation. The analytical solution for the KdV equation has been found in the form where c is the phase speed of the travelling soliton and x 0 denotes the initial phase [67]. The sine-Gordon equation, the third model equation, has been shown to have an analytical solution in the form where c is the speed at which the phase of the travelling wave propagates and x 0 denotes the initial phase [51]. Finally, the mKdV equation where α is the nonlinear coefficient and β is the dispersion coefficient, was used as the fourth model equation. The analytical solution for the mKdV equation has been found in the form where x 0 denotes the initial phase and c is the phase speed of the travelling wave [1,41,82,86].

The higher order Haar wavelet method
This section will give an overview of how the HOHWM is applied for spatial integration in case of each specific model equation within this paper. Firstly, the general steps are described and then the results specific to each model equation are provided.
With the HOHWM, according to [63], the derivative of n + 2s order, where n is the maximal spatial derivative present within the equation, is expanded into the Haar series for a fixed time moment t j . Thus, is obtained. In order to obtain the function u, equation (5.1) is integrated n + 2s times to obtain where c i are unknown coefficients. In equation (5.2), x i denotes an element wise power of the vector x. Since the equation has n spatial derivatives, these provide n conditions that can be evaluated in order to solve for the unknown coefficients c i . However, 2s additional coefficients are introduced by the HOHWM (in comparison to the HWM) which must be obtained using some supplementary information.
In the present paper s = 1 is used. This gives u (n+2)x (x, t j ) = a · H. Integrating n + 2 times yields Using the n boundary conditions as well as evaluating the differential equation at its boundaries (A and B) one obtains the n + 2 equations needed to solve for the n + 2 unknown coefficients c i . Some of these coefficients can depend on the Haar wavelet coefficients vector a. After solving for the unknown coefficients and inserting the results into expression (5.3) one obtains the result in the form of where R n+2 and S n+2 are matrices which only depend on the boundary conditions and the grid. Thus the values for these can be calculated right after the grid has been specified.

The Burgers' equation
For the given homogeneous boundary conditions Since the boundary conditions (5.4) are homogeneous u t (0, t) = u t (1, t) = 0 and (5.5) simplifies to The addition of (5.6) gives The initial condition was taken from the exact solution (4.2). It can be shown that at t = 0 the solution simplifies to which is used as the initial condition. Within the context of this paper u 0 = 1 is used. Since the boundary conditions (5.7) are homogeneous u t (0, t) = u t (1, t) = 0 and the above simplifies to

The KdV equation
The addition of (5.8) gives The initial condition was taken from the exact solution (4.4) at t = 0. This gives which is used as the initial condition with specified values for c, x 0 as well as for α and β.

The sine-Gordon equation
For given constant boundary conditions u(0, t) = 0, u(1, t) = π (5.9) the equation (4.5) is evaluated at the boundary points as Since the boundary conditions (5.9) are constant u t (0, t) = u t (1, t) = 0 and the above simplifies to u xx (0, t) = 0, u xx (1, t) = 0. (5.10) The addition of (5.10) gives The initial condition was taken from the exact solution (4.6). It is evaluated at t = 0 to give which is used as the initial condition with specified value of phase speed c.

The mKdV equation
Given homogeneous boundary conditions u(0, t) = u(1, t) = u x (1, t) = 0, (5.11) the equation (4.7) is evaluated at the boundary points as Since the boundary conditions (5.11) are homogeneous u t (0, t) = u t (1, t) = 0 and the above simplifies to The addition of (5.12) gives The initial condition was taken from the exact solution (4.8) at t = 0. This gives which is used as the initial condition with specified values for coefficients c, x 0 as well as parameters α and β. The sine-Gordon equation (4.5) was numerically solved at J=4 with parameter values varying from d=0.1, 0.11, . . . , 0.85, e=0.05, 0.1, 0.2, 0.3, 0.4, 0.5 and 0.001 < δ < 0.05. The maximum deviation from the exact solution (4.6) was calculated and compared between the various numerical experiments. Such numerical results are shown in Figure 2. Figure 2 demonstrates that the optimal values are located in the region where δ < 0.02 and 0.15 < d < 0.4 for all considered values of parameter e. According to the detailed analysis the best results were obtained for d = 0.33, e = 0.3 and δ = 0.008. With those parameter values the maximum deviation from the exact solution was ∆u = 0.033198. This makes the relative deviation 0.53%. The highest deviation from the exact solution was found on the wavefront of the solution (near x − ct − x 0 = 0). The parameter values between which the best results were obtained for the KdV, the sine-Gordon and the mKdV equations are shown in Table 1.

Discussion of numerical results
The Burgers' equation (4.1) was numerically solved for various values of parameter ν at various resolutions J. The numerical solutions were compared to the exact solution (4.2). The nonuniform grid (2.2) was used for this problem since the solution is known for creating a steep decline near its boundary x = 1. The time moment t f = max(0.5, t c ) was determined by obtaining the value of critical time moment t c which is the maximal time t for which max x,t<tc |u(x, t) − u e (x, t)| < 10 −3 . These numerical results were compared with the uniform grid results previously published in [76]. Such numerical results are shown in Table 2. Table 2.
Values of t f against the resolution J in the case of the Burgers' equation with ν = 1 100π , ν = 1 110π , 1 120π ; parameter q characterises the nonuniformity of the grid. Uniform Nonuniform Uniform Nonuniform Uniform Nonuniform As the resolutions J increases, the value of the nonuniformity parameter q at which the best results were obtained approaches 1. This is because otherwise there would be too many grid points gathered near the boundary of the domain which leads to numerically singular matrices because the grid points or the respective Haar wavelet function values are numerically equal. It is clear that for J > 6 the advantages of a nonuniform grid vanish simply because of the abundance of collocation points. In fact, none of the nonuniform grids used outperformed the uniform grid in this case, which is why such results are not  Table 3 where N = 2M denotes the number of collocation points. It must be noted that for J < 6 the HOHWM using the uniform grid was unable to successfully calculate the numerical result up to the final time of t f = 0.5. However, even at J = 6 the nonuniform grid approach easily outperforms the uniform grid version of the HOHWM. The table also shows that for resolution J > 4 the accuracy no longer increases significantly in the case of the nonuniform grid.
The KdV equation (4.3) was solved at various resolutions J and the numerical results were compared to the exact solution (4.4). The uniform grid results from [76] were used and new results with the nonuniform grid (2.3) and the adaptive grid (2.4) were calculated. The numerical result itself was used as the basis for the weight function in the adaptive grid for this model equation. The maximal deviation from the exact solution for the numerical experiments are shown in Table 4. It must be noted that the values of d, e and δ for which the best numerical accuracy was obtained vary from resolution to resolution.
When dealing with the sine-Gordon equation (4.5), the uniform grid, nonuniform grid (2.3) as well as the adaptive grid (2.4) were used. In this case, the first spatial derivative was used as the basis for the weight function. The numerical results were compared to the exact solution (4.6). The maximal deviation from the exact solution ∆u for the different grids can be seen in Table 5. Finally, the mKdV equation (4.7) was numerically solved using the uniform grid, the nonuniform grid (2.3) as well as the adaptive grid (2.4). The numerical result itself was used for the weight function in the adaptive grid for this equation. The obtained numerical results were compared to the exact solution (4.8). The maximal deviation from the exact solution ∆u for the different grids at different resolutions can be seen in Table 6.  Tables 2-6 show that the nonuniform grid always outperformed the uniform grid version of the HOHWM. This is especially clear given that in some cases the uniform and static nonuniform grid approaches in Tables 4 and 6 are unable to successfully finish the integration. Furthermore, Tables 4-6 show that the adaptive grid outperforms the static nonuniform grid as well. It is clear from Tables 2 and 3 that the same (or even better) accuracy obtained with the uniform at 256 or 512 grid points is achieved with only 32 grid points by employing the nonuniform grid (2.2). Similarly, Tables 4-6 clearly show that the same (or even better) accuracy obtained with the uniform grid at 256 grid points can be achieved with only 64 grid point by employing the adaptive grid (2.4). It must be noted that as can be seen from Table 6, a higher resolution J does not always guarantee a numerical result of higher accuracy. This is the case because the matrices involved get closer to being singular matrices as the resolution increases and solving a linear system with a near-singular matrix introduces inaccuracies.
It is clear that the usage of the nonuniform as well as the adaptive grid will also have an impact on the computational time. The computational times of the nonuniform grid approach is compared to the uniform grid approach in case of the Burgers' equation in Table 7. In Table 7 for J = 8 only the uniform grid approach is considered since a nonuniform approach was not observed to give a better result. Similar comparison is made for the other model equations in Table 8. Tables 7-8 show that a result of similar or better accuracy is able to be obtained with less CPU time using the adaptive grid approach when compared to the uniform grid approach.
An example of how the grid adapts to the moving wave can be seen in Figure 3. The figure shows that as the wave moves, the grid adapts to its position.

Conclusions
A nonuniform grid approach was developed and used alongside the HOHWM in order to account for abrupt changes within the solution. Firstly, two static nonuniform grids were introduced. For either of those grids, a prior knowledge of the shape of the solution is required in order to show significant improvement in accuracy. However, in a lot of cases it is impossible to know the shape of the solution before the start of the calculations. Thus, an adaptive nonuniform grid was proposed. It changes the grid adaptively so that the areas with the most abrupt changes always have the most grid points. The numerical results (Tables 2-3) showed that in case of the Burgers' equation a significant gain in accuracy can be obtained by the use of a static nonuniform grid when compared to the uniform grid HOHWM. In fact, the same level of accuracy was obtained with the nonuniform grid given 8-16 times fewer grid points.
The adaptive grid (Tables 4-6) also showed a significant gain in accuracy when compared to both the uniform grid approach as well as the static nonuni- Table 7: CPU times (in seconds) in case of numerical solutions of the Burgers' equation for cases in Table 2.  form grid approach. The same level of accuracy was obtained using 8 times fewer grid points than were needed for the uniform grid approach.
The most computationally expensive part of the numerical calculation is solving the linear system of equations. Therefore, the ability to obtain results with the same accuracy with considerably fewer grid points is quite an advantage. Instead of solving a 256 × 256 linear system one can solve a 32 × 32 linear system.
The use of a nonuniform or an adaptive grid alongside the HOHWM have shown to give significant advantages. However, further studies need to be carried out in order to find an optimal way of obtaining the necessary parameters for the adaptive grid. Herein constant boundary conditions were utilized. However, time dependent boundary conditions can be handled similarly to cases where differential equations include time dependent function as coefficients or right-hand side terms (similarly to handling of graded materials in [63]). This is a promising subject for future studies.