Numerical Scheme Based on the Implicit Runge-Kutta Method and Spectral Method for Calculating Nonlinear Hyperbolic Evolution Equations

A numerical scheme for nonlinear hyperbolic evolution equations is made based on the implicit Runge-Kutta method and the Fourier spectral method. The detailed discretization processes are discussed in the case of one-dimensional Klein-Gordon equations. In conclusion, a numerical scheme with third-order accuracy is presented. The order of total calculation cost is O(Nlog2N). As a benchmark, the relations between numerical accuracy and discretization unit size and that between the stability of calculation and discretization unit size are demonstrated for both linear and nonlinear cases.


Introduction
The dynamics of nonlinear hyperbolic equations are fascinating enough to be applicable to wave propagation on any scale, from elementary particles to waves on a cosmic scale. Even fundamental properties have not been fully understood for nonlinear hyperbolic problems, and it is difficult to elucidate these properties only by pure mathematical analysis. To understand the fundamental properties of nonlinear waves we employ numerical calculations. While in terms of treating nonlinear problems (e.g., various types of boundary, discontinuity such as shock propagation) it is important to individually specialize numerical schemes, we are going to establish a basic framework for calculating nonlinear hyperbolic evolution equations. In this context, much attention is paid to "universal applicability" and "reliability". In this paper, as benchmarks for linear and nonlinear hyperbolic problems, we start with reproducing solutions with a simple and general framework (i.e., spatially-continuous solutions under periodic boundary conditions).
We consider nonlinear hyperbolic evolution equations. For concrete examples of hyperbolic evolution equations, here we take one-dimensional linear and nonlinear Klein-Gordon equations (for a textbook, see [1]). The initial and boundary values problem of one-dimensional Klein-Gordon equations is written by u(x, 0) = f (x), u(0, t) = u(L, t), ∂u ∂t (x, 0) = g(x), ∂u ∂t (0, t) = ∂u ∂t (L, t), for (x, t) ∈ [0, L] × [0, T], where α, β and T are real numbers, and f (x) and g(x) are initial functions. The periodic boundary condition is imposed. Inhomogeneous term F(u) is either linear or nonlinear function of u (e.g., the polynomials and trigonometric functions of u). For constructing the numerical scheme, the equations are represented by the first-order evolution equations as for time.
, v(0, t) = v(L, t). (2) Various numerical schemes have been investigated to accurately and efficiently reproduce the nonlinear solutions of partial differential equations such as the nonlinear Klein-Gordon equations mentioned above. For nonlinear hyperbolic equations, numerical schemes well reproducing the conservation laws are required to be highly accurate, since the smoothing effect particularly associated with parabolic partial differential equations cannot be expected.
For the spatial discretization, a conventional finite difference method can be used to discretize the spatial variables in calculating Klein-Gordon equations, but it requires a very small spatial unit ∆x to keep sufficient accuracy. Furthermore, the problem of numerical dispersion is inevitable in typical finite difference methods in which numerical solutions are known to be difficult to satisfy the conservation laws with certain sufficient accuracy. Meanwhile, in the spectral method, the solution in the wavenumber space is known to be efficient to avoid the problem arising from the numerical dispersion, and rather easy to satisfy conservation laws to certain satisfactory degrees. Since the calculation cost of the spectral method is generally higher than that of the finite difference method, some approximation methods have been used to improve the feasibility of calculations (cf. pseudo-spectral method [2][3][4][5], collocation method [6,7]). Note that the spectral method is a discretization method in which the solution is represented by a linear combination of a finite number of Fourier series, and therefore the boundary condition is a periodic boundary condition. On the other hand, numerical schemes (e.g., spectral element method, spectral penalty method, etc. [8][9][10][11]) have been also developed for various boundary conditions such as Dirichlet boundary condition, Neumann boundary condition, and other boundary conditions.
For the time discretization, explicit methods have been used in calculating Klein-Gordon equations, but it requires additional treatments to obtain sufficiently accurate solutions. While linear and nonlinear solvers based on explicit methods are relatively simple with low computation cost, it is also known that there is a restriction called the Courant-Friedrichs-Lewy Condition (CFL condition, for short) on the time unit ∆t in order to obtain numerical results stably. Therefore, a numerical scheme combining explicit and implicit methods has been studied: e.g., a method for weakening the restriction on ∆t with keeping the stability of calculations [12][13][14]. In this context unconditionally-stable fully implicit method has been studied recently [15][16][17]. When it comes to a fully implicit method, the problem is generally reduced to a nonlinear type self-consistent equation, and it is necessary to apply numerical iteration. Numerical iterative methods generally require extra-ordinary high computational costs compared to the explicit method (non-iterative methods). The convergence and efficient implementation of numerical iterative methods have also come into the recent spotlight.
In this paper, we discuss how to construct a high-precision scheme for hyperbolic evolution equations. Based on the two-stage and third-order implicit Runge-Kutta method [18] and the spectral method [2,4,[19][20][21], the order of computation is confirmed to be O(N log 2 N) in the benchmark calculations. This order estimate simply shows the feasibility and the applicability of the proposed scheme. Consequently, the precision and stability of the numerical scheme are examined for linear and nonlinear cases. As a remark for the limitation of the proposed method, we will not go into a detailed discussion of boundary conditions other than the periodic boundary conditions. However, by incorporating spectral elements methods, spectral penalty methods, and so on, it is definitely possible to construct numerical schemes for the other boundary conditions.

Discretization of Space Using Spectral Method
The spectral method [2,4,[19][20][21] is employed to discretize the spatial variables. Let the solution of Equation (2) be expanded by the Fourier series.
Then substitute them into the first equation of (2). After multiplying cos 2π L lx and sin 2π L lx respectively, they are integrated for Ω = [0, L] with respect to x. It follows that Similarly, substitute Equation (3) into the second equation of (2). After multiplying cos 2π L lx and sin 2π L lx respectively, they are integrate for Ω = [0, L] with respect to x. It follows that The solution to the original Equation (2) is obtained by solving Equations (4) and (5) in which a 0 , c 0 and a l , b l , c l , d l (l = 1, · · · , N) are found.
In terms of dealing with the nonlinearity, the following integral values appearing in Equation (5) is the bottle neck of the computational cost.
To deal with these terms, we use the operator transformation method [2,4,19,20]. Its main idea is not to solve the nonlinear term in the Fourier-transformed momentum space, but in the original coordinate space by carrying out the inverse Fourier transform. This procedure remarkably reduce the computational cost; indeed, if we introduce the approximation based on the trapezoidal formula, the nonlinear terms are written by where, under the periodic boundary condition, the equidistant J division of Ω is denoted as x j for j = 0, · · · , J, and the value of u in time t at each point is denoted as u j = u(x j , t).
Here, the right-hand side of the above equation is characterized by the fact that it is expressed in the same form as the discrete Fourier transform. Therefore, when calculating the integral (6) numerically, u j is obtained by the discrete inverse Fourier transform from a 0 (t), a l (t), and b l (t), (l = 1, · · · , N), and the right side of Equation (6) is obtained by the discrete Fourier transform from F(u j ), (j = 0, · · · , J − 1). Note that if F(u) is a M-degree polynomial of u, the values of the left and right sides coincide for J ≥ (M + 1)N + 1 [2,4,19,20]. Furthermore, by using the Fast Fourier Transform (FFT) for the above discrete inverse Fourier transform and discrete Fourier transform, the computational cost of (6) becomes O(N log 2 N). Consequently, the spatially-discretized equation becomes based on the spectral method with the operator transformation treatment. For the discretization of spatial variables by the spectral method, the pseudo-spectral method and collocation methods are sometimes applied as effective numerical solution methods [2][3][4][5][6][7].
In particular, integral values are approximately obtained by weakening the condition in which the number of fractional points is taken to be J ≥ (M + 1)N + 1. In such situations, there is a risk that the conservation law may not be well reproduced due to an aliasing error arising from the overlaps between different wave number components. This error is caused by the superposition of high wave number components. If the high wave number components of the solution are sufficiently small compared to N by keeping the cut-off wave number N to be sufficiently large, it cannot be a significant problem. Therefore, pseudo-spectral method and collocation methods are sometimes preferred.

Matrix Form
As a preparation for the discretization of the time variables, Equation (7) is represented as a matrix form. Vectors a, b, c, d are defined as follows.
Here, using Equations (8) to (10), the matrix representation of Equation (7) is as shown below.
where 0 = (0, 0, · · · , 0) t . Using the notations Equation (11) is represented by where g, h depend on {u i } J−1 i=0 , which is obtained by the inverse Fourier transform of W. Therefore, after the spatial discretization, the inhomogeneous term generally holds the nonlinearity F = F(W).

Implicit Runge-Kutta Method
Following the literature [18], we introduce the implicit Runge-Kutta method of twostage and third-order for discretizing the time. Let u(t) be the solution of the initial value problem of the abstract evolution equation.
in a Hilbert space (corresponding to Equation (13) and the approximated value of unknown function u(t m ) is denoted by U m . In this case, the method for obtaining the approximate value U m+1 is called the Runge-Kutta method. More precisely, the Runge-Kutta method is represented by Here, the natural number s is called the number of steps, and a ij , b i , c i are the parameters that define the formula. It is called the implicit Runge-Kutta method when a ij = 0 (j > i).
The conditions is known as the Butcher tableau [18]. Under the assumption that U m = u(t m ), the local discretization error is defined by Note that when the local discretization error is evaluates as T n+1 = O((∆t) p ), the Runge-Kutta formula is said to be of order p. Although the number of stages and orders is arbitrary, in this paper, we adopt the Runge-Kutta method with two-stage and third-order by taking into account the balance of the accuracy and the calculation cost. The Butcher tableau for the implicit Runge-Kutta formula of two-stage and third-order is given by In general, explicit schemes have a restriction on the setting of the time spacing variables ∆t, called the CFL condition. Therefore, it is important for numerical schemes to be A-stable.
On the other hand, the Implicit Runge-Kutta Method is known to be A-stable. In this sense, the Implicit Runge-Kutta Method is preferably employed as a stable and high-order scheme [22][23][24]. However, the Implicit Runge-Kutta Method is not widely used. This is firstly due to the calculation cost; indeed the application of implicit schemes to nonlinear partial differential equations requires numerical iteration for each single time step. Another issue is that the convergence of the iterative method is not guaranteed when the degree of nonlinearity and the time increment range of the target problem are large. The numerical scheme proposed in this paper achieves a good balance between convergence (stability) and simplicity by using simple devices without applying complicated methods (cf. Section 3.4).

Application of Implicit Runge-Kutta Method and Iteration Formula
Applying the two-stage and third-order implicit Runge-Kutta method to discretize the time variables in (13), the resulting equations are shown by where, using the (N , vectors k 1 and k 2 are represented by Time evolution of solution can be found by calculating k 1 and k 2 , and substituting them into the first equation of (14). The second equation of (14) can be decomposed into four parts.
Similarly, the third equation of (14) is decomposed into four parts.
k 1 , k 2 satisfying Equations (15) and (16) are found by means of an iterative method. Let the value of the ν-th iteration be represented by Then the formula for finding the (ν + 1)-th value from the ν-th, which corresponds to the transforms (15) and (16) are summarized as This is a fully discretized equation for both time and space.

Implementation of Half Step
To solve Equation (17), take k 1 1 = W n , k 1 2 = W n as the initial value. Iteration is carried out until both k ν 1 and k ν 2 converge. In order to make the iteration process as short as possible (i.e., k ν 1 and k ν 2 preferably converge in smaller iteration numbers). we introduce k ν+1/2 1 and k ν+1/2 Accordingly Equation (17) is modified to In the following, we describe the procedure for calculating Equation (18) in detail. The calculation consists of two stages.

The First Stage
Let the l-th element of the vector k X,ν x (X = a, b, c, d; x = 1, 2) be made of k X,ν x,l (l = 0, 1, · · · , N), and let the l-th element of the vector a, b, c, d be made of a l , b l , c l , d l (l = 0, 1, · · · , N) respectively. Calculations are performed for each components.

Decision of Convergence
If is a sufficiently small positive value (e.g., 10 −12 ), and if the smaller value of the relative or absolute error between the ν-th and (ν − 1)-th iterations is smaller than , then the ν-th iteration is regarded as converged. If k ν+1 1 and k ν+1 2 are not converged, return to the first stage of calculation. If k ν+1 1 , k ν+1 2 is convergent, then W n+1 is obtained using the formula The calculation cost of this numerical scheme applying the implicit Runge-Kutta method and the spectral method is estimated to be O(N log 2 N), and the local discretization error of the time variable is estimated to be of third-order.

Comparison to Exact Solution
Let us take the initial and boundary value problem (2) with F(u) = u, α = −1, β = 1, The exact solution to this problem is given by For Equation (24), the time evolutions of the numerical solution u, v (J ≥ 2N + 1, L = 8) until t = 1 are shown in Figure 1, and those at t = 1 are shown in Figure 2. Since the exact solution is represented by Equation (25), it is possible to evaluate the accuracy depending on the spatio-temporal discretization. Furthermore, in terms of clarifying the merit of the proposed scheme, the error between the numerical results by the θ method with θ = 1/2 (for the preceding calculation using the θ method in time, see [21]) are also calculated. The comparison between the numerical and the exact solutions are shown in Figures 3 and 4. In the present paper, the error is defined as the smaller one of the relative and absolute errors. In Figure 3, numerical solutions with different time spacing unit ∆t are examined. The maximum error, which is determined by running for x j (j = 0, · · · , J), between the numerical and the exact solutions are calculated. We see that exponential dependence on the time discretization is noticed. It also shows that the high accuracy of implicit Runge-Kutta method; indeed at ∆t < 10 −4 , almost 10 −4 times accurate result can be obtained compared to the θ-method with θ = 1/2 (the Crank-Nicolson method). In Figure 4, numerical solutions with different spatial spacing parameter N are examined. The maximum error, which is determined by running for x j (j = 0, · · · , J), between the numerical and the exact solutions are calculated. We see that in the logarithmic scale, there is no significant N dependence in calculations of any kinds. Figure 3 shows ∆t-dependence of errors. As for the two-stage and third-order implicit Runge-Kutta method (+ and • in the figure), we see that the error becomes 1/8 times smaller if ∆t is taken to be 1/2 times smaller. In other words, the numerical results confirm that the scheme holds the third-order accuracy with respect to time. On the other hand, as for the θ method (× and in the figure), we see that the error becomes 1/4 times smaller if ∆t is taken to be 1/2 times smaller.

Accuracy Depending on Discretization of Time Variables
For each of the above two numerical schemes, by comparing the errors in case of N = 2 5 and N = 2 10 , the values of the errors are almost unchanged if the amplitude of ∆t is the same (note that in case of N = 2 10 , the numerical calculation does not converge and no numerical solution is obtained for large ∆t). That is, if N is taken to be sufficiently large, the error depends only on the size of ∆t and not on the size of N. If we limit ourselves to the linear case of our benchmark calculations, this advantage arises essentially from the introduction of the Fourier spectral method for the spatial direction.
In any case, the error depends only on ∆t for sufficiently large N, and the smaller ∆t results in the smaller error. Comparing two different schemes, it is concluded in the linear case of benchmark that the implicit Runge-Kutta method possibly includes 10 −4 times smaller errors compared to the θ method.     Figure 4 shows that the error does not depend on N = 2 2 , 2 3 , · · · , 2 15 (refer to + and • in Figure 4 for the errors by implicit Runge-Kutta method, and to × and for those by θ method). Regardless of the choice of N, it is confirmed in Figure 4 that the errors with ∆t = 2 −13 are smaller than those with ∆t = 2 −6 (note that in case of ∆t = 2 −6 , the numerical calculation does not converge and no numerical solution is obtained for large N). This means that the error depends only on ∆t and not on N.

Accuracy Depending on Discretization of Spatial Variables
In conclusion, most of the error arises from the discretization of time. Therefore, if we limit ourselves to the benchmark (the linear case), for example, by taking ∆t ∼ 10 −4 and N ∼ 2 5 , we can obtain high-precision numerical solutions.

Convergence/Stability of Iteration
The numerical scheme proposed in this paper is A-stable, because it employs the implicit Runge-Kutta method. However, the problem to be solved by the implicit scheme results in numerical iterations. Since the convergence of the iterative method is not guaranteed in general, the convergence of the iterative method and the speed of convergence, which indicate a kind of stability in the case of applying implicit Runge-Kutta methods, must be discussed. Here, for improving convergence and accelerating an iterative process, the intermediate step shown in Equation (18) is equipped. This iterative method is referred to as the modified iterative method in this paper. In order to examine the efficiency of the modified iterative method, we compare the average number of iterations in calculating the time evolution of Equation (24) up to t = 1 with several N and ∆t. In addition, the results of applying the iterative method with Equation (17) (normal iterative method) are also shown for comparison.
The exact solution to this problem is given by where sn, cn, dn are Jacobi's elliptic functions, and L can be expressed using complete elliptic integral of the first kind (See Appendix A and textbook [25]).
For Equation (27), the time evolutions of the solution u, v (with J ≥ 2N + 1, L = 6.7430014192503) until t = 1 are shown in Figure 5, and those at time t = 1 are shown in Figure 6. Since the exact solution is represented by Equation (28), it is possible to evaluate the accuracy depending on the spatio-temporal discretization. Furthermore, similar to the linear case, the error between the numerical results by the θ method with θ = 1/2 are also calculated. The comparison between the numerical and the exact solutions are shown in Figures 7 and 8. In Figure 7, numerical solutions with different time spacing unit ∆t are examined. The maximum error, which is determined by running for x j (j = 0, · · · , J), between the numerical and the exact solutions are calculated. The exponential dependence on the time discretization is also noticed in the nonlinear case. It also shows that the high accuracy of implicit Runge-Kutta method; indeed at ∆t < 10 −4 , almost 10 −4 times accurate result can be obtained compared to the θ-method with θ = 1/2. In Figure 8, numerical solutions with different spatial spacing parameters N are examined. The maximum error, which is determined by running for x j (j = 0, · · · , J), between the numerical and the exact solutions are calculated. We see that, in the logarithmic scale, there is no significant N dependence when N is sufficiently large.     Figure 7 shows ∆t-dependence of errors. As for the two-stage and third-order implicit Runge-Kutta method (+ and • in the figure), we see that the error becomes 1/8 times smaller if ∆t is taken to be 1/2 times smaller. In other words, even in the nonlinear case, the numerical results confirm that the scheme holds the third-order accuracy with respect to time. Also in nonlinear cases, as for the θ method (× and in the figure), we see that the error becomes 1/4 times smaller if ∆t is taken to be 1/2 times smaller.

Accuracy Depending on Discretization of Time Variables
For each of the above two numerical schemes, by comparing the errors in case of N = 2 5 and N = 2 10 , the values of the errors are almost unchanged if the amplitude of ∆t is the same (note that in case of N = 2 10 , the numerical calculation does not converge and no numerical solution is obtained for large ∆t). That is, if N is taken to be sufficiently large, the error depends only on the size of ∆t and not on the size of N. If we limit ourselves to our benchmark calculations(also in the nonlinear case), this advantage arises essentially from the introduction of the Fourier spectral method for the spatial direction.
In any case, the error depends only on ∆t for sufficiently large N, and the smaller ∆t results in the smaller error. Comparing two different schemes, it is concluded also in the nonlinear case of benchmark that the implicit Runge-Kutta method possibly includes almost 10 −4 times smaller errors compared to the θ method. Figure 8 shows that if N is larger than 2 5 , the error does not depend on the size of N for either ∆t = 2 −6 , 2 −13 (refer to +, • in Figure 8 for the errors by implicit Runge-Kutta method, and to ×, for those by θ method). In addition, the smaller ∆t leads to a smaller error. On the other hand, if N is smaller than 2 3 , the smaller N leads to the larger error in which the order of error is almost the same regardless of the choice of the scheme and the amplitude of ∆t. That is, most of the error arises from the discretization of space variables if N is smaller than 2 3 , while most of the error arises from the discretization of time variables if N is larger than 2 5 .

Accuracy Depending on Discretization of Spatial Variables
In conclusion, if N is sufficiently large, the error associated with the discretization of spatial variables becomes sufficiently small, and most of the error arises from the discretization of time. Therefore, similar to the linear case, the benchmark calculations suggest that high-precision numerical solutions can be obtained by taking ∆t ∼ 10 −4 and N ∼ 2 5 .

Convergence/Stability of Iteration
As in the linear case, in order to confirm the properties of the iterative method from the viewpoint of stability in a broad sense, we present the results of the average number of iterations in calculating the time evolution of Equation (27) up to time t = 1 with several N and ∆t. In addition, the results of applying the iterative method with Equation (17) (normal iterative method) are also shown for comparison.

Summary
In order to understand the dynamics of nonlinear hyperbolic equations, the numerical approach is one of the most efficient tools, where the smoothing effect is hopeless for generic hyperbolic equations. In this sense, accuracy is an indispensable factor in any way (for the conservation of physical quantities of the present scheme, see [26,27]). A high precision numerical scheme for nonlinear hyperbolic evolution equations is proposed, and its performance is examined by linear and nonlinear benchmark calculations. The numerical scheme consists of the implicit Runge-Kutta method and the Fourier spectral method. For a concrete example, the demonstration of scheme is carried out for one-dimensional Klein-Gordon equations. The precision due to the implicit Runge-Kutta method and spectral method is quantitatively shown by the errors of benchmark calculations in comparison to the θ method, and in comparison to the difference of time spacing variables.
As demonstrated in benchmark cases, it is confirmed that the scheme constructed in this paper has a third-order accuracy with respect to time. We have also confirmed that the truncation error associated with the spatial discretization is much smaller than the error associated with the time discretization by setting N to be sufficiently large. That is, for sufficiently large N, most of the error arises from the error associated with ∆t.
In time discretization, due to the limitation of the time spacing variables by the CFL condition, it is generally difficult to achieve both small error and low computational cost simply by the conventional explicit method. Furthermore, in spatial discretization, due to the effects of numerical dispersion, it is generally difficult for the typical finite difference method to preserve the conserved quantities to a sufficient degree. The proposed method settles these two points by simultaneously introducing the implicit Runge-Kutta method and the spectral method. The applicability of the scheme is confirmed by its computational order O(N log 2 N). In conclusion, the proposed scheme provides a generic framework for high-precision computation with relatively low computational cost.
From a future perspective, more complex problems with various boundary conditions and the behavior of solutions with discontinuities, construction of numerical scheme employing the spectral element method and/or the spectral penalty method is a promising direction. The proposed scheme should be a steadfast basis for such future works.

Conflicts of Interest:
The authors declare that there is no conflict of interest.

Appendix A. Jacobian Elliptic Functions
A set of basic elliptic functions was introduced by Carl Gustav Jacob Jacobi [28] in 1829. These functions are named the Jacobian elliptic functions after him. For any k ∈ [0, 1), we define K(k) by the incomplete elliptic integral of the first kind.
Then, for any k ∈ [0, 1) and any x ∈ [−K(k), K(k)], we define sn(x, k) by an inverse of the incomplete elliptic integral of the first kind. .