Linear barycentric rational collocation method for solving a class of generalized Boussinesq equations

1.   Introduction

In this article, we investigate numerical solutions of a class of generalized Boussinesq equations:

where $(\cdot)_x$ is the derivative for $x$, $s$ is a positive integer, $\beta_1$, $\beta_2$, $\beta_3$ are constants, and $0 < \gamma\leq 1$ is a small parameter. For solving the above equations, we take the following initial condition:

and the following boundary condition:

Then, we obtain three $(1+1)$ dimensional nonlinear PDEs. It is remarked that other forms of boundary conditions may be studied in an identical way.

The Boussinesq equations have various kinds of applications in physics, such as describing propagation of long waves in shallow water, vibrations in a nonlinear string, nonlinear lattice waves, iron sound waves in plasma, see ^[1,2,3,4,5] and references therein. In ^[4,5], Daripa et al. discussed some physical implications of (1.3) under the background of water waves. It is remarked that if $\gamma = 0$ in (1.3), then (1.3) becomes the standard Bousinesq equation (1.2).

As we know, very complicated behaviors exist in nonlinear differential equations, such as oscillation, bifurcation, and chaos, see ^[6,7,8]. It is natural to solve these equations to obtain their solutions to study their dynamical behaviors. However, most of nonlinear equations, especially nonlinear PDEs, can not be solved. Therefore, various of numerical methods for solving nonlinear equations appeared, for example, the finite element method, the finite difference method, the boundary element method, and the collocation method, see ^[9,10,11,12]. It is noted that the former three methods are classical numerical computation methods, which treat the time and space variables separately when solving PDEs. The last one is very simple but efficient numerical computation method, especially with fits for nonlinear problems. This method can treat the time and space variables equally and simultaneously, and establish collocation computation schemes for the initial and boundary problems of PDEs, see ^[9,10,11,12].

In recent years, the barycentric rational interpolation ^{[13,14,15,16]} has been proposed by researchers for the collocation method. Specifically, Floater ^[17,18] presented one linear rational interpolation formula and showed the convergence rate for equidistant partition. Wang et al. ^{[9,10,19,20,21]} used the barycentric interpolation method to study various kinds of linear and nonlinear problems. Luo et al. ^[22] has used the barycentric rational collocation method (BRCM) to solve the nonlinear parabolic PDEs. More recently, Li et al. has used the linear barycentric rational collocation method (LBRCM) to solve Volteria-differential equations ^[23,24], biharmonic equations ^[25], heat conduction and diffusion equations ^[26,27], telegraph equations ^[28], and so on. From the above results, we see that the LBRCM is very effective in solving nonlinear problems. This motivated us to study numerical solutions of the nonlinear PDEs (1.1)–(1.3) by using the LBRCM. In addition, Akinyemi et al. ^[29] also studied the numerical solutions of PDEs (1.1)–(1.3), where they used the homotopy perturbation technique method (HPTM), which consisted of the homotopy perturbation method and the Laplace method. The error between numerical and analytic solutions in their examples had very high accuracy. This showed that the HPTM is very effective. However, the convergence rate of the errors was not given. Besides, it used a power series expansion and had to compare the same power to get the coefficients of the power series, which was very complicated and even difficult to obtain. In this paper, we will apply the LBRCM to solve nonlinear PDEs (1.1)–(1.3). By using the method of direct linearization, nonlinear PDEs will be transformed into linear PDEs whose forms can be expressed as simple matrix equations. By using the LBRCM, the numerical solutions will be obtained. Furthermore, by using an error estimate for the barycentric interpolation, the convergence rate of the errors will be verified. This shows that the LBRCM does not require very complicated computations and is very simple and effective.

The paper is arranged as follows. The direct linearization of nonlinear PDEs and differentiation matrices of the corresponding discrete forms are given in Section 2. The error estimate is given in Section 3. Three examples are presented to show validity of theoretical results in Section 4. Finally, conclusions are made in Section 5.

2.   Direct linearization and differentiation matrices

We first give the linearization equations. Let $u_0$ be a known initial function. Then, the generalized Boussinesq equations (1.1)–(1.3) are transformed into linear PDEs on $\Omega: = [a, b]\times[0, T]$ by using the direct linearization method ^[10]:

