Space–time spectral collocation method for Klein–Gordon equation

Abstract

By using the Legendre–Laguerre collocation method, we can construct a spectral collocation scheme to solve the Klein–Gordon equation on the half-line. The Laguerre function collocation method (based on the Lagrange interpolation) in space and the Legendre–Gauss–Lobatto collocation method in time are used. A Newton iterative algorithm is provided. The numerical results demonstrate the high efficiency and accuracy of suggested algorithms.

Introduction

As we know, the Klein–Gordon equation plays a significant role in many scientific applications such as solid-state physics, nonlinear optics, and quantum field theory, see Rashidinia et al. ¹ and references therein. The nonlinear Klein–Gordon (NLKG) equation, which belongs to the nonlinear wave phenomena and particularly solitary waves, has the general form^2–6:

$${\mathrm{\partial }}_t^{2}u-{\alpha }^2{\mathrm{\partial }}_x^{2}u+g(u)=0$$

(1)

where $u(x,t)$ represents the displacement of wave in the direction perpendicular to the x-axis, $\alpha$ is a known constant, and $g(u)$ is the nonlinear force. Such as, when $g(u)=\sin (u)$ , equation (1) is the well-known sine-Gordon equation. Some authors considered the solutions of the generalized equation (1) for $g(u)=au-b{u}^3$ ( $a,b$ are all the known constants), see the literature.^7–10 Moreover, the numerical solutions to the nonlinear Klein–Gordon equation with various forms of $g(u)$ have received considerable attention in existing literature, see, for example, Rashidinia et al., ¹ Bratsos, ² Mittal and Bhatia, ³ Khuri and Sayfy, ⁴ Lynch, ⁵ Li and Guo, ¹¹ Lakestani and Dehghan, ¹² and references therein. A variety of second-order finite difference schemes have been presented,^5,13–15 we note that some authors use spectral and pseudospectral methods to solve the NLKG equation.^11,16 Also, a finite-element method is introduced in Mahanti. ¹⁷ Furthermore, two numerical techniques based on the finite difference and collocation methods are presented for the solution of NLKG equation in Lakestani and Dehghan, ¹² a technique to solve NLKG equation numerically has been proposed using the cubic B-spline collocation method in Mittal and Bhatia ³ and Khuri and Sayfy. ⁴ The Chebyshev cardinal functions are used for the solution of a partial differential equation with an unknown time-dependent coefficient and one-dimensional linear hyperbolic equation in Lakestani and Dehghan ¹⁸ and Dehghan and Lakestani. ¹⁹ All of the above methods have two things in common, one is a finite interval on spatial and time direction, the other is using finite difference schemes on time direction. Usually, the error order in space and time is different for the same number of nodes on spatial and time direction, so to balance, more nodes are needed in the time direction, which is adding to the work. Thus, it is interesting to use other methods to solve the NLKG equation, such as the collocation method.

Recently, more and more attention has been drawn to the problems defined on unbounded domains. In general, one uses the spectral or pseudospectral method of Laguerre (or Hermite) for the problems defined on the half-line (or the whole line).^20–25 While the pseudospectral method is more preferable in actual computation, since it is easier to be implemented and saves work to deal with nonlinear terms. ²² Spectral collocation methods are efficient techniques for solving nonlinear equations accurately. ²⁶

In this paper, we will develop a Laguerre–Legendre spectral collection method to construct a scheme to calculate the numerical solution of NLKG equation (1) defined on the half-line with initial value conditions:

$$u(x,t){|}_{t=0}=h_0(x),{\mathrm{\partial }}_{t}u(x,t){|}_{t=0}=h_1(x)$$

(2)

and boundary condition for the half-line:

$$u(x,t){|}_{x=0}=g(t)$$

(3)

The Laguerre collocation method (based on the Lagrange interpolation) in space, and the Legendre–Gauss–Lobatto collocation method in time are used. A simple nonlinear system is suggested, which can be implemented by the Newton iteration. Numerical results demonstrate the efficiency of the proposed algorithm．

