Numerical meshless solution of high-dimensional sine-Gordon equations via Fourier HDMR-HC approximation

TL;DR

This paper introduces a meshless numerical method combining HDMR and Fourier HC approximation to efficiently solve high-dimensional sine-Gordon equations with stability and accuracy.

Contribution

It proposes a novel implicit meshless scheme using Fourier HDMR-HC to handle high-dimensional SGEs with sparse matrices and large time steps.

Findings

Stable large time-step solutions achieved

Reduced number of nodes needed for accuracy

Effective simulation of high-dimensional SGEs

Abstract

In this paper, an implicit time stepping meshless scheme is proposed to find the numerical solution of high-dimensional sine-Gordon equations (SGEs) by combining the high dimensional model representation (HDMR) and the Fourier hyperbolic cross (HC) approximation. To ensure the sparseness of the relevant coefficient matrices of the implicit time stepping scheme, the whole domain is first divided into a set of subdomains, and the relevant derivatives in high-dimension can be separately approximated by the Fourier HDMR-HC approximation in each subdomain. The proposed method allows for stable large time-steps and a relatively small number of nodes with satisfactory accuracy. The numerical examples show that the proposed method is very attractive for simulating the high-dimensional SGEs.

Tables

Table 1. Table 1: Errors and condition number κ 𝜅 \kappa for ( 2 + 1 ) 2 1 (2+1) -D SGE

	$L_{\infty }$-error		RMS-error
Time($s$)	RBF	HDMR-HC	RBF	HDMR-HC	$\kappa$
1.0	0.0670	0.0326	0.0050	0.0043	$3.4\times {10}^5$
3.0	0.0834	0.0343	0.0103	0.0045	$5.0\times {10}^5$
5.0	0.1015	0.0355	0.0145	0.0045	$4.2\times {10}^5$
7.0	0.1516	0.0368	0.0187	0.0047	$6.3\times {10}^5$

Table 2. Table 2: Errors and condition number κ 𝜅 \kappa for ( 5 + 1 ) 5 1 (5+1) -D SGE

	$L_{\infty }$-error		RMS-error
Time($s$)	RBF	HDMR-HC	RBF	HDMR-HC	$\kappa$
1.0	0.2071	0.1083	0.0130	0.0076	$4.4\times {10}^6$
3.0	0.1956	0.0910	0.0143	0.0079	$3.7\times {10}^6$
5.0	0.2132	0.1205	0.0165	0.0080	$4.5\times {10}^6$
7.0	0.2203	0.1124	0.0187	0.0081	$4.9\times {10}^6$

Equations

$$u_{tt}(\mathbf{x},t)+\beta u_t(\mathbf{x},t)=\Delta u(\mathbf{x},t)-\psi (\mathbf{x})\sin (u(\mathbf{x},t)),$$

$$\mathbf{x}=(x_1,x_2,\cdots ,x_n)\in \Omega ,t>0,$$

$$u(\mathbf{x},0)=v_1(\mathbf{x}),$$

$$u_t(\mathbf{x},0)=v_2(\mathbf{x}),$$

$$\frac{\partial u}{\partial \mathbf{l}}(\mathbf{x},t)=w(\mathbf{x},t),\mathbf{x}\in \Gamma ,t>0$$

$$u(\mathbf{x},t)=4\arctan \left[C\exp \left(\sum _{i=1}^n{a}_i{x}_i-bt\right)\right]$$

$${\mathbb{T}}_{\lambda }^n:={\left[-\frac{1}{2\lambda },\frac{1}{2\lambda }\right)}^n$$

$${\mathcal{W}}_1^s({\mathbb{T}}_{\lambda }^n)=\{f\in \mathcal{L}({\mathbb{T}}_{\lambda }^n):{\partial }^{r}f\in \mathcal{L}({\mathbb{T}}_{\lambda }^n),|r{|}_{\infty }⩽s\}$$

$$\Vert f{\Vert }_{{\mathcal{W}}_1^s({\mathbb{T}}_{\lambda }^n)}=\sum _{0⩽|r{|}_{\infty }⩽s}\Vert {\partial }^{r}f{\Vert }_{\mathcal{L}({\mathbb{T}}_{\lambda }^n)}$$

$$|r{|}_{\infty }:=\underset{1⩽j⩽n}{\max }{r}_j$$

