The Spectral Method for Solving Sine-Gordon Equation Using a New Orthogonal Polynomial

Abstract

We have presented a numerical scheme for solving one-dimensional nonlinear sine-Gordon equation. We apply the spectral method with a basis of a new orthogonal polynomial which is orthogonal over the interval [0,1] with weighting function one. The results show the accuracy and efficiency of the proposed method.

1. Introduction

The sine-Gordon equation is a nonlinear hyperbolic partial differential equation (PDE). It was originally considered in the nineteenth century in the course of surface of constant negative curvature. This equation attracted a lot of attention in the 1970s due to the presence of soliton solutions [1, 2]. The sine-Gordon equation appears in a number of physical applications [3–5], including applications in the chain of coupled pendulums and modelling the propagation of transverse electromagnetic (TEM) wave on a superconductor transmission system. Consider the one-dimensional nonlinear sine-Gordon equation𝜕²𝑢(𝑥,𝑡)𝜕𝑡²=𝜕²𝑢(𝑥,𝑡)𝜕𝑥²−sin(𝑢(𝑥,𝑡)),(𝑥,𝑡)∈[0,1]×𝑡₀,𝑇,(1.1) with the initial conditions𝑢𝑥,𝑡₀=𝑓⁰(𝑥),𝜕𝑢𝜕𝑡𝑥,𝑡⁰=𝑓₁(𝑥),𝑥∈[0,1],(1.2) and the boundary conditions𝑢(0,𝑡)=𝑔₀(𝑡),𝑢(1,𝑡)=𝑔₁(𝑡),𝑡≥𝑡₀.(1.3) In the last decade, several numerical schemes have been developed for solving (1.1), for instance, high-order solution of one-dimensional sine-Gordon equation using compact finite difference and DIRKN methods [3]. The authors of [6] proposed a numerical scheme for solving (1.1) using collocation and radial basis functions. Also, the boundary integral equation approach is used in [7]. Bratsos has proposed a numerical method for solving one-dimensional sine-Gordon equation and a third-order numerical scheme for the two-dimensional sine-Gordon equation in [8, 9], respectively. Also, in [10], a numerical scheme using radial basis function for the solution of two-dimensional sine-Gordon equation is used. In addition, several authors proposed spectral methods and Fourier pseudospectral method for solving nonlinear wave equation using a discrete Fourier series and Chebyshev orthogonal polynomials [11–13].

In this paper, we have proposed a numerical scheme for solving (1.1) using spectral method with a basis of a new orthogonal polynomial. This polynomial is introduced by Chelyshkov in [14].

The outline of this paper is as follows. In Section 2, we introduce a new orthogonal polynomial. In Section 3, we apply spectral method for solving (1.1). In Section 4, we show the numerical results. Finally, a brief conclusion is given in Section 5.

2. The New Orthogonal Polynomial

In this section, we briefly introduce the new orthogonal polynomial, that is introduced in [14].

Let 𝑛 be a fixed whole number, and 𝑃_𝑛 is a sequence of polynomial, that is,𝑃_𝑛=𝑃_𝑛𝑘^0𝑘=𝑛,𝑃_𝑛𝑘≡𝑃_𝑛𝑘(𝑥)=^𝑛_𝑙=𝑘𝜏_𝑛𝑘𝑙𝑥^𝑙,(2.1) that satisfy the orthogonality relationship¹⁰𝑃_𝑛𝑘𝑃_𝑛𝑙𝑑𝑥=0,1𝑘≠𝑙,(𝑘+𝑙+1),𝑘=𝑙,(2.2) and standardizationsign𝜏_𝑛𝑘𝑛=(−1)^𝑛−𝑘.(2.3) The coefficient 𝜏_𝑛𝑘𝑙 of the polynomial 𝑃_𝑛𝑘 is defined by requirements (2.2), (2.3), and the Gram-Schmidt orthogonalization procedure (without normalization).The explicit definition of the polynomial 𝑃_𝑛𝑘 is as follows:𝑃_𝑛𝑘(𝑥)=^𝑛−𝑘_𝑗=0(−1)^𝑗⎛⎜⎜⎝𝑛−𝑗𝑘⎞⎟⎟⎠⎛⎜⎜⎝𝑛+𝑛𝑘−+𝑘1+𝑗⎞⎟⎟⎠𝑥_𝑘+𝑗,𝑘=0,1,…,𝑛.(2.4) This yields Rodrigue’s type representation as follows:𝑃_𝑛𝑘(𝑥)=1(𝑛−𝑘)!𝑥1^𝑘+1𝑑^𝑛−𝑘𝑑𝑥^𝑛−𝑘𝑥^𝑛+𝑘+1(1−𝑥)^𝑛−𝑘,𝑘=0,1,…,𝑛,(2.5)