Therefore, we can obtain the iterative formats of the above equations on $\Omega$ as the following:

where $p = 1, 2, \cdots$. In practical computation, for a given control precision $\varepsilon$, when $\|u_p(x, t)-u_{p-1}(x, t)\|\leq\varepsilon$, iterations will be stopped and an approximate solution for a nonlinear problem will be derived. This process is called the direct linearization iteration method.

Next, we will apply the barycentric rational interpolation method to derive the differentiation matrices of (2.4)–(2.6). Two forms of the partition for the domain $\Omega$ will be used in this paper. One is the equidistant partition. Let $x_i = a+hi$, $t_j = \tau j$ with $h = \frac{b-a}{m}$, $\tau = \frac{T}{n}$ for $0\leq i\leq m$, $0\leq j\leq n$. Then, $(x_i, t_j)$ divides $\Omega$ into a uniform mesh. The other is the non-equidistant partition. Here, the second kind of Chebyshev points $x_i = \cos\frac{i\pi}{m}$, $t_j = \cos\frac{j\pi}{n}$ for $0\leq i\leq m$, $0\leq j\leq n$, are used to form a nonuniform mesh. In this case, set $h = \max\limits_{0\leq i\leq m-1}(x_{i+1}-x_i)$ and $\tau = \max\limits_{0\leq j\leq n-1}(t_{i+1}-t_i)$.

Let $u_{mn}(x, t)$ in the following be the approximate function of $u(x, t)$,

where $u_{ij} = u(x_i, t_j)$,

and

$Q_{1i} = \{q\in Q_1: i-d_1\leq q\leq i\}$ for $Q_1 = \{0, 1, \cdots, m-d_1\}$ and $Q_{2j} = \{q\in Q_2: j-d_2\leq q\leq j\}$ for $Q_2 = \{0, 1, \cdots, n-d_2\}$. With a direct calculation, one can obtain another form of $u_{mn}(x, t)$ as follows:

where

The expression of (2.12) has a barycentric interpolation form and is very suitable for the program implementation in Matlab. It is obvious that $\xi_i(x_i) = 1$, $\xi_k(x_i) = 0$ for $k\neq i$, and $\eta_j(t_j) = 1$, $\eta_k(t_j) = 0$ for $k\neq j$. By computing the derivative of the basis function (2.8) at $x = x_i$, we can get the following results:

For convenience, denote the $p$-order derivative of (2.8) at $x_i$ as follows:

where $p$ is some positive integer. Then, by induction, one can get recurrence formulas for the $p$-order derivative of (2.8) at $x_i$

Similarly, if we take the derivative of the basis function (2.9) at $t_j$, we get

We also denote the $p$-order derivative of (2.9) at $t_j$ as the following:

for some positive integer $p$.

Let

and

be the matrices of order $(m+1)$ and $(n+1)$, respectively. Before introducing differentiation matrices, we need to give some other notations. Let $x = (x_1, x_2, \cdots, x_n)^T$, $y = (y_1, y_2, \cdots, y_n)^T$ be vectors, $c$ be a real number, $m$ be an integer, define the following operations:

The Kronecker product of two matrices $W = (w_{ij})_{k\times m}$, $N = (n_{ij})_{p\times q}$ is defined by

Using the above symbols, and substituting (2.7) into Eqs (2.4)–(2.6), we can get the corresponding differentiation matrices of them:

where

$I_{m+1}$ and $I_{n+1}$ are the $m+1$ and $n+1$ identity matrices, respectively.

In general, there are three methods to apply initial-boundary conditions (1.4) and (1.5), which are the elimination method, the replacement method, and the addition method, see ^[9,10] for more details. Here, we take the replacement method to apply initial-boundary conditions.

It is remarked that Eqs (1.1)–(1.3) with other types of boundary conditions, such as the Neumann boundary condition, can also be similarly considered. It is also remarked that the most important thing in our method is to obtain the differentation matrices (2.26)–(2.28). Therefore, the method can also be generalized to study Eqs (1.1)–(1.3) with the viscous term $-u_{txx}$ and other types of terms, which will be our further studying.

3.   Error estimate of numerical solutions

In the equidistant partition, the error function between numerical solution $u_{mn}(x, t)$ and the exact solution $u(x, t)$ are defined as follows, see ^[27,30]:

