In this study, we propose an unconditionally stable temporally second-order accurate scheme for a parabolic sine-Gordon equation. The proposed scheme is based on an operator splitting method. We solve linear and nonlinear equations using a Fourier spectral method and a closed-form solution, respectively. The proposed numerical method is temporally second-order accurate and unconditionally stable. To verify the superior efficiency and accuracy of the proposed scheme, we conduct various numerical tests. Computational tests validate the accuracy, efficiency, and simplicity of the proposed scheme.

The sine-Gordon (SG) equation is a nonlinear hyperbolic partial differential equation (PDE) that plays an essential role in the modeling of many interesting problems in physical sciences, including plasma physics, solid-state physics, quantum mechanics, and nonlinear optics,

2u(x,t)t2=κ2Δu(x,t)+sinu(x,t),xΩ,t>0.
(1)

Here, u(x, t) is an order parameter, ΩRd(d=1,2,3) is a domain, and κ2 is the diffusion constant.1 The SG equation(1) has been used in nonlinear optics to describe the propagation of light pulses in optical fibers, which play the most important role in transmitting signals from communication devices, such as telephones and the internet between different locations. The SG equation has been investigated by numerous researchers, including a method for finding a new exact solution to the coupled SG equation using a modified Kudryashov approach,5 a method for calculating a soliton solution to a generalized nonlinear Fokas–Lenells equation through an SG extension,19 and a modified cubic B-spline differential quadrature technique for a numerical simulation of a two-dimensional SG soliton.21 Other proposed approaches include a numerical method of a one-dimensional nonlinear SG equation by reducing the problem to a system of a first-order ordinary differential equation (ODE),22 a meshless method based on a singular boundary method for the numerical solutions to a nonlinear SG equation with Neumann boundary conditions,6 and a fourth-order energy-preserving scheme for the SG equation using the Fourier pseudo-spectral method.8 Moreover, a modified extended direct algebraic method that can handle both the SG equation and nonlinear equations in a direct and concise way,14 a numerical method for a space–time fractional SG equation,2 and an efficient numerical scheme for solving the two-dimensional space–time fractional SG equation within the non-rectangular domain18 have been proposed. Several studies have also been recently conducted on the second-order accurate methods4,20,23,24 to solve the SG equation. High-order accurate methods have also been studied.17,26

In this paper, we present an unconditionally stable temporally second-order accurate operator splitting scheme for a parabolic sine-Gordon equation (PSG),3 

u(x,t)t=κ2Δu(x,t)+sinu(x,t),xΩ,t>0.
(2)

Equation (2) can be derived from the gradient flow of the following energy functional:

E(u)=Ωκ22|u|2+cosudx.
(3)

Phase-field modeling has attracted considerable attention in the study of two-phase systems. Mass conservation of the phase field was considered.9 In particular, the PSG equation provides the global well-posedness and maximum principle of the classical solution.3 

This paper is organized as follows: In Sec. II, a numerical solution algorithm is introduced for an unconditionally stable scheme using the operator splitting method. In Sec. III, various numerical experiments are simulated to confirm the accuracy, stability, and simplicity of the proposed scheme and to show the dynamics of the PSG equation. Finally, some concluding remarks are presented in Sec. IV.

Now, we present an unconditionally stable temporally second-order accurate numerical solution algorithm for the PSG equation. We use the operator splitting method to solve the PSG equation (2). For simplicity of exposition, we describe the numerical solution algorithm for the PSG equation in the two-dimensional domain Ω = (Lx, Rx) × (Ly, Ry). Numerical solution algorithms in one- and three-dimensional spaces are described analogously. To discretize Eq. (2), let Ωh = {xi = Lx + (i − 0.5)h, yj = Ly + (j − 0.5)h)| 1 ≤ iNx, 1 ≤ jNy} be the discrete computational domain, where h = (RxLx)/Nx = (RyLy)/Ny is the uniform step size; in addition, Nx and Ny are the numbers of the grid points in the x- and y-directions, respectively. Let uijn be the numerical approximation of u(xi, yj, nΔt), where Δt is the time step. To obtain a second-order scheme, we require the following steps:11 

u(x,t+Δt)=NΔt/2LΔtNΔt/2u(x,t)+O(Δt3),xΩ,t>0,
(4)