$${\partial }^{r}f:=\frac{{\partial }^{|r{|}_1}f}{\partial x_1^{r_1}\cdots \partial x_n^{r_n}}$$

$$f(\mathbf{x})\sim \sum _{m\in {\mathbb{Z}}^n}{c}_m{\mathrm{e}}^{2\pi \mathrm{i}\lambda m\cdot \mathbf{x}},\mathbf{x}\in {\mathbb{T}}_{\lambda }^n$$

$$c_m=c_m(f)={\lambda }^n{\int }_{{\mathbb{T}}_{\lambda }^n}f(\mathbf{x}){\mathrm{e}}^{-2\pi \mathrm{i}\lambda m\cdot \mathbf{x}}\mathrm{d}\mathbf{x},m\in {\mathbb{Z}}^n.$$

$${\Lambda }_{j_1,\cdots ,j_v}:=$$

$$\text{and}{m}_l=0\forall l\notin \{j_1,\cdots ,j_v\},l\in \{1,\cdots ,n\}\}$$

$${\Lambda }_{j_1,\cdots ,j_v}^{k_{\lambda }}:=\{m\in {\Lambda }_{j_1,\cdots ,j_v}:k_{\lambda }⩽(2\pi {)}^v{\lambda }_{j_1}\cdots {\lambda }_{j_v}|m_{j_1}\cdots m_{j_v}|<k_{\lambda }+1,k_{\lambda }\in \mathbb{N}\}.$$

$${\mathbb{Z}}^n={\Lambda }_0+\sum _j{\Lambda }_j+\sum _{j_1<j_2}{\Lambda }_{j_1,j_2}+\cdots +{\Lambda }_{j_1,\cdots ,j_n}$$

$$\sum _{m\in {\mathbb{Z}}^n}{c}_m{\mathrm{e}}^{2\pi i\lambda m\cdot x}=$$

$$+\sum _{j_1<j_2}\sum _{m\in {\Lambda }_{j_1,j_2}^{k_{\lambda }}}{c}_m{\mathrm{e}}^{2\pi i\lambda m\cdot x}+\cdots +\sum _{m\in {\Lambda }_{j_1,\cdots ,j_n}^{k_{\lambda }}}{c}_m{\mathrm{e}}^{2\pi i\lambda m\cdot x})$$

$${\Lambda }^{K_{\lambda }}={\Lambda }_0+\sum _{k_{\lambda }=1}^{K_{\lambda }}\left(\sum _j{\Lambda }_j^{k_{\lambda }}+\sum _{j_1<j_2}{\Lambda }_{j_1,j_2}^{k_{\lambda }}+\cdots +{\Lambda }_{j_1,\cdots ,j_d}^{k_{\lambda }}\right),$$

$${\Lambda }^{K_{\lambda }}:=\left\{m\in {\mathbb{Z}}^n:\prod _{j=1}^n\max \{2\pi \lambda |m_j|,1\}⩽K_{\lambda }\right\},$$

$$S_{K_{\lambda }}(\mathbf{x},f)=\sum _{m\in {\Lambda }^{K_{\lambda }}}{c}_m{\mathrm{e}}^{2\pi \mathrm{i}\lambda m\cdot \mathbf{x}},\mathbf{x}\in {\mathbb{T}}_{\lambda }^n.$$

$$F_{\lambda ,m}(\mathbf{x})={\mathrm{e}}^{2\pi \mathrm{i}\lambda m\cdot \mathbf{x}}$$

$$S_{K_{\lambda }}(\mathbf{x},f)=\sum _{m=1}^M{c}_m{F}_{\lambda ,m}(\mathbf{x}),\mathbf{x}\in {\mathbb{T}}_{\lambda }^n$$

$$|{\partial }^r{S}_{K_{\lambda }}(\mathbf{x},f)-{\partial }^{r}f(\mathbf{x})|⩽C_{n,\lambda ,s,p,\gamma ,\delta }\Vert f{\Vert }_{{\mathcal{W}}_1^s({\mathbb{T}}_{\lambda }^n)}{K}_{\lambda }^{-(s-\gamma -p-1)},|r{|}_{\infty }=p$$