and the orthogonality relationship (2.2) is confirmed by applying (2.5). It also follows from (2.5) that¹⁰𝑃_𝑛𝑘(𝑥)𝑑𝑥=¹⁰𝑥^𝑛𝑑𝑥=1𝑛+1.(2.6)This polynomial can be connected to a fixed set of the Jacobi polynomials 𝑃^{𝑛(𝛼,𝛽)}(𝜉) [15], that is,𝑃_𝑛𝑘=𝑥^𝑛𝑃^{(2𝑛,0)𝑘−𝑛}(1−2𝑥).(2.7)

3. The Proposed Method

In this section, the introduced orthogonal polynomial is applied as the basis for spectral method to solve sine-Gordon equation.

The approximate solution 𝑢_𝑁(𝑥,𝑡) to the exact solution 𝑢(𝑥,𝑡) can be written in the form𝑢(𝑥,𝑡)≅𝑈^𝑁(𝑥,𝑡)=∑^𝑁_𝑛=0𝑎_𝑛(𝑡)𝑃_𝑛0(𝑥),(3.1)

where 𝑎_𝑛 is the time-dependent quantities which must be determined. In this paper, we apply orthogonal polynomial 𝑃_𝑛𝑘 for 𝑘=0. From (3.1), we get𝑈^𝑁𝑡𝑡=^𝑁_𝑛=0𝑈𝑎̈_𝑛(𝑡)𝑃_𝑛0(𝑥),_𝑁𝑥𝑥=^𝑁_𝑛=0𝑎_𝑛(𝑡)𝑃̈_𝑛0(𝑥).(3.2)

Using (3.2) and substituting in (1.1), we obtain^𝑁∑_𝑛=0𝑎̈_𝑛(𝑡)𝑃_𝑛0𝑥_𝑖=∑^𝑁_𝑛=0𝑎_𝑛(𝑡)𝑃̈_𝑛0𝑥_𝑖−sin𝑈^𝑁𝑥_𝑖,𝑡,𝑖=1,2,…,𝑁−1,(3.3)

where 𝑎̈_𝑛 and 𝑃̈_𝑛0 denote the second-order derivatives of 𝑎_𝑛 and 𝑃_𝑛0 with respect to 𝑡 and 𝑥, respectively, wherein the collocation points are as 𝑥_𝑖=𝑖/𝑁,𝑖=0,1,…,𝑁.

Equation (13) is rewritten as follows:𝑀𝐴̈(𝑡)=𝐾𝐴(𝑡)−𝑆₁.(3.4)

For displaying this equation in matrix form, for 𝑖=1,2,…,𝑁−1, and 𝑛=0,1,…,𝑁, we define 𝐴(𝑡)=𝑀𝑎=₀(𝑀𝑡),_in𝑎₁(𝑡),…,𝑎_𝑁(𝑡)^𝑇,𝐾=𝐾_in=𝑃_𝑛0𝑥_𝑖,𝑆₁=𝑃̈_𝑛0𝑥_𝑖,=sin𝑈₁𝑈,_𝑖sin𝑈₂,…,sin𝑈_𝑁−1^𝑇,=𝑈^𝑁𝑥_𝑖,𝑡.(3.5)