where NΔt and LΔt are the nonlinear and linear solution operators of the PSG equation, respectively. Moreover, the operators NΔtu(x,t) and LΔtu(x,t) indicate the solutions of ut = sin u and ut = κ2Δu at t + Δt, respectively. The stability and convergence of the Strang-type splitting method (4) were studied for the Allen–Cahn,13,16,25 Cahn–Hilliard,15 and Swift–Hohenberg11 equations. It was theoretically proven that the second-order accuracy of the splitting method (4) is guaranteed when each operator has at least a second-order accuracy.7 Because each operator in Eq. (4) has a temporally exact accuracy, the numerical scheme is temporally second-order accurate. The splitting method (4) can be written in three sub-steps using intermediate values as follows:

u*=NΔt/2un,
(5)
u**=LΔtu*,
(6)
un+1=NΔt/2u**.
(7)

Before we solve the nonlinear operator NΔt/2, we consider the following equation:

v(x,t)t=sin(v(x,t)),
(8)

which has a closed-form solution. If sin(v(x, t)) = 0, then v(x, t) = 0. Provided that sin(v(x, t)) ≠ 0, we can write Eq. (8) in the form

vsinv=t,
(9)

where we omitted the argument for simplicity of the notation. By multiplying both the denominator and the numerator on the left-hand side of Eq. (9) by sin v, we have

sinvsin2vv=t.
(10)

By using the Pythagorean trigonometric identity, we obtain

sinv1cos2vv=t.
(11)

Then, by integrating both sides after using a partial fraction expansion on the left side, we have

12sinv1+cosv+sinv1cosvv=t,
(12)

which results in

12ln1cosv1+cosv=t+C(x),
(13)

where C(x) = 0.5[ln(1 − cos v(x, 0)) − ln(1 + cos v(x, 0))]. Because the cosine function is an even function, we obtain the solution to Eq. (13) as

v(x,t)=sgn(sin(v(x,0)))cos11e2t+2C(x)1+e2t+2C(x),
(14)

where sgn(ϕ) is the sign function, which is plus one if ϕ is positive, minus one if it is negative, and zero otherwise. Now, using the analytic solution form (14), we solve the nonlinear operator NΔt/2(5) as

uij*=sgn(sin(uijn))cos11eΔt+2Cijn1+eΔt+2Cijn,
(15)

where Cijn=0.5[ln(1cosuijn)ln(1+cosuijn)].

Next, to solve the linear operator LΔt, we consider the diffusion equation

v(x,t)t=κ2Δv(x,t).
(16)

We use the Fourier spectral method in (16). To solve Eq. (16) with the homogenous Neumann boundary condition, we use the discrete cosine transform. For the given data {vijm|i=1,,Nxandj=1,,Ny}, the discrete cosine transform is defined as follows:

v̂pqm=αpβqi=1Nxj=1Nyvijmcos(ξpπxi)cos(ηqπyj).

The inverse discrete cosine transform is

vijm=p=1Nxq=1Nyαpβqv̂pqmcos(ξpπxi)cos(ηqπyj),
(17)

where

αp=1Nx(p=1),2Nx(p2),βq=1Ny(q=1)2Ny(q2),
(18)

ξp = (p − 1)/Lx, and ηq = (q − 1)/Ly. Let us assume that

v(x,y,mΔt)=p=1Nxq=1Nyαpβqv̂pqmcos(ξpπx)cos(ηqπy).
(19)

Plugging Eq. (19) in Eq. (16) yields

dv̂pqdt=κ2ξpπ2+ηqπ2v̂pq.
(20)

Therefore, we obtain the following solution from Eq. (20):

v̂pqn+1=v̂pqneΔtκ2(ξpπ)2+(ηqπ)2.
(21)

We can then solve the linear operator LΔt(6) by obtaining the numerical solution uij** using Eqs. (17) and (21), i.e.,

ûpq**=ûpq*eΔtκ2(ξpπ)2+(ηqπ)2
(22)

and

uij**=p=1Nxq=1Nyαpβqûpq**cos(ξpπxi)cos(ηqπyj).
(23)

The final sub-step (7) can also be solved analogously as (15),

