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.

## I. INTRODUCTION

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,

Here, *u*(**x**, *t*) is an order parameter, $\Omega \u2282Rd(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 domain^{18} have been proposed. Several studies have also been recently conducted on the second-order accurate methods^{4,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}

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

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.

## II. NUMERICAL SOLUTION ALGORITHM

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 Ω = (*L*_{x}, *R*_{x}) × (*L*_{y}, *R*_{y}). Numerical solution algorithms in one- and three-dimensional spaces are described analogously. To discretize Eq. (2), let Ω_{h} = {*x*_{i} = *L*_{x} + (*i* − 0.5)*h*, *y*_{j} = *L*_{y} + (*j* − 0.5)*h*)| 1 ≤ *i* ≤ *N*_{x}, 1 ≤ *j* ≤ *N*_{y}} be the discrete computational domain, where *h* = (*R*_{x} − *L*_{x})/*N*_{x} = (*R*_{y} − *L*_{y})/*N*_{y} is the uniform step size; in addition, *N*_{x} and *N*_{y} are the numbers of the grid points in the *x*- and *y*-directions, respectively. Let $uijn$ be the numerical approximation of *u*(*x*_{i}, *y*_{j}, *n*Δ*t*), where Δ*t* is the time step. To obtain a second-order scheme, we require the following steps:^{11}

where $N\Delta t$ and $L\Delta t$ are the nonlinear and linear solution operators of the PSG equation, respectively. Moreover, the operators $N\Delta tu(x,t)$ and $L\Delta tu(x,t)$ indicate the solutions of *u*_{t} = sin *u* and *u*_{t} = *κ*^{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–Hohenberg^{11} 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:

Before we solve the nonlinear operator $N\Delta t/2$, we consider the following equation:

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

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

By using the Pythagorean trigonometric identity, we obtain

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

which results in

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

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\Delta t/2$(5) as

where $Cijn=0.5[ln(1\u2212cos\u2061uijn)\u2212ln(1+cos\u2061uijn)]$.

Next, to solve the linear operator $L\Delta t$, we consider the diffusion equation

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,\u2026,Nxandj=1,\u2026,Ny}$, the discrete cosine transform is defined as follows:

The inverse discrete cosine transform is

where

*ξ*_{p} = (*p* − 1)/*L*_{x}, and *η*_{q} = (*q* − 1)/*L*_{y}. Let us assume that

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

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

and

where $Cij**=0.5[ln(1\u2212cos\u2061uij**)\u2212ln(1+cos\u2061uij**)]$. 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 *n*th solution satisfies $|uijn|\u2264\pi $. In the first sub-step (5), the solution *u** is bounded by *π* from Eq. (15) for any time step Δ*t*, i.e.,

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),

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

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*,

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.

## III. NUMERICAL EXPERIMENTS

We use the second-order operator splitting scheme (6) to solve the PSG equation. Unless otherwise mentioned, we use numerical parameters as *N*_{x} = 128, *N*_{y} = 128, and *N*_{z} = 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.

### A. Convergence and stability

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(\Vert eNxNt\Vert /\Vert eNx2Nt\Vert )$ 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*π*/*N*_{x} 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.

Case . | Δt
. | Rate . | Δt/2
. | Rate . | Δt/4
. | Rate . | Δt/8
. |
---|---|---|---|---|---|---|---|

l_{2}-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 . | Δt
. | Rate . | Δt/2
. | Rate . | Δt/4
. | Rate . | Δt/8
. |
---|---|---|---|---|---|---|---|

l_{2}-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} |

### B. Boundedness of the numerical solutions

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 *N*_{x} = 2^{8} 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.

### C. Linear stability analysis

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:

For the positive integers *K*_{i} and $x\u2208Rd$ where *d* = 1, 2, 3 is the space dimension, we can obtain

where *α*(*t*) is the amplitude of *u*(**x**, *t*). If *d* = 3, then **x** = (*x*_{1}, *x*_{2}, *x*_{3}). Substituting *u*(**x**, *t*) in Eq. (30) into Eq. (29), we have

By dividing $\u220fi=1d\u2061cos(Kixi)$ on the both sides of Eq. (31), we obtain

We can analytically solve the ordinary differential equation (32) and find the solution as follows: *α*(*t*) = *α*(0)*e*^{λt}, where $\lambda =1\u2212\kappa 2\u2211j=1dKj2$ is the analytic growth rate. We then define the numerical growth rate as

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 $\lambda \u0303$ 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.

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.

### D. Growth simulation

Let us consider the following initial condition:

which implies that the initial values are between 0 and *π*. In this simulation, we consider the numerical parameters as Δ*t* = 0.025, *N*_{x} = 128, *N*_{y} = 128, and *N*_{z} = 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.

## IV. CONCLUSIONS

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}

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts of interest to disclose.

## DATA AVAILABILITY

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