The matrices 𝑀 and 𝐾 are positive definite matrices of order (𝑁−1)×(𝑁−1). The system given in (3.4) consists of (𝑁−1) equations in (𝑁+1) unknowns. To obtain a unique solution for this system, we need two additional equations. For this, we add two boundary conditions, which are given by (1.3), in the following forms:𝑈^𝑁(0,𝑡)=^𝑁_𝑛=0𝑈_𝑁𝑎_𝑛(𝑡)𝑃_𝑛0(0)=𝑔₀(𝑡),𝑡≥𝑡₀,(1,𝑡)=^𝑁_𝑛=0𝑎_𝑛(𝑡)𝑃_𝑛0(1)=𝑔₁(𝑡),𝑡≥𝑡₀.(3.6)

From (3.6), we get^𝑁_𝑛=0𝑎̈_𝑛(𝑡)𝑃_𝑛0(0)=𝑑²𝑔₀(𝑡)𝑑𝑡²,𝑡≥𝑡₀,_𝑁_𝑛=0𝑎̈_𝑛(𝑡)𝑃_𝑛0(1)=𝑑²𝑔₁(𝑡)𝑑𝑡²,𝑡≥𝑡₀.(3.7)

Therefore, by using (3.6), and (3.7), a new system is obtained as follows:𝑀𝐴̈(𝑡)=𝐹(𝑡,𝐴(𝑡)),𝐹(𝑡,𝐴(𝑡))=𝐾𝐴(𝑡)−𝑆+𝐺(𝑡),(3.8)

where 𝐴(𝑡), 𝐹(𝑡,𝐴(𝑡)), 𝑆, and 𝐺(𝑡) are vectors of order (𝑁+1) with (𝑁+1) components. The matrices 𝑀 and 𝐾 are of order (𝑁+1)×(𝑁+1), which are defined as follows:𝑀=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝𝑝⁰⁰𝑝(0)𝑝₁₀(0)⋯𝑝_𝑁0(0)₀₀𝑝⋮𝑥₁𝑝₁₀⋮𝑥₁⋱⋯𝑝_𝑁0⋮₀₀𝑥₁𝑝₀₀𝑥_𝑁−1𝑝₁₀𝑥_𝑁−1⋯𝑝_𝑁0𝑥_𝑁−1𝐾=⎛⎜⎜⎜⎜⎜⎜⎜⎜⎝𝑝̈₀₀0(1)𝑝0₁₀(1)⋯⋯𝑝0_𝑁0(1)⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠,(3.9)𝑝̈₀₀⋮𝑥₁𝑝̈₁₀⋮𝑥₁⋱⋯𝑝̈_𝑁0⋮𝑥₁𝑥0_𝑁−1𝑝̈₁₀𝑥0_𝑁−1⋯⋯𝑝̈_𝑁00𝑥_𝑁−1⎞⎟⎟⎟⎟⎟⎟⎟⎟⎠.(3.10) Also, the vectors 𝑆 and 𝐺(𝑡) are defined as follows:𝑆=𝐺0,(s𝑡i)n=𝑈₁𝑑₂,sin𝑈₂,…,sin𝑈_𝑁−1,0^𝑇,𝑔₀(𝑡)𝑑𝑡²,0,…,0,𝑑²𝑔₁(𝑡)𝑑𝑡²^𝑇.(3.11) The nonlinear second-order system of ODEs (3.8) can be solved numerically using the fourth-order A-stable DIRKN method [3, 16]. From the initial conditions (1.2), we determine the initial vectors 𝐴₀ and 𝐴̇₀. Using (3.1), we get𝑈^𝑁𝑥_𝑖,0=^𝑁_𝑛=0𝜕𝑈_𝑁𝑎_𝑛𝑡₀𝑃_𝑛0𝑥_𝑖,𝑖=0,1,…,𝑁,𝜕𝑡𝑥^𝑖,0=^𝑁_𝑛=0𝑎̇_𝑛𝑡₀𝑃_𝑛0𝑥_𝑖,𝑖=0,1,…,𝑁.(3.12)

