# A Mesh-Free Collocation Method Based on RBFs for the Numerical Solution of Hunter–Saxton and Gardner Equations

Authors: Rahman Ullah; Ihteram Ali; Sadaf Shaheen; Tahir Khan; Faizullah Faiz; Hamid Rahman
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, the radial basis function (RBF) collocation method is applied to solve nonlinear partial differential equations (PDEs). First, the given equation is reduced to time-discrete form using Ө-weighted scheme. Then, with the help of RBFs, the given PDEs are transformed into a system of algebraic equations that is easy to solve. The proposed technique is verified by solving Hunter–Saxton and Gardner equations. For the solution of these problems, three types of RBFs, including multiquadric (MQ), inverse multiquadric (IMQ), and Gaussian (GA) have been used. The obtained solution is validated with the help of absolute error,L2, and L∞ error norms. The involved shape parameter is selected by the automatic algorithm that lies in the stable region of the method. The stability of the scheme is discussed by the spectral matrix method and validated computationally.


---

## Body

## 1. Introduction

Partial differential equations (PDEs) are used to model various physical, chemical, and biological processes. Due to their considerable application, researchers are continuously looking for the robust numerical method to find their approximate solution. Consider the one-dimensional nonlinear Hunter–Saxton (HS) equation given as follows:(1)wvvt+2wvwvv+wwvvv=0,which was modeled by Hunter and Saxton [1] to describe the wave propagation in the director field of nematic liquid crystal, where v is a variable with the moving frame of reference and t is the slow time variable. Many researchers have used various numerical techniques for approximate solutions of the Hunter–Saxton equation. For instance, Baxter et al. [2] studied the model and found separable as well as self-similar solutions to the problem. Behzadi [3] also studied the equation and obtained series solutions through analytical methods. Arbabi et al. [4] computed semianalytical solution by using Haars and Adomian decomposition method.Also, we consider the Gardner (KdV + mKdV) equation which was developed by Clifford Gardner in 1968 and is of the form as follows:(2)wt−η1wwv+η2w2wv+η3wvvv=0,where η1,η2,andη3 are arbitrary constant. This model is used to study the propagation of negative ion-acoustic plasma waves and also explain internal waves with large amplitude [5]. Various techniques determine the exact solution of the Gardner equation. Zuntao et al. [6] found solitary wave solution using projective Riccati equation. Wazwaz [7] obtained the kink solution by applying the tanh approach. The approximate solution of the Gardner equation is computed using various numerical techniques. The propagation of one soliton solution was approximated by a conservative finite difference scheme in [8]. Restrictive Taylor’s technique and extended cubic spline approaches were also applied to find approximate solution [9].In the present work, we use the mesh-free Radial basis functions (RBF) method to investigate the approximate solution of the above models. The RBF meshless method is the convenient and most powerful technique to solve multivariate problems. Owing to its fast convergence, ease implementation and low computational cost RBFs technique have given preference over the traditional methods. In 1971, Hardy used multiquadrics to interpolate two-dimensional data and obtain accurate results. In 1982, Franke [10] conducted a number of experiments using different techniques for scattered data approximation and concluded that MQ-RBF is the best one. Kansa (1990) extended the study by using the radial basis function method to find a numerical solution of partial differential equation [11]. Later on, most of the researchers have used the RBF method to solve different kinds of PDEs, e.g., Korteweg–de Vries equation [12], Coupled Korteweg–de Vries equations [13], Kuramoto–Sivashinsky equation [14], RLW equation [15], and Burger equation [16]. Most recently, researchers are using the local meshless method for approximate solution of various classes of differential equations. In [17, 18], Pouria Assari and his co-worker investigated the meshless local discrete collocation and discrete Galerkin method for boundary integral equations. The same author studied the solution of logarithmic singular boundary integral equations which arise from the Laplace equation using meshless Galerkin scheme based on radial basis function [19, 20]. Similarly, other numerical techniques are used for approximate solution of different class of PDEs in [21–24] and the reference there in.The convergence of the method mostly depends on shape parameterc, which needs to be defined by the user. The importance of shape parameter c in RBF was discussed in detail by Tarwater [25]. Many algorithms are available in the literature for the selection of c. Hardy (1971) introduced an algorithm c=0.815d where d=1/N∑di. Fasshauer (2002) chose optimal parameter as c=2/N. Similarly, Rippa [26] introduced the other way for selection of optimal c. However, all of the above tools have no generality in selecting shape parameters. In the current discussion, we used an automatic algorithm for selection of shape parameter c, which choose value of c in the region of stability which corresponds to the convergence of the solution. For computation, we use MATLAB 2013, using Intel core-i7 computer having 4 GB Ram. Rest of the paper is organized as follows. In Section 2, general scheme, stability of the method, and quasi-linearization are discussed, while Section 3 presents error norms and some test problems. Finally, the paper is concluded in Section 4.