uijn+1=sgn(sin(uij**))cos11eΔt+2Cij**1+eΔt+2Cij**,
(24)

where Cij**=0.5[ln(1cosuij**)ln(1+cosuij**)]. More details regarding the procedure and the definitions of the notions are provided in Ref. 12. We note that the proposed numerical algorithm is an unconditionally stable scheme. First, we suppose that the nth solution satisfies |uijn|π. In the first sub-step (5), the solution u* is bounded by π from Eq. (15) for any time step Δt, i.e.,

|uij*|π.
(25)

The stability condition of a semi-analytical Fourier spectral method for the Allen–Cahn equation was presented in Ref. 12. In the second sub-step (6), substituting Eq. (22) in Eq. (23),

uij**=p=1Nxq=1Nyαpβqûpq*eΔtκ2(ξpπ)2+(ηqπ)2×cos(ξpπxi)cos(ηqπyj).
(26)

The uij** in Eq. (26) is a solution of the heat equation for an initial condition uij*. Because it is well known that the heat equation with homogeneous Neumann boundary condition satisfies the maximum principle, we have

maxij|uij**|maxij|uij*|π.
(27)

Finally, in the last sub-step (7), the following inequality is established because |uij**| is bounded by π from Eq. (15) for any time step Δt,

|un+1|π.
(28)

Therefore, the proposed method is unconditionally stable regardless of the time step Δt and the maximum norm of each solution is bounded by π, which implies the boundedness of the numerical solutions.

We use the second-order operator splitting scheme (6) to solve the PSG equation. Unless otherwise mentioned, we use numerical parameters as Nx = 128, Ny = 128, and Nz = 128 for grid points of the spatial discretization, κ = 0.2 and Ω = (−π, π)d, where d = 1, 2, 3. Each nonlinear and linear operator is solved using (15) and (23), respectively.

We investigate the convergence and stability of the proposed scheme with respect to the time step. For the test, we take ϕ(x, 0) = cos(x) as the initial condition. First, we evaluate the errors and compare with the reference solution and convergence rate. Let log2(eNxNt/eNx2Nt) be the temporal convergence rate. Table I lists the errors with a reference solution and temporal convergence rates for the proposed method in the one-dimensional space with various values of Δt and a fixed space step size h = 2π/Nx at t = 1. According to Table I, we confirm that this scheme guarantees a second-order convergence with respect to time step. Figure 1 shows the evolution of the energy functional (3) with reference energy Eref (solid line) and different time steps Δt = 1 (circle), 2−2 (triangle), and 2−4 (dot). The energy monotonically decreases with time, and the numerical energies fit well with the reference energy for all values of Δt. In Fig. 2, to check the stability condition, we use random perturbations between −0.1 and 0.1 as the initial condition in each case. Figures 2(a) and 2(b) show the dynamics of the 1D PSG equation with Δt = 0.1 and Δt = 10, respectively. In Fig. 2(b), the solution does not blow up and shows similar dynamics as the solution illustrated in Fig. 2(a), which uses a relatively small time step compared to Fig. 2(b). From the stable results using a time step 100 times larger than 0.1, we can confirm the unconditional stability of the proposed scheme.

TABLE I.

l2-norm errors and convergence rates for the proposed method with fixed h and various Δt. Here, h = 2π/210, Δt = 2−4, and Δtref = 2−15 are used with t = 1.

CaseΔtRateΔt/2RateΔt/4RateΔt/8
l2-error 1.036 88 × 10−5 1.9997 2.592 80 × 10−6 1.9999 6.482 35 × 10−7 2.0000 1.620 57 × 10−7 
CaseΔtRateΔt/2RateΔt/4RateΔt/8
l2-error 1.036 88 × 10−5 1.9997 2.592 80 × 10−6 1.9999 6.482 35 × 10−7 2.0000 1.620 57 × 10−7 
FIG. 1.

Energy evolution with different time steps.

FIG. 1.

Energy evolution with different time steps.

Close modal
FIG. 2.

Temporal dynamics of the one-dimensional PSG equation with random perturbations between −0.1 and 0.1 for the initial condition, with different time steps (a) Δt = 0.1 and (b) Δt = 10.

FIG. 2.