$$u_j(\mathbf{x})=J_{\lambda }(\mathbf{x}-{\mathbf{x}}_j)u(\mathbf{x})=J_{j}u(\mathbf{x}),\mathbf{x}\in {\Omega }_j$$

$${\Omega }_j\subset \text{supp}{J}_{\lambda }(\mathbf{x}-{\mathbf{x}}_j)\text{and}\text{supp}{J}_{\lambda }(\mathbf{x})\subset {\mathbb{T}}_{\lambda }^n$$

$$\sum _{j=1}^D{J}_j=1,\forall \mathbf{x}\in \Omega \subseteq {\mathbb{R}}^n$$

$${\hat{u}}_j(\mathbf{x})=\sum _{m_j=1}^{M_j}{c}_{m_j}{F}_{\lambda ,m_j}(\mathbf{x}-{\mathbf{x}}_j)={{\mathbf{F}}_{\mathbf{\lambda }}}_j^{\mathrm{T}}(\mathbf{x}){\mathbf{c}}_j$$

$${\mathbf{U}}_j={\mathbf{F}}_j{\mathbf{c}}_j$$

$${\mathbf{U}}_j=[u_1,u_2,\cdots ,u_{M_j}{]}^{\mathrm{T}}$$

$${\mathbf{c}}_j=[c_1,c_2,\cdots ,c_{M_j}{]}^{\mathrm{T}}$$

Full text

††journal: J Math Chem

1 Introduction

The sine-Gordon equation (SGE) is a nonlinear hyperbolic partial differential equation (PDE) involving the d’Alembert operator and the sine of the unknown function, and the SGE plays an important role in many mathematical physics applications. It was originally introduced by Bour (1862) and rediscovered by Frenkel and Kontorova (1939). Further details about the background and applications of the SGE can be found in Aero et al. (2009); Drazin and Johnson (1989); Scott et al. (1973); Taleei and Dehghan (2014). In this paper, we focus on developing an effective means for the numerical solution of the SGE in an arbitrary number of dimensions. A $(n+1)$-dimensional SGE generally takes the form:

[TABLE]

where $n$ is a positive integer, $\Delta$ is the Laplacian operator in $n$ spatial dimensions, $\Omega=[a_{1},b_{1}]\times[a_{2},b_{2}]\times\cdots\times[a_{n},b_{n}]\subseteq\mathbb{R}^{n}$. The initial conditions associated with Eq.(1) are given by

[TABLE]

and the Neumann boundary conditions are

[TABLE]

where $\boldsymbol{l}$ denotes the (typically exterior) normal to the boundary of the domain, and $\Gamma$ is the boundary of $\Omega$, i.e $\Gamma=\partial\Omega$. The real parameter $\beta\geqslant 0$ weights the dissipative term. When $\beta=0$, Eq.(1) reduces to an undamped SGE in $n$ spatial variables, while when $\beta>0$, the damped SGE is obtained. The function $\psi(\boldsymbol{x})$ can be interpreted as a Josephson current density, while $v_{1}$ and $v_{2}$ in Eqs.(2) and (3) represent wave modes or the kink and velocity, respectively.

The $(1+1)$-dimensional SGE first appeared in a strictly mathematical context in differential geometry regarding the theory of surfaces of constant curvature Lamb (1971). Moreover, it is well known that the $(1+1)$-dimensional SGE arises in many important systems, such as the Thirring model, the Coulomb gas system and the ferromagnetic $XY$ model, etc. Minnhagen (1985, 1987); Ni et al. (1990). Because of its wide applications, the $(1+1)$-dimensional SGE has been studied with a variety of numerical methods, including finite difference methods (FDM) and finite element methods (FEM), etc. Argyris and Haase (1987). Recently, additional solution methods have been proposed including collocation Dehghan and Shokri (2008a); Lakestani and Dehghan (2010), the boundary integral approach Dehghan and Mirzaei (2008a); Dehghan and Ghesmati (2010a), and a combination of the finite difference with the the diagonally implicit Runge-Kutta-Nyström (DIRKN) method Mohebbi and Dehghan (2010), etc.