where

and $u[x_i, \cdots, x_{i+d_1}, x, t]$ represents the divided difference of $u$ at $x_i, \cdots, x_{i+d_1}, x, t$, and the other two terms of (3.1) have the same meaning. It is remarked that the third term in (3.1) is a high order infinitesimal of the other two terms, see Theorem 3.1 of ^[30]. Therefore, the third term in (3.1) can be ignored when the convergence rate of the error is studied. Li ^[25] gave the following error estimate of $e(x, t)$, which was motivated by the error estimate for the univariate barycentric rational interpolation in ^[18].

Lemma 3.1. ^[25] For the error function $e(x, t)$ given by the first two terms of (3.1), assume $u(x, t)\in C^{d_1+k_1+2}[a, b]\times C^{d_2+k_2+2}[0, T]$, then

where $C$ is a positive constant and $k_1$, $k_2 = 0, 1\cdots$.

Assume $u(x_m, t_n)$ and $u(x, t)$ are numerical and exact solutions of equation $Lu(x, t) = 0$, respectively, where $L$ is a bounded operator, then we have

and

For Eq (2.1), we have the error estimate as follows.

Theorem 3.1. Assume that $u(x_m, t_n)$ satisfies $Lu(x, t) = 0$, $u(x, t)\in C^{d_1+4}[a, b]\times C^{d_2+4}[0, T]$, $L$ is a bounded operator, then

where

$d_1\geq 1$, $d_2\geq 1$, and $C$ is a positive constant.

Proof. For some given initial value $u_0$, we get

where

For $e_1(x, t)$, it follows from Lemma 3.1 that

where $\mu_i(x)$ and $\lambda_j(t)$ are defined in (3.2) and (3.3), $C_1$ is a positive constant.

For $e_2(x, t)$, by Lemma 3.1, we also get

where $C_2$ is a positive constant.

By the same way, we get that

and

where $C_3$ and $C_4$ are positive constants.

Then, combining (3.7)–(3.10), it completes the proof. □

Similarly to Theorem 3.1, we also obtain the error estimate for Eqs (2.2) and (2.3). For simplicity, we omit their proofs and only state them as follows.

Theorem 3.2. Assume that $u(x_m, t_n)$ satisfies $Lu(x, t) = 0$, $u(x, t)\in C^{d_1+6}[a, b]\times C^{d_2+4}[0, T]$, $L$ is a bounded operator, then

where

$d_1\geq 3$, $d_2\geq 1$, $C > 0$ is a constant.

Theorem 3.3. Assume that $u(x_m, t_n)$ satisfies $Lu(x, t) = 0$, $u(x, t)\in C^{d_1+8}[a, b]\times C^{d_2+4}[0, T]$, $L$ is a bounded operator, then

where

$d_1\geq 5$, $d_2\geq 1$, $C > 0$ is a constant.

Remark 3.1. From Theorems 3.1–3.3, we can see the convergence rate of errors for Eqs (2.1)–(2.3) can reach $O(h^{d_1-1}+ \tau^{d_2-1})$, $O(h^{d_1-3}+ \tau^{d_2-1})$ and $O(h^{d_1-5}+ \tau^{d_2-1})$ as $h\to 0$, $\tau\to 0$, respectively.

4.   Numerical examples

In this section, we will provide three examples using the LBRCM to solve nonlinear PDEs (1.1)–(1.3), and show the validity of the theorems. These examples are done with Matlab R2013a on a PC (Configuration: Intel(R) Core(TM) i5-6200U CPU @ 2.30GHz 2.40 GHz).

The $L_{\infty}$ norm is used for maximum absolute error, that is,

where $u^c$ and $u^e$ are the approximate and analytic solutions of equations, respectively. The convergence orders for space or time variables are defined as:

where $E_i$ represents the error during the $i$-th mesh and $h_i$, $\tau_i$ represent the maximum stride length of partition corresponding to the space or time variables in the $i$-th mesh, respectively.

Example 4.1. Consider the generalized Boussinesq equation

for the cases of $s = 1$ and $s = 2$, where $\beta_1$, $\beta_2$, $\beta_3$ are constants.

When $s = 1$, it has the initial condition