Temporal dynamics of the one-dimensional PSG equation with random perturbations between −0.1 and 0.1 for the initial condition, with different time steps (a) Δt = 0.1 and (b) Δt = 10.

Close modal

To ensure that the PSG equation satisfies the boundedness of the numerical solutions, we consider the random perturbation between −1 and 1 as an initial condition with Nx = 28 within the 1D domain Ω ∈ (−π, π). Here, Δt = 0.05 is used. Figure 3(a) shows the temporal dynamics of the 1D PSG equation with a randomly perturbed initial condition. Figure 3(b) shows the maximum and minimum values of the solution over time t = 0 to t = 10. In Figs. 3(a) and 3(b), we can see that the values of the solution are bounded by π. This shows that the PSG equation satisfies the boundedness of the numerical solutions.

FIG. 3.

(a) Temporal dynamics of the one-dimensional PSG equation with random perturbation between −1 and 1 and (b) maximum and minimum values of the solution with Nx = 28 and Δt = 0.05.

FIG. 3.

(a) Temporal dynamics of the one-dimensional PSG equation with random perturbation between −1 and 1 and (b) maximum and minimum values of the solution with Nx = 28 and Δt = 0.05.

Close modal

We consider a linear stability analysis of the PSG equation. At u(x, t) = 0, the nonlinear term sin u(x, t) can be linearized using the Taylor expansion as sin u(x, t) ≈ u(x, t). Therefore, the linearized PSG equation is as follows:

u(x,t)t=κ2Δu(x,t)+u(x,t).
(29)

For the positive integers Ki and xRd where d = 1, 2, 3 is the space dimension, we can obtain

u(x,t)=α(t)i=1dcos(Kixi),
(30)

where α(t) is the amplitude of u(x, t). If d = 3, then x = (x1, x2, x3). Substituting u(x, t) in Eq. (30) into Eq. (29), we have

α(t)i=1dcos(Kixi)=κ2α(t)j=1dKj2i=1dcos(Kixi)+α(t)i=1dcos(Kixi).
(31)

By dividing i=1dcos(Kixi) on the both sides of Eq. (31), we obtain

α(t)=1κ2j=1dKj2α(t).
(32)

We can analytically solve the ordinary differential equation (32) and find the solution as follows: α(t) = α(0)eλt, where λ=1κ2j=1dKj2 is the analytic growth rate. We then define the numerical growth rate as

λ̃=1TloguNtα(0).
(33)

In a one-dimensional space Ω = (−π, π), we test the linear stability experiments with the initial condition u(x, 0) = 0.1 cos(Kx), where K = 0, 1, …, 10. As shown in Fig. 4, numerical growth rates λ̃ approximate well the analytic growth rates λ for each K. As mode K increases, the growth rate decreases. In particular, the growth rates have positive values only when K is less than or equal to 4.

FIG. 4.

Analytical and numerical growth rates of the PSG equation for different mode numbers K with κ=0.05, α(0) = 0.1, Δt = h2, and T = 500Δt.

FIG. 4.

Analytical and numerical growth rates of the PSG equation for different mode numbers K with κ=0.05, α(0) = 0.1, Δt = h2, and T = 500Δt.

Close modal

In Fig. 5, we simulate the 1D PSG equation to check how the numerical solutions grow depending on the mode number K. Only one-period solutions are shown. As it can be seen from the results in Figs. 4 and 5, the growth rate of the solution decreases as K increases. In particular, when mode K is greater than 4, the growth rates become negative, and the solutions are damped.

FIG. 5.

Snapshots of the evolution of the solutions for the one-dimensional PSG equation with the initial condition u(x, 0) = 0.1 cos(Kx). The values of K are shown in each figure (a) K = 1, (b) K = 4, (c) K = 5, and (d) K = 6.

FIG. 5.

Snapshots of the evolution of the solutions for the one-dimensional PSG equation with the initial condition u(x, 0) = 0.1 cos(Kx). The values of K are shown in each figure (a) K = 1, (b) K = 4, (c) K = 5, and (d) K = 6.

Close modal

Let us consider the following initial condition:

u(x,0)=π21+tanh1|x|0.1,x(10,10)d(d=1,2,3),
(34)