There is recent interest in the SGE in higher dimensions. As Barone et al. (1971) pointed out, Eq.(1) also has been applied in many branches of physics for the $n=2$ and $3$ cases. The exact solutions for the undamped SGE in higher dimensions have been obtained by Hirota’s method Hirota (1973), Lamb’s method Zagrodzinsky (1979), the Bäcklund transformation Christiansen and Olsen (1979) and Painlevé transcendents Kaliappan and Lakshmanan (1979), etc. Moreover, numerical solutions for the $(2+1)$-dimensional undamped SGE have been proposed by Christiansen and Lomdahl (1981) using a generalized leapfrog method, Guo et al. (1986) using two finite difference schemes, Argyris et al. (1991) using finite elements. Xin (2000) studied the SGE as an asymptotic reduction of the two level dissipationless Maxwell-Bloch system, Sheng et al. (2005) presented a numerical method with a split cosine scheme, and Bratsos (2005) used a three-time level fourth-order explicit finite difference scheme to solve the undamped SGE. Following a similar approach, Bratsos (2007) transformed the SGE to a second-order initial value problem with the help of the method of lines. Numerical approaches for the damped SGE were proposed by Nakajima et al. (1974) who considered dimensionless loss factors and unitless normalized bias, and Gorria et al. (2004) investigated nonlinear wave propagation in a planar wave guide consisting of two rectangular regions joined by a bent domain of constant curvature using as a model of the kink solution to the SGE. Additionally, Dehghan and Mirzaei (2008b) developed the dual reciprocity boundary element method for both the undamped and damped $(2+1)$-dimensional SGE, and Jiwari et al. (2012b) obtained a numerical scheme based on a polynomial differential quadrature method.

Although the SGE is nonintegrable except for $n=1$, some properties and exact solutions for the $(n+1)$-dimensional SGE have been obtained by various methods. Kobayashi and Izutsu (1976) extensively studied the exact traveling wave solutions of the SGE in the field of theoretical physics. Many additional mathematical methods have been proposed for finding traveling wave solutions of the SGE. Feng (2004) applied the Painlevé analysis to the study of an approximate SGE and its traveling solitary wave solution in $(n+1)$-dimensional space. With the help of exact solutions to the cubic nonlinear Klein-Gordon fields, Lou et al. (2005) studied the exact solutions for the $(n+1)$-dimensional SGE. Adopting $\beta=0$ and $\psi(\boldsymbol{x})=1$, de la Hoz and Vadillo (2012) generalized the exact soliton solution for the $(n+1)$-dimensional SGE:

[TABLE]

where $\sum_{i=1}^{n}a_{i}^{2}=1+b^{2}$. By adopting the proper ansatz, more general solutions can be obtained for the multi-dimensional SGE, including the three-dimensional case allowing for non-constant $C$ Aero et al. (2009). Obtaining the exact solution for the general SGE would be ideal, but unfortunately it is very difficult for practical engineering problems that are usually complex in nature. Despite numerical methods commonly used in many types of linear and nonlinear PDEs, de la Hoz and Vadillo (2012) remarked that there were no references to the numerical treatment of the SGE for dimensions larger than three, which motivated them to propose a numerical method for the $n$-dimensional SGE based on using operational matrices.

Many standard numerical methods for solving PDEs are widely used in engineering, but they usually require the construction and update of a mesh, which is an inherent disadvantage. In order to overcome these difficulties, recently the meshless numerical method has attracted attention. This method can establish a system of algebraic equations over the entire problem domain without using a predefined mesh. Rather, a set of scattered nodes, called field nodes, are used within the problem domain as well as on the boundaries of the domain Liu and Gu (2005). The meshless method does not require a priori information about the relationship between the nodes for the interpolation or approximation of unknown functions over the field of variables Liu and Gu (2005). Employing the meshless method, Dehghan and Shokri (2008a) studied the one-dimensional nonlinear SGE and used Thin Plate Spline Radial Basis Functions (TPS-RBF) to approximate the solution, and they also applied the TPS-RBF method to both the Klein-Gordon equation Dehghan and Shokri (2009) and the two-dimensional SGE Dehghan and Shokri (2008b). A series of meshless approaches have been presented Asgari and Hosseini (2013); Dehghan and Ghesmati (2010b); Jiang and Wang (2012); Jiwari et al. (2012a); Karamanli and Mugan (2013); Mirzaei and Dehghan (2009, 2010); Pekmen and Tezer-Sezgin (2012); Shao and Wu (2014); Taleei and Dehghan (2014). Moreover, since the nodal distribution for most existing meshless methods is preassigned, Xu et al. (2015) proposed a numerical two-step meshless method for soliton-like structures based on the optimal sampling density of kernel interpolation.