This paper is arranged as follows. In section “Lagrange interpolation polynomials based on Gauss nodes,” first, we recall the basis results related to the differential matric of Lagrange interpolation polynomial based on the Legendre–Gauss–Lobatto nodes and derive the scaled differential matrix of Lagrange interpolation polynomial based on the Laguerre–Gauss nodes. In section “Space–time spectral collocation method for Klein–Gordon equation,” we applied the Legendre–Laguerre spectral collocation method to solve NLKG equations (1) to (3). In section “Numerical result,” numerical experiments are presented to demonstrate the viability and the efficiency of the proposed method computationally. The final section is a brief conclusion for this paper.

Lagrange interpolation polynomials base on Gauss nodes

Interpolation basis functions based on Legendre–Gauss–Lobatto nodes and differential matrix

Let $I=(-1,1)$ . The Legendre polynomial of degree m is (cf. Guo ²⁷ and Jie and Tao ²⁸ )

$$L_m(\xi )={\displaystyle \frac{{(-1)}^m}{2^{m}m!}{\mathrm{\partial }}_{\xi }^m(1-{\xi }^2{)}^m,\xi \in I}$$

It is the mth eigenfunction of the singular Sturm–Liouville problem

$${\mathrm{\partial }}_{\xi }((1-{\xi }^2){\mathrm{\partial }}_{\xi }u(\xi ))+\lambda u(\xi )=0,\xi \in I$$

(4)

Related to mth eigenvalue ${\lambda }_m=m(m+1)$ .

Let ${\xi }_0=-1,{\xi }_M=1,$ and ${\xi }_j(1\le j\le M-1)$ be the roots of ${\mathrm{\partial }}_{\xi }{L}_M(\xi )=0$ . Take $w_{M+1}(\xi )=(\xi -{\xi }_0)(\xi -{\xi }_1)\cdots (\xi -{\xi }_M)$ . The Lagrange interpolation basis functions associated with the Legendre–Gauss–Lobatto nodes $\{{\xi }_j{\}}_{j=0}^M$ are as follows:

$${\varphi }_j(\xi )={\displaystyle \frac{w_{M+1}(\xi )}{(\xi -{\xi }_j){\mathrm{\partial }}_{\xi }{w}_{M+1}({\xi }_j)},j=0,1,\dots ,M}$$

(5)

Obviously

$$w_{M+1}(\xi )=c({\xi }^2-1){\mathrm{\partial }}_{\xi }{L}_M(\xi )$$

By using (4), (5) reads

$${\varphi }_j(\xi )={\displaystyle \frac{({\xi }^2-1){\mathrm{\partial }}_{\xi }{L}_M^{}(\xi )}{(\xi -{\xi }_j)L_M({\xi }_j)},j=0,1,\dots ,M}$$

(6)

Then, for $u(\xi )\in C(\overline{I})$ , the Lagrange interpolation polynomial can be written as

$$p_M(\xi )=\sum _{j=0}^{N}u({\xi }_j){\varphi }_j(\xi ),\xi \in [-1,1]$$

Differentiating the above equations m times, and setting $\xi ={\xi }_k$ gets that

$${\mathrm{\partial }}_{\xi }^m{p}_M({\xi }_k)=\sum _{j=0}^{N}u({\xi }_j){\mathrm{\partial }}_{\xi }^m{\varphi }_j({\xi }_k)$$

Let

$$D^{(m)}=({\mathrm{\partial }}_{\xi }^m{\varphi }_j^{}({\xi }_k){)}_{0\le j,k\le M},D=D^{(1)}=(d_{kj}={\mathrm{\partial }}_{\xi }^{}{\varphi }_j^{}({\xi }_k){)}_{0\le j,k\le M}$$

The entries of the first-order differentiation matrix $D=(d_{kj})$ are determined by (cf. Shen et al. ²⁹ )