which implies that the initial values are between 0 and π. In this simulation, we consider the numerical parameters as Δt = 0.025, Nx = 128, Ny = 128, and Nz = 128. Figures 6(a)6(c) show the snapshots of the evolution of solutions of the PSG equation with the initial condition (34) in one-, two-, and three-dimensional spaces, respectively.

FIG. 6.

(a)–(c) are snapshots of the evolution of the solutions of the PSG equation with the initial condition (34) for 1D, 2D, and 3D spaces, respectively. Here, the isosurface plots are at a level of u = 0.1.

FIG. 6.

(a)–(c) are snapshots of the evolution of the solutions of the PSG equation with the initial condition (34) for 1D, 2D, and 3D spaces, respectively. Here, the isosurface plots are at a level of u = 0.1.

Close modal

Little research has been conducted on the specificity and dynamics of the parabolic sine-Gordon (PSG) equation. The main novelty of this study is to develop a fast and accurate numerical method for the PSG equation. We presented an unconditionally stable second-order accurate numerical method for the PSG equation. We solved linear and nonlinear equations using a Fourier spectral method and a closed-form solution. Because each split operator has a temporally exact accuracy, it can be extended to a high-order accurate method. To demonstrate the superior performance of the proposed scheme, we performed several computational experiments. We showed the second-order accuracy and unconditionally stability of the proposed scheme. The growth rate was analyzed using linear stability analysis, and it was confirmed that it was in good agreement with the simulation results. Through the computational tests, we confirmed the accuracy, efficiency, and simplicity of the proposed scheme. In future work, the proposed second-order operator splitting scheme will be extended to numerically solve the spatial fractional-order parabolic-type sine-Gordon equation.10,27

J. Kim was supported by the National Research Foundation (NRF), Korea, under project BK21 FOUR. The authors would like to thank the reviewers for their constructive and helpful comments regarding the revision of this paper.

The authors have no conflicts of interest to disclose.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