In dealing with high-dimensional PDEs, obtaining good quality approximate solutions is a difficult problem because of the so-called ‘curse of dimensionality’. High dimensional model representation (HDMR) Li et al. (2001); Rabitz and Alis (1999, 2000); Rabitz et al. (1998); Sobol (1993) provides a viable approach based on the fact that high-dimensional functions often can be efficiently expressed as sums of low-dimensional functions. The HDMR decomposition is also well known in statistics as the ANOVA (analysis of variance) decomposition Efron and Stein (1981); Fisher (1925); Griebel and Holtz (2010); Stone (1994). In recent years, the HDMR decomposition has been under rapid development becoming an important tool for understanding high-dimensional functions Griebel and Holtz (2010); Griebel et al. (2010, 2013); Li et al. (2001); Luo et al. (2014a, b); Rabitz and Alis (1999, 2000); Rabitz et al. (1998); Sobol (1993). In this paper, we will use a HDMR decomposition in conjunction with the Fourier hyperbolic cross (HC) approximation Luo et al. (2017).

The remainder of the paper is organized as follows: Section 2 presents a function approximation method using HDMR-HC. Then a new meshless numerical scheme is proposed in Section 3 for solving the $(n+1)$-dimensional SGE using the HDMR-HC approximation. In Section 4, we provide several examples with a comparative numerical error analysis. Section 5 summarizes the relevant results.

2 Function approximation using HDMR-HC and the partition of unity

2.1 HDMR-HC approximation

Let $n\geqslant 2$, $\lambda\in\mathbb{R}_{+}$ and $f(\boldsymbol{x})\in\mathcal{W}_{1}^{s}(\mathbb{T}^{n}_{\lambda}),s\in\mathbb{N}_{0}$ be a $n$-variate function which is $\frac{1}{\lambda}$-periodic in each variable, where $\mathbb{T}^{n}_{\lambda}$ is the $n$-torus given by Luo et al. (2017)

[TABLE]

and the function space $\mathcal{W}_{1}^{s}(\mathbb{T}^{n}_{\lambda})$ is defined by

[TABLE]

with the norm

[TABLE]

where $r\in\mathbb{N}_{0}^{n}$ denotes a $n$-dimensional multi-index with the norm

[TABLE]

and the $r$-th order mixed derivative is given by

[TABLE]

We consider the multivariate Fourier series of $f(\boldsymbol{x})$

[TABLE]

where the Fourier coefficients $c_{m}$ are defined by

[TABLE]

For a nonempty set $\{j_{1},\cdots,j_{v}\}\subset\{1,\cdots,n\}$, let

[TABLE]

and

[TABLE]

Note that $\Lambda_{0}$ is the set consisting of the $n$-dimensional zero vector. Then $\mathbb{Z}^{n}$ can be decomposed into the following form

[TABLE]

and we refer to this as a HDMR decomposition of $\mathbb{Z}^{n}$; then a multiple Fourier series can be decomposed with an HDMR structure

[TABLE]

If let

[TABLE]

that is,

[TABLE]

then we define the Fourier HDMR-HC partial sum up to $K_{\lambda}$-th order Luo et al. (2017)

[TABLE]

Suppose $M$ is the number of points $\boldsymbol{x}\in\mathbb{T}_{\lambda}^{n}$, which depends on both $\lambda$ and $n$, and

[TABLE]

then Eq.(14) can be rewritten as

[TABLE]

Moreover, it follows that, for any $f\in\mathcal{W}_{1}^{s}(\mathbb{T}^{n}_{\lambda})$, $s>\gamma+p+1$ and $K_{\lambda}\geqslant\delta$, the bound

[TABLE]

holds pointwise almost everywhere, where the constant $C_{n,\lambda,s,p,\gamma,\delta}$ depends on $n,\lambda,s,p$ and $\delta$.

