1. Nonlinear Hyperbolic Evolution Equations

The dynamics of nonlinear hyperbolic equations are fascinating enough to be applicable to wave propagation in any scales from elementary particles to waves in 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. So that to understand fundamental properties of nonlinear waves, it is useful to employ numerical calculations. Note that, while in terms of treating various nonlinear problems (e.g., various type of boundary, discontinuity such as shock propagation) it is practical to specialize numerical schemes individually, a basic and general framework for precisely calculating nonlinear hyperbolic evolution equations had not been well established until Takei-Iwata ^[1]. Here some benchmark calculations of high-precision numerical scheme for calculating nonlinear hyperbolic evolution equations are demonstrated. Much attention is paid to “universal applicability” and “reliability”. Among others, demonstration movies are presented in this article for the first time.

For a concrete problem, we consider one-dimensional nonlinear Klein-Gordon equations (for a textbook, see ^[2]). The initial and boundary values problem of one-dimensional Klein-Gordon equations is written by

(1)

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 or trigonometric functions of u). For constructing the numerical scheme, the equations are written by a system of first-order evolution equations:

(2)

2. Background: Challenges in Numerical Computation

Various numerical schemes have been investigated to reproduce the solutions of nonlinear partial differential equations accurately and efficiently. For nonlinear hyperbolic equations, numerical schemes well reproducing the conservation laws is required to be highly accurate, since the smoothing effect particularly associated with parabolic equations cannot be expected. Here, when typical numerical schemes such as the difference method are applied, we face various problems in the discretization of spatial and time variables as shown below.

2.1. Spatial Discretization

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 ensure the sufficient accuracy. Furthermore, the numerical dispersion inevitably appears in typical finite difference methods in which numerical solutions are known to be difficult to satisfy the conservation laws with certain sufficient accuracy.

2.2. Time Discretization

For the time discretization, explicit methods have been usually employed 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 ^[3][4][5].

2.3. High-Precision Numerical Scheme

For a high-precision numerical scheme, here we present spectral method for spatial discretization and Runge-Kutta method for time discretization ^[1]. By the spectral method, the solution in the wavenumber space is efficient to avoid the numerical dispersion, and to satisfy conservation laws to certain satisfactory degrees. By the two-stage and third-order Runge-Kutta method as unconditionally-stable fully implicit method, the numerical solution with high accuracy can be necessarily obtained stably. The calculation cost of numerical scheme is small enough to be at the order of O(N log₂ N). The precision of the scheme can be found in Figure 1.

Figure 1. Δt (left panel) and Δx (right panel) dependence of error, where Δx^-1is parametrized by N. The proposed scheme corresponds to RK MH u_{err, t13} (for details of notation, see ^[1]).

3. Details of numerical scheme

3.1. Fourier Spectrum Method

The spectral method ^{[6][7][8][9][10][11]} is employed for discretizing the spatial variable. Let the solution of Equation (2) be expanded by the Fourier series

(3)

Then substitute them into the first equation of (2). After multiplying cos(2πkx/L) and sin(2πkx/L) respectively, they are integrated over Ω = [0, L] with respect to x. Similarly, substitute Equation (3) into the second equation of (2). After multiplying cos(2πkx/L) and sin(2πkx/L) respectively, they are integrated over Ω = [0, L] with respect to x. The unknowns are formulated as

(4)

The equations are written by

(5)

In each time step, the solution to the original Equation (2) is obtained by solving Equations. (4) and (5) in which a₀, c₀ and a_l , b_l , c_l , d_l (l = 1, ・ ・ ・ , N) are calculated.

3.2. Implicit Runge-Kutta Method

The implicit Runge-Kutta method ^[12] of two-stage and third-order is employed for discretizing the time variable. Let u(t) be the solution of the initial value problem of the abstract evolution equation

in a Hilbert space. The time interval (α, β) is divided into equally-discretized M segments being incremented by Δt = (β - α)/M. Let the discrete time sequence {t_m} be represented by

and the approximated value of unknown function u(t_m) be denoted by U_m. In this situation, the following method for obtaining the approximate value U_m+1 is called the Runge-Kutta method

Here, the natural number s is called the number of steps, and a_{i j}, b_i, c_i are the parameters that define the formula. It is called the implicit Runge-Kutta method when a_{i j}≠0 (j > i) exists. The conditions

are imposed on the parameters a_{i j}, c_i.

4. Demonstration of high-precision numerical scheme

4.1. Example: Linear Case (Video)

Let us take the initial and boundary value problem (2) with F(u) = u, α = -1, β = 1, Ω = [0, L].

(6)

The initial function is shown in Figure 2, and the calculated movie is shown in Video 1 and Video 2.

Figure 2. Initial function u and v.

4.2. Example: Nonlinear Case (Video)

Let us take the initial and boundary value problem (2) with F(u) = sin u, α = -1, β = 1, Ω = [0, L], L= 6.7430014192503. The nonlinear wave equation in this case is known as the Sine Gordon equation.

(7)

The initial function is shown in Figure 3, and the calculated movie is shown in Video 3 and Video 4.

Figure 3. Initial function u and v.

4.3. Example: Other Cases with High Complexity

In addition, the calculation results of more complicated solutions such as the breather solution for nonlinear hyperbolic equations can be found in Refs. ^[13][14].