1.
K. K.
Ali
,
M. S.
Osman
, and
M.
Abdel-Aty
, “
New optical solitary wave solutions of Fokas–Lenells equation in optical fiber via sine-Gordon expansion method
,”
Alexandria Eng. J.
59
(
3
),
1191
1196
(
2020
).
2.
M.
Ahmadinia
and
Z.
Safari
, “
Analysis of local discontinuous Galerkin method for time-space fractional sine-Gordon equations
,”
Appl. Numer. Math.
148
,
1
17
(
2020
).
3.
X.
Cheng
,
D.
Li
,
C.
Quan
, and
W.
Yang
, “
On a parabolic sine-Gordon model
,”
Numer. Math. Theory Methods Appl.
14
(
4
),
1068
1084
(
2021
).
4.
Y.
Fu
,
W.
Cai
, and
Y.
Wang
, “
A linearly implicit structure-preserving scheme for the fractional sine-Gordon equation based on the IEQ approach
,”
Appl. Numer. Math.
160
,
368
385
(
2021
).
5.
K.
Hosseini
,
P.
Mayeli
, and
D.
Kumar
, “
New exact solutions of the coupled sine-Gordon equations in nonlinear optics using the modified Kudryashov method
,”
J. Mod. Opt.
65
(
3
),
361
364
(
2018
).
6.
Y.
Ji
, “
Meshless singular boundary method for nonlinear sine-Gordon equation
,”
Math. Probl. Eng.
2018
,
6460480
.
7.
T.
Jahnke
and
C.
Lubich
, “
Error bounds for exponential operator splittings
,”
BIT Numer. Math.
40
(
4
),
735
744
(
2000
).
8.
C.
Jiang
,
J.
Sun
,
H.
Li
, and
Y.
Wang
, “
A fourth-order AVF method for the numerical integration of sine-Gordon equation
,”
Appl. Math. Comput.
313
,
144
158
(
2017
).
9.
M.
Kwakkel
,
M.
Fernandino
, and
C. A.
Dorao
, “
A redefined energy functional to prevent mass loss in phase-field methods
,”
AIP Adv.
10
(
6
),
065124
(
2020
).
10.
H. G.
Lee
, “
A second-order operator splitting Fourier spectral method for fractional-in-space reaction-diffusion equations
,”
J. Comput. Appl. Math.
333
,
395
403
(
2018
).
11.
H. G.
Lee
, “
A semi-analytical Fourier spectral method for the Swift–Hohenberg equation
,”
Comput. Math. Appl.
74
(
8
),
1885
1896
(
2017
).
12.
H. G.
Lee
and
J.-Y.
Lee
, “
A semi-analytical Fourier spectral method for the Allen–Cahn equation
,”
Comput. Math. Appl.
68
(
3
),
174
184
(
2014
).
13.
H. G.
Lee
and
J. Y.
Lee
, “
A second order operator splitting method for Allen–Cahn type equations with nonlinear source terms
,”
Physica A
432
,
24
34
(
2015
).
14.
D.
Lu
,
A. R.
Seadawy
, and
M.
Arshad
, “
Elliptic function solutions and travelling wave solutions of nonlinear Dodd–Bullough–Mikhailov, two-dimensional Sine-Gordon and coupled Schrödinger–KdV dynamical models
,”
Results Phys.
10
,
995
1005
(
2018
).
15.
D.
Li
and
C.
Quan
, “
On the energy stability of Strang-splitting for Cahn–Hilliard
,” arXiv:2107.05349 (
2021
).
16.
D.
Li
,
C.
Quan
, and
T.
Tang
, “
Energy dissipation of Strang splitting for Allen–Cahn
,” arXiv:2108.05214 (
2021
).
17.
F.
Martin-Vergara
,
F.
Rus
, and
F. R.
Villatoro
, “
Padé schemes with Richardson extrapolation for the sine-Gordon equation
,”
Commun. Nonlinear Sci. Numer. Simul.
85
,
105243
(
2020
).
18.
F.
Mirzaee
,
S.
Rezaei
, and
N.
Samadyar
, “
Numerical solution of two-dimensional stochastic time-fractional Sine–Gordon equation on non-rectangular domains using finite difference and meshfree methods
,”
Eng. Anal. Boundary Elem.
127
,
53
63
(
2021
).
19.
A. K.
Mittal
, “
A stable time-space Jacobi pseudospectral method for two-dimensional sine-Gordon equation
,”
J. Appl. Math. Comput.
63
,
239
264
(
2020
).
20.
A. H.
Msmali
,
M.
Tamsir
, and
A. A. H.
Ahmadini
, “
Crank–Nicolson-DQM based on cubic exponential B-splines for the approximation of nonlinear Sine-Gordon equation
,”
Ain Shams Eng. J.
12
,
4091
4097
(
2021
).
21.
H. S.
Shukla
,
M.
Tamsir
, and
V. K.
Srivastava
, “
Numerical simulation of two dimensional sine-Gordon solitons using modified cubic B-spline differential quadrature method
,”
AIP Adv.
5
(
1
),
017121
(
2015
).
22.
H. S.
Shukla
and
M.
Tamsir
, “
Numerical solution of nonlinear sine–Gordon equation by using the modified cubic B-spline differential quadrature method
,”
Beni-Seuf Univ. J. Basic Appl. Sci.
7
(
4
),
359
366
(
2018
).
23.
F.
Ureña
,
L.
Gavete
,
A.
García
,
J. J.
Benito
, and
A. M.
Vargas
, “
Solving second order non-linear hyperbolic PDEs using generalized finite difference method (GFDM)
,”
J. Comput. Appl. Math.
363
,
1
21
(
2020
).
24.
J.-Y.
Wang
and
Q.-A.
Huang
, “
A family of effective structure-preserving schemes with second-order accuracy for the undamped sine–Gordon equation
,”
Comput. Math. Appl.
90
,
38
45
(
2021
).
25.
X.
Yang
, “
Error analysis of stabilized semi-implicit method of Allen–Cahn equation
,”
Discrete Contin. Dyn. Syst. B
11
(
4
),
1057
(
2009
).
26.
Z.
Xing
,
L.
Wen
, and
W.
Wang
, “
An explicit fourth-order energy-preserving difference scheme for the Riesz space-fractional Sine–Gordon equations
,”
Math. Comput. Simul.
181
,
624
641
(
2021
).
27.
Y.
Zhou
and
Z.
Luo
, “
An optimized Crank–Nicolson finite difference extrapolating model for the fractional-order parabolic-type sine-Gordon equation
,”
Adv. Differ. Equations
2019
,
1
.