From the definition of $\mathbb{T}^{n}_{\lambda}$, we find that $\mathbb{T}^{n}_{\lambda}$ is smaller as $\lambda$ becomes larger. The convergence rate of the Fourier HDMR-HC partial sums $S_{K_{\lambda}}$ is closely related to the value of $\lambda$. In particular, for a given accuracy, when $\lambda$ is large, $S_{K_{\lambda}}$ generally is well approximated by a low order truncated HDMR, and this is the basic starting point of the Fourier HDMR-HC approximation. For example, if there is a function of $n=5$ variables, and suppose $\lambda=1/\pi$ and the desired accuracy is $\varepsilon=\mathbf{O}(K_{\lambda}^{-(s-p-\gamma-1)})$, then a $K_{\lambda}$-th order Fourier HDMR-HC partial sum $S_{K_{\lambda}}$ is just a $2$nd order truncated HDMR of the function when $K_{\lambda}\leqslant 31$. Therefore, we expect that a low order truncated HDMR-HC can be used to effectively capture the behavior of a high-dimensional function and its derivatives.

2.2 Partition of unity

In this subsection, we will discuss how to approximate a function via HDMR-HC. Suppose $u(\boldsymbol{x})$ is a function defined on $\mathbb{R}$, where $\boldsymbol{x}\in\Omega\subseteq\mathbb{R}^{n}$. First, we divide the domain $\Omega$ into $D\in\mathbb{N}$ subdomains, denoted as $\Omega_{j},~{}(j=1,2,\cdots,D)$, thus $\Omega=\bigcup_{j=1}^{D}\Omega_{j}$. In each subdomain $\Omega_{j}$, let $\boldsymbol{x}_{j}$ be the centre of the region, and let $\boldsymbol{\chi}_{j}$ denote all $M_{j}$ nodes in $\Omega_{j}$, i.e., $\sum_{j=1}^{D}M_{j}=N$. Further, suppose $u(\boldsymbol{x}),\boldsymbol{x}\in\Omega_{j}$ can be represented as $u_{j}$, then we have

[TABLE]

where

[TABLE]

is the characteristic function of $u_{j}$ satisfying

[TABLE]

From Eq.(15), $u_{j}$ at any point $\boldsymbol{x}\in\Omega_{j}$ can be approximated as

[TABLE]

where $\boldsymbol{F_{\lambda}}_{j}^{\mathrm{T}}(\boldsymbol{x})=[F_{\lambda,1},F_{\lambda,2},\cdots,F_{\lambda,M_{j}}]$ with the same $\lambda$ for all nodes and $c_{m_{j}}$ is the $m_{j}$-th unknown coefficient. Utilizing the values of these $M_{j}$ nodes, there are $M_{j}$ equations with one for each node, then we have following matrix form

[TABLE]

where

[TABLE]

is the vector of function values at the $M_{j}$ nodes, and

[TABLE]

is the vector of undetermined coefficients with

[TABLE]

where $F_{\lambda,im_{j}}=F_{\lambda,m_{j}}(\boldsymbol{\chi}_{i})$. Suppose $\boldsymbol{F}_{j}^{-1}$ exists (i.e., this condition can always be satisfied Powell (1992); Wendland (1998)), then $\boldsymbol{c}_{j}$ can be obtained by solving Eq.(20), i.e.

[TABLE]

From Eq.(16), we have

[TABLE]

where $\boldsymbol{U}=\bigcup_{j=1}^{D}\boldsymbol{U}_{j}$. Then

[TABLE]

Substitute both Eq.(21) and Eq.(22) back into Eq.(19), we have an approximation of function $u_{j}$:

[TABLE]

Furthermore, from Eqs.(16) and (18) we have

[TABLE]

Then $u(\boldsymbol{x})$ at any point $\boldsymbol{x}\in\Omega$ can be formally approximated as

[TABLE]

3 Solution for $(n+1)$-dimensional SGEs

Now, we present the numerical scheme for solving the $(n+1)$-dimensional SGE based on using the Fourier HDMR-HC approximation. Suppose the approximated function of the field function $u(\boldsymbol{x})$, ($\boldsymbol{x}=\{x_{1},x_{2},\cdots,x_{n}\}$) is formally denoted as

[TABLE]

where $M_{j}$ is the number of field nodes used in the selected domain, and $\boldsymbol{U_{j}}$ is the vector that collects the true nodal function values for these $M_{j}$ field nodes, and $\boldsymbol{x}_{j}$ is the centre of this selected region. Further, the derivatives of $u(\boldsymbol{x})$ at any point $\boldsymbol{x}$ can be approximated as