Equation (3.12) is system of (𝑁+1) equations in (𝑁+1) unknowns. These equations can be written in the matrix form 𝑀𝐴₀=𝑏₁ and 𝑀𝐴̇₀=𝑏₂, where 𝑏𝑏₁₂=𝑓₀𝑥₀,…,𝑓₀𝑥_𝑁^𝑇,=𝑓₁𝑥₀,…,𝑓₁𝑥_𝑁^𝑇.(3.13)

3.1. The DIRKN Method

For solving (3.8), we use DIRKN method [3, 16]:𝑀𝑄_𝑖=𝑀𝑀𝐴𝐴_𝑛𝑛+1+𝑐_𝑖Δ𝑡𝑀𝐴̇_𝑛+Δ𝑡∑^2𝑠_𝑗=1𝑎_𝑖𝑗𝐹𝑡_𝑛+𝑐_𝑗Δ𝑡,𝑄_𝑗,𝑖=1,…,𝑠,=𝑀𝐴_𝑛+𝑐_𝑖Δ𝑡𝑀𝐴̇_𝑛+Δ𝑡^2𝑠_𝑗=1𝑀𝐴̇_𝑛+1𝑏_𝑖𝐹𝑡_𝑛+𝑐_𝑗Δ𝑡,𝑄_𝑗,=𝑀𝐴̇_𝑛+Δ𝑡^𝑠_𝑗=1𝑏_𝑖𝐹𝑡_𝑛+𝑐_𝑗Δ𝑡,𝑄_𝑗,(3.14)

where 𝑄_𝑖 is the internal stages, 𝐴_𝑛+1 and 𝐴̇_𝑛+1 represent approximations to 𝐴(𝑡_𝑛+1) and 𝐴(̇𝑡_𝑛+1), respectively. Also we defineΔ𝑡=𝑡₀+𝑗Δ𝑡,𝑗=0,1,…,𝑀,(3.15)where Δ𝑡=(𝑇−𝑡₀)/𝑀. The RKN method can be denoted in 𝐵𝑢𝑡𝑐ℎ𝑒𝑟^,𝑠 notation by the table of coefficients:𝑐₁𝑎₁₁⋯𝑎_1𝑠⋮𝑐⋮⋱⋮_𝑠𝑎_𝑠1⋯𝑎_𝑠𝑠𝑏₁⋯𝑏_𝑠𝑏₁⋯𝑏_𝑠(3.16)

In the DIRKN method, 𝑎_𝑖𝑗=0 as (𝑖<𝑗) and 𝑎_𝑖𝑖 are equal. Also, in the above table of coefficients, 𝑠 is the number of stages. We use the four-stage A-stable DIRKN method [16], with algebraic fourth order, which is defined by the following table of coefficients for solving all experiments.

Using Table 1, we can calculate 𝑄_𝑖, 𝑖=1,…,4, separately from (3.14) with Newton-Raphson iteration method. We iterate until the following criterion is satisfied:‖‖𝑄^𝑖𝑘+1−𝑄^𝑘𝑖‖‖_∞<10⁻¹²,(3.17) where 𝑄^𝑘𝑖 is the value of 𝑄_𝑖 at the kth iteration of the Newton-Raphson method.

Table 1 

Four-stage A-stable DIRKN method.

4. Numerical Experiments

In this section, we present the results of numerical experiments using the method introduced in Section 3. The 𝐿₂ and 𝐿_∞ error norms are defined as follows:𝐿₂=‖‖𝑢−𝑈^𝑁‖‖=^𝑁_𝑖=0𝑢_𝑖−𝑈^𝑁𝑖²^1/2𝐿_∞,=‖‖𝑢−𝑈^𝑁‖‖=max_{0≤𝑖≤𝑁}||𝑢_𝑖−U^𝑁𝑖||.(4.1)

We choose Δ𝑡=1/1000, for solving all experiments.

The computations associated with the experiments discussed above were performed in Maple 13 on a PC with a CPU of 2.4 GHZ.