## 2. Interpretation of the Proposed Method

To explain the basic concept of the method, consider a general one-dimensional time-dependent PDE over a spatial domainΓ and time t∈0,T.(3)∂w∂tv,t+£wv,t=hv,tv∈Γ,t∈0,T.With initial and boundary conditions,(4)wv,0=wov,v∈Γ,ßwv,t=gv,t,v∈∂Γ,where \pounds is any differential operator in spatial domain v∈Γ, ß is boundary operator, hv,t is given smooth function, and Γ and ∂Γ are interior and boundary of the spatial domain. Applying θ−weighted scheme, equation (3) can be written as(5)1δtwn+1v−wnv+θ£wn+1v+1−θ£wnv=hv,tn+1.Now, approximate the spatial part of (3) by RBFs as(6)wnv=∑j=1N⋋jnΥrijrij=vi−vj,v∈Γandi=1,2,….,N.where . is the Euclidean norm between collocation points. The solution wv,tn is represented by wnv which can be written in matrix form as(7)wn=A⋋n,where ⋋n are the unknown coefficient vector which can be determined and A is N×N matrix whose elements are in the form Aij=Υvi−vj. Now, using (6) in (5), we have(8)1δt∑j=1N⋋jn+1Υrij−∑j=1N⋋jnΥrij+θ∑j=1N⋋jn+1£Υrij+1−θ∑j=1N⋋jn£Υrij=hv,tn+1.After rearranging, we have(9)∑j=1N⋋jn+1Υrij+δtθ∑j=1N⋋jn+1£Υrij=∑j=1N⋋jnΥrij+δtθ−1∑j=1N⋋jn£Υrij+δthv,tn+1.Similarly, the boundary condition can be approximated as follows:(10)∑j=1N⋋jn+1ßΥrij=gv,tn+1.The matrix form of the above equations can be written as follows:(11)H⋋n+1=F⋋n+Gn+1.Forj=1,…,N, the components of these matrices are as follows:(12)H=Υrij+δtθ£Υrij,i=2,…,N−1ßΥrij,i=1,N,(13)F=Υrij+δt1−θ£Υrij,i=2,…,N−10,i=1,N,(14)G=δthn+1v,i=2,…,N−1gn+1v,i=1,N,where H and F are N×N matrices. Gand⋋ are N×1 column vectors. Using the above matrices in (11), the unknown coefficient ⋋ can be found easily, and then solution can be obtained via (7).


### 2.1. Stability Analysis

To check the stability of system (11), we use a spectral radius of the amplification matrix. Let w denote approximate solution while u represents the exact solution, and then the error can be defined as follows:(15)En=un−wn.From (7), we know wn=Aλn⇒λn=A−1wn; by putting the value of λ for any n in (11), we have(16)HA−1wn+1=FA−1wn+Gn+1.After rearranging, we get(17)wn+1=AH−1FA−1wn+AH−1Gn+1.The above difference equation gives approximate solution at any time level′n′ for all nodal points. Similarly, if u is the exact solution of PDE, then it must satisfy the difference equation such that [27](18)un+1=AH−1FA−1un+AH−1Gn+1.Substituting the values from (11) and (18) in (15), we have(19)En+1=un+1−wn+1=AH−1FA−1En=MEn,where M=AH−1FA−1 is the amplification matrix. Scheme (17) is stable when the spectral radius ρM of matrix M is such that ρM≤1, where ρM=maxκii=1N and κ are the eigenvalue of matrix M. The validation of this condition is described with the help of the following example.


### 2.2. Quasi-linearization

Quasi-linearization technique is generalization of the Newton–Raphson method for functional equations. It converges quadratically. Consider a nonlinearmth order differential equation as(20)ψmv=gv,ψ,ψ′,ψ″,…,ψm−1,wherev∈Γ.Using quasi-linearization technique, (20) reduces to following form:(21)ψs+1m=gv,ψ,ψ′,ψ″,…,ψm−1+∑i=0m−1ψs+1i−ψsigψiv,ψs,ψs′,ψs″,…,ψsm−1,which is mth order linear differential equation in iterated form and s denotes number of iterations. From equation (21), one can easily compute ψ at s+1th iteration when it is known at sth iteration. For better understanding, we consider a second-order nonlinear differential equation of the form as follows:(22)ψ″v=gv,ψ,ψ′,forv∈Γ,where ψ′ is another function, and then by using (21), we linearize (22) as(23)ψs+1″v=gv,ψs,ψs′+ψs+1−ψsgψv,ψs,ψs′+ψs+1′−ψs′gψ′v,ψs,ψs′.Equation (23) is always a linear differential equation and can be solved recursively, where ψsv is known and one can use for obtaining ψs+1v.