$$d_{kj}=\{\begin{array}{cc}-{\displaystyle \frac{M(M+1)}{4},}& k=j=0\\ {\displaystyle \frac{L_M({\xi }_k)}{({\xi }_k-{\xi }_j)L_M({\xi }_j)},}& k\ne j\\ {\displaystyle \frac{M(M+1)}{4},}& k=j=M\\ 0,& k=j,(k,j=1,2,\dots ,N-1)\end{array}$$

(7)

Moreover, $D^{(m)}=D^m$ .

Interpolation basis functions based on Laguerre–Gauss nodes and differential matrix

Let ${\mathbb{R}}^{+}=(0,\mathrm{\infty })$ . The Laguerre polynomial of degree n is as follows^27,28:

$${\mathcal{L}}_n(y)={\displaystyle \frac{1}{n!}{\mathrm{e}}^y{\mathrm{\partial }}_y^n(y^n{\mathrm{e}}^{-y}),y\in {\mathbb{R}}^{+}}$$

It is the nth eigenfunction of the singular Sturm–Liouville problem:

$${\mathrm{\partial }}_y(y{\mathrm{e}}^{-y}{\mathrm{\partial }}_{y}u(y))+\lambda {\mathrm{e}}^{-y}u(y)=0$$

(8)

Related to the nth eigenvalue ${\lambda }_n=n.$

The Laguerre function is

$${\tilde{\mathcal{L}}}_n(y)={\mathrm{e}}^{-y/2}{\mathcal{L}}_n(y)$$

Let $y_j$ are the zeros of $y{\mathrm{\partial }}_y{\mathcal{L}}_{N+1}(y)=0$ . We take $y_j$ as the interpolation notes of the Lagrange base function, similarly, the Lagrange base function is

$${\phi }_j^{}(y)={\displaystyle \frac{{\mathrm{e}}^{-y/2}}{(N+1){\mathrm{e}}^{-y_j/2}}{\displaystyle \frac{y{\mathrm{\partial }}_y{\mathcal{L}}_{N+1}(y)}{(y-y_j){\mathrm{\partial }}_y{\mathcal{L}}_{N+1}(y_j)},j=0,1,2,\dots ,N.}}$$

Then $u(y)\in C({\overline{\mathbb{R}}}^{+})$ , the Lagrange interpolation function can be written as

$$p_N(y)=\sum _{j=0}^{N}u(y_j){\phi }_j(y),y\in \{0\}\cup {\mathbb{R}}^{+}$$

Differentiating the above equations m times, and setting $y=y_k$ gets that

$${\mathrm{\partial }}_y^m{p}_N(y_k)=\sum _{j=0}^{N}u(y_j){\mathrm{\partial }}_y^m{\phi }_j(y_k)$$

Let

$${\mathcal{D}}^{(m)}=({\mathrm{\partial }}_y^m{\phi }_j(y_k){)}_{0\le j,k\le N},\mathcal{D}={\mathcal{D}}^{(1)}=({\tilde{d}}_{kj}={\mathrm{\partial }}_y^{}{\phi }_j(y_k){)}_{0\le j,k\le N}$$

The entries of the matrix $\mathcal{D}=({\tilde{d}}_{kj}{)}_{0\le k,j\le N}$ are determined by (cf. Shen et al. ²⁹ )

$${\tilde{d}}_{k,j}^{}=\{\begin{array}{cc}{\displaystyle \frac{{\tilde{\mathcal{L}}}_{N+1}(y_k)}{(y_k-y_j){\tilde{\mathcal{L}}}_{N+1}(y_j)},}& k\ne j\\ -{\displaystyle \frac{N+1}{2},}& k=j=0\\ 0,& 1\le k=j\le N\end{array}$$

(9)

To match the asymptotic behaviors of the exact solutions as $x\to +\mathrm{\infty }$ , we introduce the scaled generalized Laguerre interpolation functions. To do this, let $\beta$ be a positive constant, we define the following Laguerre function:

$${\mathcal{L}}_n^{(\beta )}(x)={\mathcal{L}}_n(\beta x),{\tilde{\mathcal{L}}}_n^{(\beta )}(x)={\mathrm{e}}^{-\beta x/2}{\mathcal{L}}_n^{(\beta )}(x)$$

The corresponding Lagrange base function is

$${\phi }_j^{(\beta )}(x)={\displaystyle \frac{{\mathrm{e}}^{-\beta x/2}}{\beta (N+1){\mathrm{e}}^{-\beta x_j/2}}{\displaystyle \frac{x{\mathrm{\partial }}_x{\mathcal{L}}_{N+1}^{(\beta )}(x)}{(x-x_j){\mathrm{\partial }}_x{\mathcal{L}}_{N+1}^{(\beta )}(x_j)},j=0,1,2,\dots ,N}}$$

Let

$${\mathcal{D}}_{\beta }^{(m)}=({\mathrm{\partial }}_x^m{\phi }_j^{(\beta )}(x_k){)}_{0\le j,k\le N},{\mathcal{D}}_{\beta }={\mathcal{D}}_{\beta }^{(1)}=({\tilde{d}}_{kj}^{(\beta )}={\mathrm{\partial }}_x^{}{\phi }_j^{(\beta )}(x_k){)}_{0\le j,k\le N}$$

The entries of the matrix ${\mathcal{D}}_{\beta }=({\tilde{d}}_{kj}^{(\beta )}{)}_{0\le k,j\le N}$ are determined by

$${\tilde{d}}_{k,j}^{(\beta )}=\{\begin{array}{cc}{\displaystyle \frac{{\tilde{\mathcal{L}}}_{N+1}^{(\beta )}(x_k)}{(x_k-x_j){\tilde{\mathcal{L}}}_{N+1}^{(\beta )}(x_j)},}& k\ne j\\ -{\displaystyle \frac{\beta (N+1)}{2},}& k=j=0\\ 0,& 1\le k=j\le N\end{array}$$

(10)

Also, ${\mathcal{D}}_{\beta }^{(m)}=-{\mathcal{D}}_{\beta }^m.$

Space–time spectral collocation method for the Klein–Gordon equation

We consider initial-boundary value problems of the NLKG equation as following ² :

$$\{\begin{array}{ccc}{\mathrm{\partial }}_t^{2}v+pv-q{v}^3-c_0^2{\mathrm{\partial }}_x^{2}v=0,& x\in R^{+},& 0<t\le T\\ v(x,t)\to 0,& x\to \mathrm{\infty },& 0\le t\le T\\ v(x,0)=h_0(x),& x\ge 0& \\ {\mathrm{\partial }}_{t}v(x,0)=h_1(x),& x\ge 0& \\ v(0,t)=g(t),& & 0\le t\le T\end{array}$$

(11)

where $p,q$ , and $c_0$ are some given constants. We need the variable transformation for time t. Let $t=(\xi +1)T/2,\xi \in [-1,1]$ , then we can obtain a new equation:

$$\{\begin{array}{ccc}{\displaystyle \frac{4}{T}{\mathrm{\partial }}_{\xi }^{2}u+pu-q{u}^3-c_0^2{\mathrm{\partial }}_{\xi }^{2}u=0,}& x\in R^{+},& -1<\xi \le 1\\ u(x,\xi )\to 0,& x\to \mathrm{\infty },& -1\le \xi \le 1\\ u(x,-1)=h_0(x),& x\ge 0& \\ 2{\mathrm{\partial }}_{t}u(x,-1)=h_1(x),& x\ge 0& \\ u(0,\xi )=g({\displaystyle \frac{(\xi +1)T}{2}),}& & -1\le \xi \le 1\end{array}$$

(12)

we expand the numerical solution as

$$u_{M,N}(x,\xi )=\sum _{i=0}^M\sum _{j=0}^N{\hat{u}}_{i,j}{\phi }_j^{(\beta )}(x){\varphi }_i(\xi ),u_{i,j}=u(x_j,{\xi }_i)$$

to approximate the exact solution of (12). Take $x=x_l,\xi ={\xi }_k$ in (12), then (12) can be transformed into nonlinear systems

$$\{\begin{array}{c}{\displaystyle \frac{4}{T^2}\sum _{i=0}^M\sum _{j=0}^N{\hat{u}}_{i,j}{\phi }_j^{(\beta )}(x_l){\mathrm{\partial }}_{\xi }^2{\varphi }_i({\xi }_k)+p\sum _{i=0}^M\sum _{j=0}^N{\hat{u}}_{i,j}{\phi }_j^{(\beta )}(x_l){\varphi }_i({\xi }_k)}\\ -q(\sum _{i=0}^M\sum _{j=0}^N{\hat{u}}_{i,j}{\phi }_j^{(\beta )}(x_l){\varphi }_i({\xi }_k){)}^3-c_0^2\sum _{i=0}^M\sum _{j=0}^N{\hat{u}}_{i,j}{\mathrm{\partial }}_x^2{\phi }_j^{(\beta )}(x_l){\varphi }_i({\xi }_k)\\ =0,\begin{array}{cc}& \end{array}k=2,3,\dots ,M;l=1,2,\dots ,N,\\ \sum _{i=0}^M\sum _{j=0}^N{\hat{u}}_{i,j}{\phi }_j^{(\beta )}(x_l){\varphi }_i({\xi }_0)=h_0(x_l),l=0,1,2,\dots ,N,\\ {\displaystyle \frac{2}{T}\sum _{i=0}^M\sum _{j=0}^N{\hat{u}}_{i,j}{\phi }_j^{(\beta )}(x_l){\mathrm{\partial }}_{\xi }^{}{\varphi }_i({\xi }_0)=h_1(x_l),l=0,1,2,\dots ,N,}\\ \sum _{i=0}^M\sum _{j=0}^N{\hat{u}}_{i,j}{\phi }_j^{(\beta )}(x_0){\varphi }_i({\xi }_k)=g({\displaystyle \frac{({\xi }_k+1)T}{2}),k=0,1,2,\dots ,M.}\end{array}$$

(13)

For convenience, let

$${\hat{d}}_{k,i}={\mathrm{\partial }}_{\xi }^2{\varphi }_i({\xi }_k),0\le k,i\le M,{\hat{d}}_{l,j}^{(\beta )}={\mathrm{\partial }}_x^2{\phi }_j(x_l),0\le l,j\le N$$

Then (13) reads

$$\{\begin{array}{c}{\displaystyle \frac{4}{T}\sum _{i=0}^M{\hat{u}}_{i,l}{\hat{d}}_{k,i}+p{\hat{u}}_{k,l}-q({\hat{u}}_{k,l}{)}^3-c_0^2\sum _{j=0}^N{\hat{u}}_{k,j}{\hat{d}}_{l,j}^{(\beta )}=0}\\ k=2,3,\dots ,M;l=1,2,\dots ,N\\ {\hat{u}}_{0,l}=h_0(x_l),{\displaystyle \frac{2}{T}\sum _{i=0}^M{\hat{u}}_{i,l}{d}_{0,i}=h_1(x_l),l=,0,1,\dots ,N}\\ {\hat{u}}_{k,0}=g\left({\displaystyle \frac{({\xi }_k+1)T}{2}}\right),k=0,1,\dots ,M\end{array}$$

(14)

By virtue of the initial value conditions

$${\hat{u}}_{1,l}=\left({\displaystyle \frac{T{h}_1(x_l)}{2}-h_0(x_l)d_{0,0}-{\hat{u}}_{2,l}{d}_{0,2}-{\hat{u}}_{3,l}{d}_{0,3}-\cdots -{\hat{u}}_{M,l}{d}_{0,M}}\right)/d_{0,1},l=0,1,\dots ,M,$$

Taking the above formula into the first equation of (14) and using the boundary condition, we can have that

$$\begin{array}{rl}& {\displaystyle \frac{4}{T^2}({\hat{u}}_{2,l}{\hat{d}}_{k,2}+{\hat{u}}_{3,l}{\hat{d}}_{k,3}+\dots +{\hat{u}}_{M,l}{\hat{d}}_{k,M})}\\ & -{\displaystyle \frac{4}{T^2{d}_{0,1}}(\hat{u}{}_{2,l}{d}_{0,2}+\hat{u}{}_{3,l}{d}_{0,3}+\cdots +\hat{u}{}_{M,l}{d}_{0,M}){\hat{d}}_{k,1}}\\ & -c_0^2({\hat{u}}_{k,1}{\hat{d}}_{l,1}^{(\beta )}+{\hat{u}}_{k,2}{\hat{d}}_{l,2}^{(\beta )}+\cdots +{\hat{u}}_{k,N}{\hat{d}}_{l,N}^{(\beta )})+p\hat{u}{}_{k,l}-q{\hat{u}}_{k,l}^3\\ & =c_0^{2}g({\displaystyle \frac{({\xi }_k+1)T}{2}){\hat{d}}_{l,0}^{(\beta )}+{\displaystyle \frac{4d_{0,0}}{T^2{d}_{0,1}}{\hat{d}}_{k,1}{h}_0(x_l)}}\\ & -{\displaystyle \frac{4}{T^2}{h}_0(x_l){\hat{d}}_{k,0}-{\displaystyle \frac{2}{T{d}_{0,1}}{\hat{d}}_{k,1}{h}_1(x_l)}}\\ & l=1,2,\dots ,N,k=2,3,\dots ,M\end{array}$$

Let

$$X=\left(\begin{array}{cccc}{\hat{u}}_{2,1}& {\hat{u}}_{2,2}& \cdots & {\hat{u}}_{2,N}\\ {\hat{u}}_{3,1}& {\hat{u}}_{3,2}& \cdots & {\hat{u}}_{3,N}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{u}}_{M,1}& {\hat{u}}_{M,2}& \cdots & {\hat{u}}_{M,N}\end{array}\right)$$

We can express the above equations by a matrix form:

$$AX-c_0^{2}X{B}^{\mathrm{T}}=-pX+qX.*X.*X-F-H+Q+c_0^{2}C$$

(15)

where

$$A={\displaystyle \frac{4}{T^2}\left(\left(\begin{array}{cccc}{\hat{d}}_{2,2}& {\hat{d}}_{2,3}& \cdots & {\hat{d}}_{2,M}\\ {\hat{d}}_{3,2}& {\hat{d}}_{3,3}& \cdots & {\hat{d}}_{3,M}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{d}}_{M,2}& {\hat{d}}_{M,3}& \cdots & {\hat{d}}_{M,M}\end{array}\right)-{\displaystyle \frac{1}{d_{0,1}}\left(\begin{array}{c}{\hat{d}}_{2,1}\\ {\hat{d}}_{3,1}\\ ⋮\\ {\hat{d}}_{M,1}\end{array}\right)\left(\begin{array}{cccc}{d}_{0,2}& d_{0,3}& \cdots & d_{0,M}\end{array}\right)}\right)}$$

$$\begin{array}{r}B=\left(\begin{array}{cccc}{\hat{d}}_{1,1}^{(\beta )}& {\hat{d}}_{1,2}^{(\beta )}& \cdots & {\hat{d}}_{1,N}^{(\beta )}\\ {\hat{d}}_{2,1}^{(\beta )}& {\hat{d}}_{2,2}^{(\beta )}& \cdots & {\hat{d}}_{2,N}^{(\beta )}\\ ⋮& ⋮& \ddots & ⋮\\ {\hat{d}}_{N,1}^{(\beta )}& {\hat{d}}_{N,2}^{(\beta )}& \cdots & {\hat{d}}_{N,N}^{(\beta )}\end{array}\right),C=\left(\begin{array}{c}g({\displaystyle \frac{({\xi }_2+1)T}{2})}\\ g({\displaystyle \frac{({\xi }_3+1)T}{2})}\\ ⋮\\ g({\displaystyle \frac{({\xi }_M+1)T}{2})}\end{array}\right)(\begin{array}{cc}{\hat{d}}_{1,0}^{(\beta )}& {\hat{d}}_{1,2}^{(\beta )}\end{array}\cdots {\hat{d}}_{1,N}^{(\beta )})\end{array}$$

$$F={\displaystyle \frac{4}{T^2}\left(\begin{array}{c}{\hat{d}}_{2,0}\\ {\hat{d}}_{3,0}\\ ⋮\\ {\hat{d}}_{M,0}\end{array}\right)(h_0(x_1)h_0(x_2)\cdots h_0(x_N)),H={\displaystyle \frac{2}{T{d}_{0,1}}\left(\begin{array}{c}{\hat{d}}_{2,1}\\ {\hat{d}}_{3,1}\\ ⋮\\ {\hat{d}}_{M,1}\end{array}\right)(\begin{array}{cc}{h}_1(x_1)& h_1(x_2)\end{array}\cdots h_1(x_N))}}$$

$$Q={\displaystyle \frac{4d_{0,0}}{T^2{d}_{0,1}}\left(\begin{array}{c}{\hat{d}}_{2,1}\\ {\hat{d}}_{3,1}\\ ⋮\\ {\hat{d}}_{M,1}\end{array}\right)(h_0(x_1)h_0(x_2)\cdots h_0(x_N))}$$

“.*” means the times of corresponding entries of matrixes A and B.

Moreover, let $Y=\mathrm{vec}(X)={\left(\begin{array}{c}{\hat{u}}_{2,1}{\hat{u}}_{2,2}\cdots {\hat{u}}_{2,N}{\hat{u}}_{3,1}{\hat{u}}_{3,2}\cdots {\hat{u}}_{3,N}\cdots {\hat{u}}_{M,1}{\hat{u}}_{M,2}\cdots {\hat{u}}_{M,N}\end{array}\right)}^{\mathrm{T}}$ be $(M-1)$ N-dimensional vector, then (15) is equivalent of the following formulation:

$$((A+p{E}_{M-1})\otimes E_N-c_0^2{E}_{M-1}\otimes B)Y=qY.*Y.*Y+\mathrm{vec}(Q+c_0^{2}C-F-H)$$

(16)

where “E_n” stands for n-order identity matrix, “ $\otimes$ ” means the times of Kronecker and $\mathrm{vec}(Q+c_0^{2}C-F-H)$ is similar to $\mathrm{vec}(X).$

Numerical result

Obviously, formula (15) is a nonlinear system, we use the Newton iteration to solve it numerically. Let

$$S(Y)=\left(\begin{array}{ccc}{y}_1^2& & 0\\ & \ddots & \\ 0& & y_{(M-1)N}^2\end{array}\right)$$

$$R(Y)=((A+p{E}_{M-1})\otimes E_N-c_0^2{E}_{M-1}\otimes B)Y-qY.*Y.*Y-\mathrm{vec}(Q+c_0^{2}C-F-H)$$

$$D(Y)=(A+p{E}_{M-1})\otimes E_N-c_0^2{E}_{M-1}\otimes B-3qS(Y)$$

the Newton iteration scheme is as follows:

$$Y_{n+1}=Y_n-(D(Y_n){)}^{-1}R(Y_n),n=0,1,2,\dots$$

In actual computation, we take

$$Y_0={\left(\begin{array}{c}{h}_0(x_1)h_0(x_2)\cdots h_0(x_N)h_0(x_1)h_0(x_2)\cdots h_0(x_N)\cdots \\ h_0(x_1)h_0(x_2)\cdots h_0(x_N)\end{array}\right)}^{\mathrm{T}}$$

and the end condition of iteration: $\mathrm{\forall }\epsilon >0,\Vert X_{n+1}-X_n{\Vert }_{\mathrm{\infty }}<\epsilon .$

The exact solution of (11) is as follows ² :

$$v(x,t)=\sqrt{{\displaystyle \frac{2p}{q}}}sech\sqrt{{\displaystyle \frac{p}{c_0^2-c^2}}}(x-ct)$$

with initial value conditions

$$\begin{array}{rl}& h_0(x)=\sqrt{{\displaystyle \frac{2p}{q}}}sech\sqrt{{\displaystyle \frac{p}{c_0^2-c^2}}}x\\ & h_1(x)=c\sqrt{{\displaystyle \frac{2p}{q}}}sech\left(\sqrt{{\displaystyle \frac{p}{c_0^2-c^2}}}x\right)\\ & \tanh \left(\sqrt{{\displaystyle \frac{p}{c_0^2-c^2}}}x\right)\end{array}$$

and boundary conditions

$$g(t)=\sqrt{{\displaystyle \frac{2p}{q}}}sech\sqrt{{\displaystyle \frac{p}{c_0^2-c^2}}}(-ct)$$

We transform it into

$$u(x,\xi )=\sqrt{{\displaystyle \frac{2p}{q}}}sech\sqrt{{\displaystyle \frac{p}{c_0^2-c^2}}}\left(x-c{\displaystyle \frac{(\xi +1)T}{2}}\right)$$

(17)

For the description of the numerical accuracy, we introduce the maximum absolute error

$$E_{M,N}=\underset{2\le k\le M,1\le l\le N}{max}|u_{M,N}(x_l,{\xi }_k)-u(x_l,{\xi }_k)|$$

Let $\beta =5$ in (10) and $T=1$ in (11). In Figures 1 and 2, we plot the errors ${\log }_{10}{E}_{M,N}$ with $p=1,c_0=1,q=40$ , $c=1/50$ , and $c=1/100$ in (17), respectively, versus ${\log }_{10}N$ . Obviously, we can obtain the high-accuracy numerical solution even for smaller nodes $k\le 6\times 15=90$ , which indicates that the Klein–Gordon equation on a half-line for the x-direction can be used to effectively simulate.

Figure 1.

Error ${\log }_{10}{E}_{M,N}$ with $\beta =5,p=1,q=40,c=1/50,c_0=1,M=6$ .

Figure 2.

Error ${\log }_{10}{E}_{M,N}$ with $\beta =5,p=1,q=40,c=1/100,c_0=1,M=6$ .

Let $\beta =5.175$ in (10) and $T=5$ in (11). In Figure 3, we plot the errors ${\log }_{10}{E}_{M,N}$ with $p=1,c_0=1,q=40$ , $c=1/100$ in (17) versus N, which show the scheme (16) can also give better numerical results, even for larger time $T.$

Figure 3.

Error ${\log }_{10}{E}_{M,N}$ with $\beta =5.175,p=1,q=40,c=1/100,c_0=1,M=6$ .

Let $T=1$ in (11). In Figure 4, we plot the errors ${\log }_{10}{E}_{M,N}$ with $p=1,c_0=1,q=40$ , $c=1/50$ in (17) versus N and differential $\beta$ . It seems that the larger the $\beta$ , the smaller the error. Until today, we did not know how to choose an appropriate value of $\beta$ .

Figure 4.

Error ${\log }_{10}{E}_{M,N}$ with $p=1,q=40,c=1/50,c_0=1,M=6$ .

Let $T=1$ in (11). In Figures 5 and 6, we plot the exact solution $v(x,t)$ with $p=1,c_0=1,q=40$ , and $c=1/50$ in (17) versus the numerical solution with $M=10,N=15$ , and $\beta =5$ . It shows the numerical solution fits the exact solution very well.

Figure 5.

The exact solution with $p=1,q=40,c=1/50,c_0=1$ of (11).

Figure 6.

The numerical solution with $p=1,q=40,c=1/50,c_0=1,M=10,N=15$ of (11).

Conclusions

In this paper, we use the Legendre–Laguerre spectral collection method to construct a spectral collection scheme to solve the NLKG equation. We consider the NLKG equation defined on the half-line in the spatial direction, which is different from many existing literatures that only consider the NLKG equation defined on a finite interval. Furthermore, we need only a few space and time nodes to get high-accuracy numerical results. Although we only consider the NLKG equation in this paper, the suggested method will also be applicable to many other nonlinear problems.