for some constant $\delta$ with $\frac{ \beta_1}{4 \beta_3{\delta}^2+1} < 0$ and the boundary condition

The exact solution is

When $s = 2$, it has the initial condition

for some constant $\delta$ with $\frac{ \beta_3}{ \beta_2} > 0$, $\frac{ \beta_1}{ \beta_3{\delta}^2+1} < 0$ and the boundary condition

The exact solution is

In this example, we always take the parameters $\beta_1 = \beta_3 = -1$, $\beta_2 = -6$, and $\delta = 0.1$ for simulation. Figures 1 and 2 give errors of the LBRCM in two kinds of meshes for $s = 1, 2$, $T = 1$, $[a, b] = [0, 1]$, interpolation parameters $d_1 = d_2 = 5$ and the number of interpolation nodes $M\times N = 20\times 20$, where $M: = m+1$ and $N: = n+1$ will be used in the rest of the paper.

In Tables 1 and 2, the maximum errors of the LBRCM in two kinds of meshes are given, where the numbers of partition nodes are $M\times N = 10\times 10$ or $20\times 20$, and the interpolation parameters are $d_1 = d_2 = 3$ or $5$, respectively. In Table 1, we fix the range of the time variable with $T = 0.5$ and let the range of the space variable changed. In Table 2, we fix the range of the space variable with $[a, b] = [0, 1]$ and let range of the time variable changed. From two tables, we see the following results: the maximum errors are small enough even though the space variable and the time variable are bigger; the maximum errors decrease with an increase in the number of nodes in both uniform and nonuniform meshes; the nonuniform mesh has a higher precision than the uniform mesh under the same condition.

In Tables 3 and 4, the maximum errors of the LBRCM in the uniform mesh are given for $T = 1$, $[a, b] = [0, 1]$. In Table 3, if the time interpolation parameter is set to $d_2 = 5$, the convergence rate of space variable can reach $O(h^{d_1+1})$ for both $s = 1$ and 2. In Table 4, if the space interpolation parameter is set to $d_1 = 5$, the convergence rate of time variable can reach $O(\tau^{d_2})$ for both $s = 1$ and 2.

In Tables 5 and 6, the maximum errors of the LBRCM in nonuniform mesh are given for $T = 1$, $[a, b] = [0, 1]$. In Table 5, if the time interpolation parameter is set to $d_2 = 5$, the convergence rate of space variable can reach $O(h^{2d_1})$ for both $s = 1$ and 2. In Table 6, if the space interpolation parameter is set to $d_1 = 5$, the convergence rate of time variable can reach $O(\tau^{2d_2})$ for both $s = 1$ and 2.

Example 4.2. Consider the generalized Boussinesq equation

for the cases of $s = 1$ and $s = 2$, where $\beta_1 < 0$, $\beta_2 < 0$, $\beta_3 < 0$.

When $s = 1$, it has the initial condition

for some constant $\delta$ and the boundary condition

The exact solution is

When $s = 2$, it has the initial condition

for some constant $\delta$ and the boundary condition

The exact solution is

In this example, we always take the parameters $\beta_1 = \beta_3 = -1$, $\beta_2 = -6$ and $\delta = 0.1$ for simulation. Figures 3 and 4 give the errors of the LBRCM in two kinds of meshes for $s = 1, 2$, $T = 1$, $[a, b] = [0, 1]$, interpolation parameters $d_1 = d_2 = 6$ and the number of interpolation nodes $M\times N = 10\times 10$. Figure 5 gives the error of the LBRCM in a uniform mesh for $s = 1, 2$, $T = 3$, $[a, b] = [0, 3]$, interpolation parameters $d_1 = d_2 = 6$, and the number of interpolation nodes $M\times N = 10\times 10$. This shows that even for large ranges of $T$ and $[a, b]$, it can also get high error precision with small interpolation parameters and small number of interpolation nodes.

In Tables 7 and 8, the maximum errors of the LBRCM in the uniform mesh are given for $T = 1$, $[a, b] = [0, 1]$. In Table 7, if we take the time interpolation parameter of $d_2 = 6$, the maximum error almost retains the same precision for small $M$, $N$, and different $d_1$ as $s = 1$ or 2. In Table 8, if we take the time interpolation parameter of $d_1 = 5$, the maximum error retains the same precision for small $M$, $N$, and different $d_2$ as $s = 1$ or 2.