## 2.1. Stability Analysis

To check the stability of system (11), we use a spectral radius of the amplification matrix. Let w denote approximate solution while u represents the exact solution, and then the error can be defined as follows:(15)En=un−wn.From (7), we know wn=Aλn⇒λn=A−1wn; by putting the value of λ for any n in (11), we have(16)HA−1wn+1=FA−1wn+Gn+1.After rearranging, we get(17)wn+1=AH−1FA−1wn+AH−1Gn+1.The above difference equation gives approximate solution at any time level′n′ for all nodal points. Similarly, if u is the exact solution of PDE, then it must satisfy the difference equation such that [27](18)un+1=AH−1FA−1un+AH−1Gn+1.Substituting the values from (11) and (18) in (15), we have(19)En+1=un+1−wn+1=AH−1FA−1En=MEn,where M=AH−1FA−1 is the amplification matrix. Scheme (17) is stable when the spectral radius ρM of matrix M is such that ρM≤1, where ρM=maxκii=1N and κ are the eigenvalue of matrix M. The validation of this condition is described with the help of the following example.


## 2.2. Quasi-linearization

Quasi-linearization technique is generalization of the Newton–Raphson method for functional equations. It converges quadratically. Consider a nonlinearmth order differential equation as(20)ψmv=gv,ψ,ψ′,ψ″,…,ψm−1,wherev∈Γ.Using quasi-linearization technique, (20) reduces to following form:(21)ψs+1m=gv,ψ,ψ′,ψ″,…,ψm−1+∑i=0m−1ψs+1i−ψsigψiv,ψs,ψs′,ψs″,…,ψsm−1,which is mth order linear differential equation in iterated form and s denotes number of iterations. From equation (21), one can easily compute ψ at s+1th iteration when it is known at sth iteration. For better understanding, we consider a second-order nonlinear differential equation of the form as follows:(22)ψ″v=gv,ψ,ψ′,forv∈Γ,where ψ′ is another function, and then by using (21), we linearize (22) as(23)ψs+1″v=gv,ψs,ψs′+ψs+1−ψsgψv,ψs,ψs′+ψs+1′−ψs′gψ′v,ψs,ψs′.Equation (23) is always a linear differential equation and can be solved recursively, where ψsv is known and one can use for obtaining ψs+1v.


## 3. Numerical Experiment

