1. Introduction
For concrete examples of hyperbolic evolution equations, onedimensional linear and nonlinear KleinGordon equations are taken. The initial and boundary values problem of onedimensional KleinGordon equations
(1) 
is considered for , where , and are real numbers, and and are initial functions. The inhomogeneous term is either linear or nonlinear function of . This problem is also written by
(2) 
In this manner, the first order evolution problem of hyperbolic type is obtained. This equation is regarded as wave equations. Several conservative quantities in association with the wave propagation is utilized to confirm the precision of scheme.
In this article, for the initial and boundary values problem (2), a precise numerical scheme is proposed. The Fourier spectral method is implemented for the spatial direction, and the scheme is introduced for the time direction. Note that is practically adopted in numerical benchmark tests, which is known as the CrankNicolson method. Consequently the obtained numerical scheme with wellcontrolled precision is used for the finitedimensional representation of infinitedimensional dynamical systems.
2. Discretization
2.1. Discretization of space
The spectral method is employed [1]. The solution of (2) is assumed to be expanded by the Fourier series:
Let us terminate the expansion by the th term. Here, according to the terminology used in the mathematical physics, the space spanned by is referred to the coordinate space, and that spanned by to the momentum space. After substituting them to the 1st equation of (2), multiplying and , and integrating by for ,
(3) 
follow, where . In the same manner, after substituting them to the 2nd equation of (2), multiplying and , and integrating by for , we obtain
(4) 
By solving Eqs. (3) and (4), the values of , , , , , and are calculated.
In terms of dealing with the nonlinearity, we pay attention to the integrals
in the righthand side of (4). The operatorconversion method [1] is employed. At first, by introducing the Fourier inverse transformation, the nonlinear terms are separately calculated in the original coordinate space spanned by . In the second, the nonlinear and linear terms are calculated in the momentum space spanned by . This twostep treatment substantially reduces the calculation costs arising from the nonlinearity. In fact, using the CrankNicolson type numerical scheme, those integrals are approximated by
(5) 
Under the periodic boundary condition, the spatial interval is equally discretized by with . After the spatial discretization, unknown function at is denoted by
. The similarity of the representation between the righthand sides and the Fourier transform simplifies the calculations. Indeed,
is calculated by the discrete Fourier transform (DFT), and , , and are calculated in the next.In the present method, if is th order polynomial, the lefthand sides and righthand sides of (5) coincide for . By utilizing the Fast Fourier transform (FFT), the total calculation cost is reduced to . Consequently the original problem is spatially discretized as
(6) 
where is the discretization in the coordinate space, and is the discretization in the momentum space. Discretizations in two different spatial directions are mixed in this proposed formalism. In this sense, the spatial precision depends on both and . Such a treatment is introduced for maintaining both the precision and the calculation cost. Note that only momentum space discretization is enough for the linear cases.
2.2. Discretization of time
The scheme is employed [2]. For a sufficiently small and natural numbers , we set . At first, using notations and , the time discretization is given by
where is a real number satisfying . This is equivalent to
and similarly
with a notation and . In the second, using notations , and , the fifth equation of (6) become
(7) 
which is equivalent to
(8) 
and in the same way
(9) 
Consequently the original problem is fully discretized as
(10) 
in terms of time and space. Solving (10) means calculating unknown functions , , , , , from the known functions , , , , , . In general, the unknown functions are existing in the righthand side of (9), so that it is necessary to be solved by the iterative scheme.
3. Iterative scheme
The scheme solver is generally an implicit method. That is, the iterative treatment is required. Here we introduce additional parameter for the implicit treatment. Using a new parameter , the intermediate stage is introduced in which , , , , , are separated into , , , , , and , , , , , , and it follows that
(11) 
Using the recurrence relation (11), the th values , , , , , are obtained by the th values , , , , , . This process is sequentially repeated until the convergence. The converged ones satisfy (10) in a numerical sense, then they are renamed as , , , , , . Here and are obtained by
(12) 
In particular being necessary at the intermediate stage is obtained by the DFT in which , , , and are utilized.
The iteration process is the bottle neck of the calculation. Consequently the total calculation cost is mostly dominated by the integration of (10) in the iteration process, whose order is . The order of cost of our total calculation is concluded to be .
4. Benchmark calculations with numerical error tests
4.1. Linear case
For the initial and boundary values problem (2),
The corresponding problem is written by
(13) 
An exact solution of this problem is given by
(14) 
and the corresponding numerical solution is shown in Figs. 1 and 2. In the actual numerical calculations, we further assume
The error comparing the exact and numerical solutions are shown in Figs. 3 and 4, where and cases with , , ,
cases are examined. In order to obtain the error estimates, the difference between the exact and numerical solutions is calculated at each point
.After calculating absolute and relative errors at each point, the minimum error between absolute and relative errors are obtained at each point . By denoting this minimum as , our error function used in the plots is defined by
where, in terms of picking up the error without having a overflow of number representation, the minimum is taken before taking the maximum. Two kinds of errors are used for calculating the relative error and the absolute error. The relative error is usually adopted, and the absolute error is adopted only when the absolute value of is too small for the value of relative error to be handled in the computer. Finally the maximum of error is searched over all the discretized points.
For the error arising from the time discretization, if we apply half of , it results in the quartered error. It simply corresponds to the fact that the present scheme is 2nd order scheme for time direction. That is, it is possible to reduce/control the error rather easily by taking sufficiently small . This information is practical to determine the value of .
Let us compare = case to case. In this comparison, there is no significant difference in error values at least if the same value is chosen for . Note here that the convergence was turned out to be false when with larger was applied. This remarkable property is an advantage of spectral method in which the error does not depend on but on the . This fact is also true in the other plots (Fig. 4). In Fig. 4, the error included in the numerical solution at is shown, and the independence of error is clearly seen. Indeed, in the present numerical scheme, remarkable properties