[TABLE]

where $x_{l}$ denotes one element of $\boldsymbol{x}=\{x_{1},x_{2},\cdots,x_{n}\}$.

In this paper, the time derivatives are approximated by the time-stepping method and we have the following approximation:

[TABLE]

[TABLE]

where $\tau$ is the time step, and $u^{(k)}$ is the approximate value of $u(\boldsymbol{x},t)$ at $(\boldsymbol{x},t_{k})$, $t_{k}=k\tau$. Moreover the Crank-Nicolson scheme is used to approximate $u$ at three respective times as

[TABLE]

To manage the nonlinearity, a Quasilinearization Method (QLM) is adopted. The QLM is very effective for dealing with the nonlinear aspects of the SGE and other PDEs. In this fashion the nonlinear term in Eq. (1) can be represented as

[TABLE]

Thus, Eq.(1) can be discretized as

[TABLE]

where $\mu=\tau^{-2}$, $\eta=\frac{\beta}{2\tau}$ and $\varphi^{(k)}=\psi(\boldsymbol{x})\cos(u^{(k)})$.

Suppose $N$ field nodes are denoted as $\boldsymbol{\chi}=\{\boldsymbol{\chi}\}_{1}^{N}$, where $N$ is determined by both $D$ and $M_{j}$. Then from Eq.(26) we have an approximation of the field variable $u$ according to the HDMR-HC approximation

[TABLE]

where $\phi_{j}(\boldsymbol{x})$ depends on $J_{j},\boldsymbol{F_{\lambda}}_{j}^{\mathrm{T}}(\boldsymbol{x})$ and $\boldsymbol{F}_{j}^{-1}$. The derivatives of $\hat{u}$ can be approximated as

[TABLE]

Therefore, for any point $\boldsymbol{x}_{i}$, the approximation in Eq.(33) can be written as

[TABLE]

Let $A_{ij}=\frac{1}{3}\Delta\phi_{j}(\boldsymbol{x}_{i})$, $B_{ij}=\phi_{j}(\boldsymbol{x}_{i})$, and $E_{ij}^{(k)}=A_{ij}-(\eta+\mu+\varphi^{(k)})B_{ij}$, $G_{ij}^{(k)}=-A_{ij}-(2\mu+\varphi^{(k)})B_{ij}$, $H_{ij}=A_{ij}-(\eta-\mu)B_{ij}$, $C_{i}^{(k)}=\psi(\boldsymbol{x}_{i})\sin(u^{(k)})$, then Eq.(36) can be re-written as

[TABLE]

For all $N$ field nodes $\boldsymbol{\chi}$ we have following matrix form:

[TABLE]

where $\boldsymbol{C}^{{}^{\prime}(k)}=[C_{1}^{(k)},C_{2}^{(k)},\cdots,C_{N}^{(k)}]^{\mathrm{T}}$.

In using a meshless strong method to solve the PDE, the solution can be unstable if there is a derivative boundary condition, so the fictitious points method is used to impose derivative boundary conditions Liu and Gu (2005). Suppose there are $N_{b}$ nodes on the boundary, then along the derivative boundaries, another $N_{b}$ fictitious points are added outside of the domain. Two sets of equations are established at each derivative boundary node: one for the derivative boundary condition, and the other for the governing equation. With the $N_{b}$ additional degrees of freedom, $\{u_{(N+1)},u_{(N+2)},\cdots,u_{(N+N_{b})}\}$, added into the system, then Eq.(36) can be re-written as

[TABLE]

and for a node at $\boldsymbol{x}_{i_{b}}$, that is, on the derivative boundary, the derivative boundary conditions have the form

[TABLE]

Assembling Eqs.(39) and (40) for the corresponding nodes, the discretized global system equation becomes

[TABLE]

where $\boldsymbol{E}$, $\boldsymbol{G}$, $\boldsymbol{B}$ and $\boldsymbol{H}$ are ${(N+N_{b})\times(N+N_{b})}$ matrices, $\boldsymbol{\hat{u}}=[\hat{u}_{1},\hat{u}_{2},\cdots,\hat{u}_{(N+N_{b})}]^{\mathrm{T}}$, $\boldsymbol{C}^{(k)}=[C_{1}^{(k)},C_{2}^{(k)},\cdots,C_{N}^{(k)},3w(\boldsymbol{x}_{N+1},t),3w(\boldsymbol{x}_{N+2},t),\cdots,3w(\boldsymbol{x}_{N+N_{b}},t)]^{\mathrm{T}}$. At the first time level, i.e. $k=0$, we adopt the following:

[TABLE]

and

[TABLE]

where $\boldsymbol{v_{1}}$ and $\boldsymbol{v_{2}}$ are the initial conditions for all nodes $\boldsymbol{\chi}$ introduced in Eqs.(2) and (3).

4 Numerical experiments

In this section, the proposed meshless numerical scheme is applied to several examples to show the efficiency and accuracy for the $(n+1)$-dimensional SGE. As mentioned in the previous section, to approximate the time derivatives we use a finite difference method, so an iterative scheme is employed to reach the final time $t$. In order to test the performance of the numerical solution, we use the $L_{\infty}$ error and root-mean-square (RMS) error norms defined as

[TABLE]

and

[TABLE]

where $N$ is the number of nodes, $f(x_{i})$ is the exact solution, and $\hat{f}(x_{i})$ is the numerical solution. To assess both the stability and the solution accuracy, we compute the condition number of the system matrix, which is defined as

[TABLE]

where $\kappa(\boldsymbol{F})$ depends on the parameter $\lambda$ and the number of nodes $N$.

4.1 Test problem for a $(2+1)$-dimensional SGE

The test problem for a $(2+1)$-dimensional SGE has the following form Djidjeli et al. (1995); Dehghan and Shokri (2008b); Jiwari et al. (2012b)

[TABLE]

where $\Omega=[-7,7]\times[-7,7]$, and the initial conditions are

[TABLE]

and the Neumann boundary condition is

[TABLE]

The analytic solution of this problem is:

[TABLE]

In this example, both the proposed HDMR-HC meshless method and the radial basis function (RBF) method in Dehghan and Shokri (2008b) are used to numerically solve the equation. Since the field nodes of the RBF method in Dehghan and Shokri (2008b) is the Sobol sequence with $N=3249$, then we adopt the same total number of field nodes for the HDMR-HC meshless method ($D=49$). The time step $\tau$ is set to $0.001$. The results of the two different measures of error are presented in Table 1. We see that with the same number of field nodes, the errors of the proposed HDMR-HC meshless method are smaller than those of the RBF method. The condition numbers at particular times are also listed in Table 1.

4.2 Test problem for a $(5+1)$-dimensional SGE

To further test the proposed HDMR-HC scheme, we choose a $(5+1)$-dimensional example, which involves all the implementation issues explained in the previous subsection. The exact solution has the form of Eq.(5)

[TABLE]

where $C=1$, $b=1$, $\Omega=[-6,6]\times[-6,6]\times\cdots\times[-6,6]\subseteq\mathbb{R}^{5}$ and

[TABLE]

The initial conditions are

[TABLE]

and the Neumann boundary conditions are

[TABLE]

Both the proposed HDMR-HC method and the RBF method are used to solve the equation. The Sobol sequence is chosen as the field nodes with a total number of $N=2^{16}$ for the HDMR-HC method ($D=243$ and we use the HDMR approximation up to order $3$) and $N^{\prime}=2^{18}$ for the RBF method. In this case the time step is chosen as $\tau=0.2$. Table 2 presents $L_{\infty}$, RMS errors and the condition number $\kappa$ at some selected times $t$.

5 Conclusions

In this paper, we propose a new meshless solution method for high-dimensional sine-Gordon equations. First, we present a function approximation using the HDMR-HC decomposition. Then we divided the whole domain into several subdomains with the help of the partition of unity, and obtain a function approximation at any random point in each subdomain. Hence, we develop a numerical procedure for the high-dimensional SGEs by a meshless strong solution method. The time-stepping method is used to approximate the time derivatives of SGEs, and a quasilinearization scheme is performed to treat the nonlinearity of the equation. Finally, to demonstrate the accuracy of the proposed method with two numerical experiments. The examples suggest that the proposed procedure is attractive for solving high-dimensional SGEs.

Acknowledgements.

The authors X.X. and X.L. acknowledge support from the National Science Foundation (Grant No. CHE-1763198), and H.R. acknowledges support from the Templeton Foundation (Grant No. 52265).