In this example, we need to compute the term of $\frac{\partial^4 u}{\partial x^4}$ by $C^{(4)}$, even if there is an error from numerical calculation for the elements of $C^{(4)}$. With fewer less points $M$ = $N$ = 8, the precision can reach $10^{-07}$ for $s = 1$ or 2. By the theorem, as the convergence rate is $O(h^{d_1-3}+ \tau^{d_2-1})$, we take $d_1 > 3$ and $d_{1} = 4, 5, \cdots,$ Lack of the compute conditionfewer, we did not obtain the convergence rate.

Example 4.3. Consider the generalized Boussinesq equation

for the cases of $s = 1$ and $s = 2$, where $\beta_1$, $\beta_2$, $\beta_3$ are negative constants, $0 < \gamma\leq 1$.

When $s = 1$, it has the initial condition

and the boundary condition

The exact solution is

When $s = 2$, it has the initial condition

and the boundary condition

The exact solution is

In this example, we always take the parameters $\beta_1 = \beta_2 = \beta_3 = -0.1$, $\gamma = 0.8$, and the positive case as $s = 2$ for simulation. Figures 6 and 7 give the errors of the LBRCM in two kinds of meshes for $s = 1, 2$, $T = 1$, $[a, b] = [0, 2]$, interpolation parameters $d_1 = d_2 = 9$, and the number of interpolation nodes $M\times N = 10\times 10$. Figure 8 gives the error of the LBRCM in uniform mesh for $s = 1, 2$, $T = 4$, $[a, b] = [0, 4]$, interpolation parameters $d_1 = d_2 = 9$, and the number of interpolation nodes $M\times N = 10\times 10$. This shows that even for large ranges of $T$ and $[a, b]$, it can also get high error precision with small interpolation parameters and small number of interpolation nodes.

In Tables 9 and 10, the maximum errors of the LBRCM in uniform mesh are given for $T = 1$, $[a, b] = [-5, 5]$. In Table 9, if we take the time interpolation parameter $d_2 = 9$, the maximum error almost retains the same precision for small $M$, $N$ and different $d_1$ as $s = 1$ or 2. In Table 10, if we take the time interpolation parameter of $d_1 = 9$, the maximum error almost retains the same precision for small $M$, $N$, and different $d_2$ as $s = 1$ or 2.

In this example, we need to compute the term of $\frac{\partial^6 u}{\partial x^6}$ by $C^{(6)}$, if we take there is error from numerical calculation for the elements of $C^{(6)}$. With fewer points $M$ = $N$ = 10, the precision can reach $10^{-07}$ or $10^{-06}$ for $s = 1$ or 2, respectively. By the theorem, if the convergence rate is $O(h^{d_1-5}+ \tau^{d_2-1})$, we take $d_1 > 5$ and $d_{1} = 6, 7, \cdots,$ Due to a lack of the compute conditionfewer, we have not get the convergence rate.

5.   Conclusions

In this article, the LBRCM is applied to solve a class of generalized Boussinesq shallow-water wave equations. These nonlinear PDEs are first transformed into linear PDEs through direct linearization method. Then, the differentiation matrices of their discretization are given for computer calculation. Based on the error estimate of the barycentric interpolation, the rates of convergence for numerical solutions of those equations are also obtained for the equidistant partition. For Eq (2.1), the convergence rate can reach $O(h^{d_1-1}+ \tau^{d_2-1})$, which is confirmed by the computer simulation of Example 4.1. For Eqs (2.2) and (2.3), the convergence rats were proven to reach $O(h^{d_1-3}+ \tau^{d_2-1})$ and $O(h^{d_1-5}+ \tau^{d_2-1})$, respectively. However, in Examples 4.2 and 4.3, there are $4$-th and $6$-th partial derivatives with respect to the space variable and the errors almost retain the same precision, which makes it difficult to test the convergence rate with a low precision computer. In addition, for the nonuniform partition, the convergence rate is not proven, which will be studied later.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The research is funded by Natural Science Foundation of Shandong Province (Grant ZR2022MA003).

Conflicts of interest

The authors declare that they have no conflicts of interest.