if is thorder polynomial of , the left and righthand sides of Eq. (4) being used in the calculation of are exactly the same if is satisfied;

for a sufficiently large , the linear solution with a chosen initial value is exactly represented by the Fourier expansion;
are confirmed. These contribute to the independence of possible errors. In conclusion, the error of the present scheme can be controlled by the time discretization.
4.2. Nonlinear case
For initial and boundary values problem (2), we assume
The master equation is known as the SineGordon equation. Using Jacobi’s elliptic functions sn, cn, dn [3], the corresponding problem is written by
(15) 
An exact solution of this problem is written by
(16) 
Using the perfect elliptic integral calculation, can be written by
(17) 
The corresponding numerical solution is shown in Figs. 7 and 8. In the actual numerical calculations, we further assume
which are exactly the same as the linear case. The error comparing the exact and numerical solutions are shown in Figs. 7 and 8, where and cases with cases are examined. In order to obtain the error estimates, the difference between the exact and numerical solutions is calculated at each point . The definition of error function follows from the linear case. For the error arising from the time discretization, if we apply half of , it results in the quartered error. This property is common to both linear and nonlinear problems.
For the error arising from the space discretization, result in Fig. 8 shows nonlinear aspect. Indeed, error decreases depending on at first, while it can be a constant in the next. This tendency is explained by the larger effect of truncation approximation (by ) of Fourier expansion only in smaller cases. Indeed, the error becomes smaller if we change the value of from to . However, if we take sufficiently large , the error values are almost constant. In conclusion, even in a nonlinear case, the error of the present scheme arises mostly from the time discretization if we take sufficiently large .
Let us compare to cases. In this comparison, there is no significant difference in error values at least if the same value for is applied. Note here that the convergence was turned out to be false when with larger was applied. That is, the same statement as the linear case.
5. Summary
A precise numerical scheme for nonlinear hyperbolic evolution equations is proposed; indeed, the value of Error function (corresponding to the relative error) is roughly at the order of for a sufficiently large (Figs. 4 and 8). It ensures the 9digit correctness in those benchmarks.
In terms of providing benchmark results showing the precision, the numerical solutions are compared to the exact solutions for both linear and nonlinear cases. The relation between the numerical precision and the discretization parameters are demonstrated. A high precision has been confirmed to be preserved by taking sufficiently large , while the total calculation cost is only at the order of . The error control parameters such as and
are heuristically found.
Such a precise calculation would be preferably used in the wave propagation and soliton propagation in future studies. For some application results including the finitedimensional representation of infinitedimensional dynamical systems, see Ref. [4]. For the preceding works treating solitons in subatomic physics, see Refs.[57]. The scheme used in [57] is based on the finitedifference methods in which the precision is cared only at the level of obtaining convergence.
Acknowledgement
Numerical calculations have been carried out at supercomputer at YITP, Kyoto University, and workstations at Kansai University and Tokyo Institute of Technology. This work was partially supported by JSPS KAKENHI Grant No. 17K05440.
References

[1]
K. Ishioka, “Introduction to numerical calculations using spectral methods” (in Japanese). Univ. Tokyo Press, 2004.

[2]
R. Dautray and J.L Lions,
“Mathematical Analysis and Numerical Methods for Science and Technology: Volume 5: Evolution Problems I (English Edition)”, SpringerVerlag, 2000.

[3]
M. Ohmiya, “Classical analysis of nonlinear waves” (in Japanese), Morikita Publishers, 2008.

[4]
Y. Iwata, and Y. Takei, “Finitedimensional representation of infinitedimensional dynamical systems”, to appear in the proceedings of SNA+MC 2020.

[5]
Y. Iwata, “Energydependent existence of soliton in the synthesis of chemical elements”, Mod. Phys. Lett. A 30 (2015) 155008.

[6]
Y. Iwata, and P. Stevenson, “Conditional recovery of timereversal symmetry in manynucleus systems”, New J. Phys. 21 (2019) 043010.

[7]
Y. Iwata, “Solitons in nuclear timedependent density functional theory”, Front. Phys. 8:154, 2020.