In this section, the proposed method is applied to find the approximate solution of the above models and compare the results with exact and existing numerical solutions for different time levelstj and collocation points vj with specific spatial and time step size δvandδt, respectively. The performance of the method can be examined in terms of L2 and L∞ and root mean square Lrms norms, which are defined as follows:(24)L2=δv∑j=1NUj−wj21/2,L∞=max1≤j≤NUj−wj,(25)Lrms=1/N∑j=1NUj−wj2,where U is the exact and w is the approximate solution of the considered problems.Example 1.
Consider the Hunter–Saxton equation [4].(26)wvvt+2wvwvv+wwvvv=0,with initial condition wv,0=2v and boundary conditions w0,t=0,w1,t=2/1+t.
The exact solution of the problem is(27)wv,t=2v1+t.
By comparing (26) with (3), we have(28)∂w∂t↦∂wvv∂t,and£=2wvv∂∂v+wvvv.
Using the methodology discussed in Section2, (26) reduces to(29)wvvn+1+δtθ2wvwvvn+1+wwvvvn+1=wvvn+δtθ−12wvwvvn+wwvvvn.
Apply quasi-lineariation technique to linearize the nonlinear terms as(30)wvwvvn+1=wvwvvn+wvn+1−wvnwvwvvwvn+wvvn+1−wvvnwvwvvwvvn,=wvwvvn+wvn+1−wvnwvvn+wvvn+1−wvvnwvn.
After simplification, we get(31)wvwvvn+1=wvn+1wvvn+wvnwvvn+1−wvwvvn.
Using similar procedure forwwvvv, we have(32)wwvvvn+1=wn+1wvvvn+wnwvvvn+1−wwvvvn.
Now, putting the above linear terms in (29), we have(33)wvvn+1+δtθ2wvn+1wvvn+wvnwvvn+1+wn+1wvvvn+wnwvvvn+1=wvvn+δtθ−12wvwvvn+wwvvvn+δtθ2wvwvvn+wwvvvn.
After using RBFs approximation for the function and their derivatives, the elements of the matrices in (11) can be obtained from the following:(34)H=Υ′′rij+δtθ2wvvnΥ′rij+wvnΥ′′rij+wvvvnΥrij+wnΥ′′′rij,i=2,…,N−1Υrij,i=1,N,(35)F=wvvn+δtθ−12wvwvvn+wwvvvn+δtθwvwvvn+wwvvvn,i=2,…,N−10,i=1,N,(36)G=0,i=2,…,N−1gn+1v,i=1,N.
For comparison, the solution of the problem is calculated using three different RBFs, which are Multiquadric (MQ), Inverse multiquadric (IMQ), and Gaussian (GA). The values of corresponding shape parameterc=1e3,390,1e−3 are, respectively, selected experimentally. The solution is computed over the domain v∈0,1 for the collocation points N=128, time T=1, and time step dt=0.001. The performance of the method is examined in terms of absolute error. The computed results by three RBFs are compared with the results produced by the Haar wavelet method [4] which are presented in Table 1. From the table, it is clear that the RBFs approximation is more accurate and reliable than [4]. The table also indicates that the absolute error produced by Gaussian is better than MQ, whereas MQ is better than IMQ RBFs. In Figure 1, surface plots of exact and approximate solutions are displayed, whereas absolute errors obtained by MQ and Gaussian are shown in Figure 2. The shape parameters versus ρM is plotted in Figure 3, which shows that the method should be stable when the value of c is selected from the stable region.Table 1
Comparison of absolute errors atT=1 and N=128 for Example 1.
iMQIMQGaussian[4]11.984e−089.406e−081.446e−087.36e−0939.762e−091.667e−081.408e−085.21e−0853.011e−086.351e−085.753e−086.91e−0874.017e−088.968e−087.116e−081.10e−0795.023e−081.773e−079.967e−081.32e−07595.077e−071.541e−066.888e−078.88e−07615.176e−071.543e−067.170e−079.16e−07635.275e−071.609e−067.153e−079.46e−07655.671e−071.588e−067.286e−079.76e−07675.471e−071.629e−067.269e−071.01e−06695.570e−071.647e−067.550e−071.04e−061192.040e−075.577e−072.759e−071.79e−061211.242e−073.976e−072.142e−071.82e−061231.041e−073.027e−071.525e−071.85e−061255.411e−081.789e−071.206e−071.88e−061274.108e−093.095e−085.896e−081.91e−06Figure 1
Exact and approximate solution profile atT = 1 of Example 1.
(a)(b)Figure 2
Error plots of Example1.
(a)(b)Figure 3
Spectral radius vs. shape parameter for IMQ (a) and GA (b).
(a)(b)Example 2.
Consider the nonlinear problem [3].(37)wt+2wv−wvvt+wwv−2wvwvv−wwvvv=0,with initial condition(38)wv,0=expv,and corresponding boundary conditions(39)w0,t=exp−3t,w0.5,t=exp0.5−3t.
The exact solution to the problem is given by(40)wv,t=expv−3t.
The approximate solution of this problem is calculated using MQ RBFs with the value ofc=0.3 at 0≤v≤0.5. The computed solutions are compared with the exact solutions at T=0.01, 0.1, 1, v=0.02, 0.16, 0.36, 0.48, δt=0.001, and N=50. The results are shown in Table 2 along with absolute errors. The table clarifies that RBFs approximation gives reasonable accuracy while using less number of collocation points. Also, it is clear from the table that as the final time T increases, the error decreases. Figure 4 presents exact, numerical, and absolute error plots for T=0.1 which shows a good agreement between exact and numerical solutions.Table 2
Absolute errors of Example2 when N = 50.
TvExactApproximateError0.010.020.990040.990013.0709e−050.161.138821.138875.0677e−050.361.390961.390781.8553e−040.481.568311.568332.3003e−050.10.020.755780.755762.1255e−050.160.869350.869393.3433e−050.361.061831.061691.3860e−040.481.197211.197221.2046e−0510.020.050790.050783.8145e−060.160.058420.058411.2941e−050.360.071360.071332.3126e−050.480.080450.080454.7831e−06Figure 4
Exact, approximate solutions, absolute errors, and spectral radius graphs for Example2.
(a)(b)(c)(d)Example 3.
Consider the general Gardner (2) given as(41)wt−η1wwv+η2w2wv+η3wvvv=0,ηi>0,i=1,2,3.
Forη1=1,η2=−5, and η3=1 in (41), we have(42)wt+wwv−5w2wv+wvvv=0,v∈−80,80,andt>0.
The exact solution [28] is(43)wv,t=1101−tanh3060v−t30.
The initial and boundary condition is extracted from the exact solution. To solve this problem, we implemented MQ, IMQ, and GA RBFs with the value of shape parameterc=0.01,1e−8,20. The solution is computed for T=0.1,1,2 with δt=0.01 and collocation point N=160. The obtained results in terms of L2, L∞, and Lrms error norms along with the spectral radius ρM of the amplification matrix are presented in Table 3. Also, the CPU time computed in seconds is shown in the same table. The table shows that the RBF collocation method is applicable even in a large domain and gives a stable and accurate solution. In Figure 5, exact and numerical solutions are shown. In Figure 6, absolute errors obtained by MQ and IMQ are plotted which indicate that the absolute error produced by MQ is better than that by IMQ. The graph of spectral radius vs. shape parameter for MQ and GA is shown in Figure 7.Table 3
Computed error norms whenN = 160 of Example 3.
RBFsTL2L∞LrmsρMCPUMQc=0.010.11.3323e−043.4912e−051.0500e−051.0000050.0453911.3954e−043.6556e−051.0997e−051.0000050.1097522.7977e−037.3277e−042.2049e−041.0000050.30672IMQc=1e−80.11.1047e−042.8906e−058.7063e−061.0000010.0251911.1570e−033.0275e−049.1187e−051.0000010.0915722.3198e−036.0702e−041.8283e−041.0000010.27502GAc=200.11.0421e−042.6589e−058.2129e−061.0000020.0190511.0916e−032.7848e−048.6030e−041.0000020.1906222.1889e−035.5835e−041.7251e−041.0000020.28174Figure 5
Solution profile of Example3 at T = 2.
(a)(b)(c)Figure 6
Absolute error plots for MQ and IMQ RBFs for Example3.
(a)(b)Figure 7
Spectral radius vs. shape parameter graphs for MQ and IMQ RBFs of Example3.
(a)(b)Example 4.
Considering the value ofη1=4,η2=−3, and η3=1 in (41), we have(44)wt+4wwv−3w2wv+wvvv=0,t≥0.
The exact solution in this case is given in [7] as(45)wv,t=212+314cosh−v/3+5/3+t/27.
The initial and boundary condition is obtained from the exact solution. This problem is solved in the spatial domain−20≤v≤30 at time T=1. The results were produced by two different RBFs like MQ and IMQ with c=0.01,1e−6, respectively. In Table 4, we compare exact and numerical solution when time step δt=0.01 and nodal points N=100, where i represents index of nodal points vi=−20+i−1dv with step size dv=0.5. From the table, it is clear that the solution produced by the RBFs method is found in good agreement with the exact solution. The error in both cases is comparable with each other. In Figure 8, exact and numerical solutions are shown. In Figure 9, absolute errors obtained by MQ and IMQ are plotted. The figure indicates that the absolute error produced by MQ is better than IMQ. Spectral radius vs. shape parameters for MQ and IMQ RBFs is plotted in Figure 10.Table 4
Comparison of absolute errors for Example4.
MQIMQiExactApproximateErrorApproximateError30.0001150.0001194.349e−060.0001194.343e−06130.0006070.0006302.294e−050.0006302.285e−05210.0022810.0023678.566e−050.0023668.493e−05430.0563040.0575071.202e−030.0575241.219e−03510.0860850.0857363.487e−040.0861324.651e−05630.0379130.0369389.755e−040.0368921.021e−03710.0122080.011813.940e−040.0118024.069e−04830.0017660.001706.248e−050.0017036.289e−05990.0001230.0001194.436e−060.0001194.456e−06Figure 8
Exact and approximate solutions for Example4 at T = 1.
(a)(b)(c)Figure 9
Absolute errors for MQ and IMQ.
(a)(b)Figure 10
Spectral radius vs. shape parameter of MQ and IMQ RBFs.
(a)(b)


## 4. Conclusion

A mesh-free collocation method based on radial basis functions (RBFs) is applied to approximate time-dependent PDEs. The approximate solutions of the Hunter–Saxton and Gardner equations have been computed using MQ, IMQ, and GA RBFs. The involved shape parameter was selected using an automatic algorithm. The method’s stability is discussed using the matrix method and verified computationally. Quasi-linearization technique is used for linearizing the nonlinear term involving in the models. Performance of the method is examined with the help of various error norms. The computed results are compared with other available techniques in the literature. The tabulated results clarify that the proposed method gives more accurate solution than existing techniques.



---

*Source: 2152565-2022-05-18.xml*