An efficient algorithm combining an optimized hybrid block method and the differential quadrature method for solving Hunter–Saxton equation

In this paper, the numerical solution of a mixed derivative type Hunter–Saxton equation is addressed. A given equation is discretized transforming it into a system of ODEs with the use of a cubic trigonometric B-splines based differential quadrature method. The system is further solved using a fifth-order optimized one-step hybrid block method. Three numerical illustrations validate the efficiency of the proposed scheme and show its better performance through very accurate results. Stability and convergence analysis are also performed.


Introduction
The non-linear Hunter-Saxton equation (HS) is a PDE which was firstly considered by Hunter and Saxton [5], given by where t and x are the time and space variables, and u corresponds to splay waves that propagate throughout liquid crystal molecules. This equation describes the orientation waves that cause disclinations throughout the crystal structure. The HS equation modeled as high frequency limit of the shallow water equation, the Camassa-Holm equation, models soliton as well as breaking waves [11]. This equation describes the geodesic flow on the diffeomorphism group of the unit circle [17]. Different real-life applications of this equation entice researchers to find its analytical and numerical solutions. Numerous numerical methods available in the literature have been applied to the HS equation. Rouhparvar [20] used a semi-analytical differential transform method. A finite difference scheme has been performed by Holden [13], and Arbabi [4] applied a Haar wavelet quasilinearization approach along with a finite difference scheme. Behzadi [7] provided approximate solutions using the Adomian decomposition method, while Aratyn [6] used the Padé approximation technique. Other schemes using collocation based on cubic trigonometric B-splines [12], quintic B-splines [14], Laguerre wavelets [24] and Barycentric interpolation [26] have been implemented to give approximate solutions of the HS equation.
A well-known simple and cost-efficient approach considered by researchers is the differential quadrature method (DQM). DQM is a numerical technique used to approximate derivatives of the solution and was firstly introduced by Bellman et al. [8]. Later, Quan and Chang [18] improved the DQM by introducing weighting coefficients. A generalized formula for higher order derivatives was given by Shu [22]. To find better approximations, several types of trial functions, like Lagrange interpolation polynomial, B-spline basis functions and many others have been considered [9,15,16,25]. To the best of our knowledge, DQM has never been applied to the HS equation. Moreover, to overcome the numerical instability, DQM based on cubic trigonometric B-splines can be a rescuer in this situation.
Along with DQM, several methods, for instance, finite difference schemes or multi-stage methods such as SSP-RK can be used for the temporal dimension. A major shortcoming of these schemes is the requirement for a finer grid size to attain acceptable accuracy and stability. To overcome the limitations of the above-mentioned schemes, hybrid block methods come into the scene. Hybrid methods assess the information at off-step points and usually leads to zero-stability and better stability characteristics. Additionally, hybrid methods have the ability to change the step size and violate the Dahlquist's barrier [10]. Apart of the property of being self-starting, hybrid block methods have been developed which claimed to be A-stable. Therefore, these methods are qualified to deal with stiff-problems [23] and are cost-efficient [19] compared to other methods. Table 1 Values of TB l (x) and its derivatives at the nodal points Hybrid block methods have recently been used to solve non-linear parabolic PDEs in [2,3]. The behavior of such methods has not been studied with mixed derivatives yet. Furthermore, to study the effect of the proposed algorithm, the results obtained were compared with exact solutions and gave better accuracy using a considerably small number of grid points in comparison to existing methods in the literature.
Further, cubic TB-basis functions are modified as [25] MTB Now, assuming that the function u(x, t) is sufficiently smooth on its entire solution domain, the approximate values of the derivatives with respect to x at x i can be approximated by where q denotes the order of the derivative and a (q) i j represent the weighting coefficients for the q-th order derivative.

Evaluation of weighting coefficients
To determine the weighting coefficients, consider the MTBs as the test functions in the DQM Eq. (3) using q = 1 for the approximation of first order derivatives, which gives This results in an algebraic system of equations that can be written in matrix form as , and the columns of S are By using MATLAB 2017b, we solve each tridiagonal system Now, considering Shu's recurrence formula [21], the coefficients in (3) for approximating the second order derivatives can be evaluated as follows Hence, we have obtained all weighting coefficients and using them, approximations to the partial derivatives can be obtained using the formula in (3). Note that the Hunter-Saxton Eq. (1) can also be written as and initial condition Initially, the DQM is applied to the Hunter-Saxton Eq. (6) using the approximations given in equation (3). This transforms the partial differential Eq. into an initial value problem system (IVP) given by: for i = 1, 2, . . . , N − 1 and, where the boundary conditions are incorporated through To solve the obtained IVP system of ODEs, we will consider a block method, which is described in the next section.