Experiment 1. In this experiment [3, 6, 10], we consider (1.1) with 𝑡₀=0 and the initial conditions 𝑓𝑓₀(𝑥)=0,₁(𝑥)=4sech(𝑥),(4.2) with the boundary conditions (1.3), which can be obtained from the following exact solution: 𝑢(𝑥,𝑡)=4arctan(sech(𝑥)t).(4.3) In Table 2, the approximate solution for several final times (𝑇) with different locations is shown. Figure 1 shows two type plots at 𝑡=0.5 of Experiment 1.

Table 2 

The approximate solution for 𝑁=16 of Experiment 1.

(a)

(b)

Figure 1 

(a) The exact and approximate solutions for𝑁=16 are shown. (b) The 𝐿∞ error norm for 𝑁=16 is shown.

Experiment 2. In this experiment [3, 8], we consider (1.1) with 𝑡₀=0 and the initial conditions 𝑓𝑓₀(𝑥)=0,(4.4)₁(𝑥)=4√1+𝑐²sech𝑥√1+𝑐²,(4.5) with the boundary conditions (1.3), which can be obtained from the following exact solution: 𝑢(𝑥,𝑡)=4arctan𝑐⁻¹sin𝑐𝑡√1+𝑐²sech√1𝑥+𝑐²,𝑐=1.(4.6)
In Tables 3 and 4, we show the 𝐿₂ and 𝐿_∞ error norms, for several final times (𝑇). As we see from Tables 3 and 4, for the given values of 𝑁, the 𝐿_∞ error norm is less than the 𝐿₂ error norm. Figure 2 shows four type plots at 𝑡=0.5 of Experiment 2.

Table 3 

The 𝐿₂ error norm for different values of 𝑁 of Experiment 2.

Table 4 

The 𝐿_∞ error norm for different values of 𝑁 of Experiment 2.

(a)

(b)

(c)

(d)

Figure 2 

(a) The exact and approximate solutions for 𝑁=16 are shown. (b) The 𝐿∞ error norm for 𝑁=16 is shown. (c) The 𝐿2 error norm for different values of 𝑁 is shown. (d) The 𝐿∞ error norm for different values of 𝑁 is shown.

Experiment 3. In experiment [3, 9], we consider (1.1) with 𝑡₀=0 and the initial conditions 𝑓₀(𝑥)=4arctan𝑐sinh𝑥𝑓√1−𝑐²,(4.7)₁(𝑥)=0,(4.8) with the boundary conditions (1.3), which can be obtained from the following exact solution: 𝑢(𝑥,𝑡)=4arctan⎛⎜⎜⎝𝑐sinh𝑥/√1−𝑐²cosh𝑐𝑡/√1−𝑐²⎞⎟⎟⎠,𝑐=0.5.(4.9)
In Table 5, we represent the 𝐿₂ and 𝐿_∞ error norms for different values of 𝑁 at 𝑡=0.5. With increasing 𝑁, we obtain better results, and these error norms decrease five orders in magnitude. Figure 3 shows four-type plots at 𝑡=0.5 of Experiment 3.

Table 5 

Error norms for different values of 𝑁 of Experiment 3.

(a)

(b)

(c)

(d)

Figure 3 

(a) The exact and approximate solutions for 𝑁=16 are shown. (b) The 𝐿∞ error norm for 𝑁=16 is shown. (c) The 𝐿2 error norm for different values of 𝑁 is shown. (d) The 𝐿∞ error norm for different values of 𝑁 is shown.

5. Conclusion

In this paper, we have presented the spectral method for solving one-dimensional nonlinear sine-Gordon equation using a new orthogonal polynomial. Numerical experiments show that the spectral method is an efficient one, and the results for value of 𝑁=16 is more accurate than the other values of 𝑁. The proposed method has higher order of accuracy even for small values of 𝑁 (i.e., 𝑁=4). We obtain better results when we choose small value of Δ𝑡, in comparison to those obtained in the literature.