Hybrid block method for the temporal variable
Here, we present the time-stepping scheme considered for the HS equation, which is a sixth-order hybrid block method. Using the concept of interpolation and collocation, an optimized one-step hybrid block method for solving a first order IVP of the form {y = f (t, y); y(t 0 ) = y 0 } has been derived in [23], which was proven to be convergent and A-stable. The four equations of the block method in [23] to be applied on an interval [t i , t i+1 ], with k = t i+1 − t i , are as follows: where μ 1 = 3 6 − √ 3 6 , μ 2 = 1 2 and μ 3 = 3 6 + √ 3 6 . The above method comprising four equations, is further applied to the system of ODEs (7) and the resulting system of equations on each iteration is solved using the "fsolve" function of the Matlab software. Figure 1 shows the flowchart for the algorithm DQM-OHB. The hybrid block method in [23] is convergent and A-stable. If the system of ODEs obtained after the use of the DQM is stable, then the algorithm DQM-OHB will be unconditionally stable. In order to prove the stability of the non-linear system (7), linearize the non-linear terms by assuming u to be locally constant [1]. Therefore, we get the following linear system of equations Now, rewriting it in matrix form, one gets . . , du N dt and F consists of the nonhomogeneous part, which is related to the boundary conditions. If the real part of all the eigenvalues of the coefficient matrix have negative values then the system of ODEs is stable. We have plotted the eigenvalues of the coefficient matrix (K ) −1 L of the linear system (8) for different values of N in Figure 2. Since all the eigenvalues have negative real part, this implies that the linearized system of ODEs is stable. As there is no restriction concerning the grid size, therefore, the algorithm DQM-OHB is unconditionally stable.

Numerical results
In this section, three test problems have been considered to validate the efficiency of the presented method. The numerical results were obtained using Matlab 2017b on a personal computer with Intel Core i5-4210U processor at 2.40 GHz. For error determination, the following formula is used:  , t) where t will specified on each case. The formula to calculate the approximate order of convergence in terms of the L ∞ errors with respect to time is given below: where k i , k i+1 are two temporal step sizes. Test problem 1 Consider the exact solution for the Hunter-Saxton Eq. (6) as follows: By extracting the initial condition and boundary conditions, the problem is solved using DQM-OHB taking N = 79. Table 2 lists the absolute error at t = 1 taking k = 0.1. The approximate orders of convergence are given in Table 3. We can see that the orders of convergence are close to six, as expected. The 3D-plot of the obtained numerical solution is shown in Fig. 3.

Test problem 2
We assume the following initial condition for Eq. (6):    and boundary conditions to be Its exact solution is as follows: Comparison of the results with the numerical methods in [4,12] and that of the algorithm DQM-OHB for t = 0.1 is shown in Table 4 considering N = 128 taking five time-steps. The results show the advantage of higher accuracy of DQM over other collocation methods. Table 5 displays the absolute errors for t = 0.01 taking five timesteps and Table 6 shows the results for t = 0.001 taking one time-step. These results demonstrate the efficiency and accuracy of the proposed algorithm over existing ones. Figure 5 represents a 3-D plot generated by DQM-OHB by taking N = 128 and k = 0.1 up to time t = 1. Convergence analysis performed for DQM-OHB is shown in Fig.  4 for different number of steps. Table 7 shows L ∞ errors computed for larger values of t. Test problem 3 Consider the Eq. (6), where x ∈ [1,2] and assuming the initial condition as       2 .
In Table 8, absolute errors for t = 0.1 are tabulated by taking k = 0.1 and N = 127 for DQM-OHB. The absolute errors calculated for t = 0.01 and t = 0.001 are given in Table 9 and Table 10, respectively. Each of the tables lists the results obtained in a single time-step.

Conclusion
Numerical experimentations of the presented algorithm for the Hunter-Saxton equation highlight the efficiency of DQM over collocation methods performed so far in the literature. The order of convergence with respect to time for DQM-OHB is six and unconditional stability is proved. To solve non-linear mixed derivative PDEs, the proposed method does not require linearization approaches or any initial guesses and hence, is a direct method. DQM-OHB can obtain higher temporal as well as spatial resolution. This work extends the scope of hybrid block methods in the field of partial differential equations.
Author Contributions A.K. wrote the draft manuscript text and prepared figures. V.K. and H.R. addressed the formal analysis and supervised the manuscript. All authors reviewed the manuscript.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature. Anurag Kaur is supported financially by the funding agency, the University Grants Commission (UGC), New Delhi, India under the scheme of UGC-CSIR NET-SRF with reference id 403645.
Data availibility Not applicable.

Declarations
Competing interests The authors declare that they have no conflict of interest.
Ethical Approval Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.