DeepAI

# Numerical scheme based on the spectral method for calculating nonlinear hyperbolic evolution equations

High-precision numerical scheme for nonlinear hyperbolic evolution equations is proposed based on the spectral method. The detail discretization processes are discussed in case of one-dimensional Klein-Gordon equations. In conclusion, a numerical scheme with the order of total calculation cost O(N log 2N) is proposed. As benchmark results, the relation between the numerical precision and the discretization unit size are demonstrated.

• 1 publication
• 1 publication
10/29/2021

### Quasi-linear analysis of dispersion relation preservation for nonlinear schemes

In numerical simulations of complex flows with discontinuities, it is ne...
05/18/2021

### The BDF2-Maruyama Scheme for Stochastic Evolution Equations with Monotone Drift

We study the numerical approximation of stochastic evolution equations w...
08/24/2020

### Inertia-driven and elastoinertial viscoelastic turbulent channel flow simulated with a hybrid pseudo-spectral/finite-difference numerical scheme

Numerical simulation of viscoelastic flows is challenging because of the...
06/10/2021

### Simulation of viscoelastic Cosserat rods based on the geometrically exact dynamics of special Euclidean strands

We propose a method for the description and simulation of the nonlinear ...
05/25/2020

### Hyperbolic Discretization via Riemann Invariants

We are interested in numerical schemes for the simulation of large scale...
11/26/2019

### Preserving the accuracy of numerical methods discretizing anisotropic elliptic problems

In this paper we study the loss of precision of numerical methods discre...
09/21/2020

### Multidomain spectral approach with Sommerfeld condition for the Maxwell equations

We present a multidomain spectral approach with an exterior compactified...

## 1. Introduction

For concrete examples of hyperbolic evolution equations, one-dimensional linear and nonlinear Klein-Gordon equations are taken. The initial and boundary values problem of one-dimensional Klein-Gordon equations

 (1) ∂2u∂t2+α∂2u∂x2+βF(u)=0, u(x,0)=f(x), ∂u∂t(x,0)=g(x), u(0,t)=u(L,t), ∂u∂t(0,t)=∂u∂t(L,t)

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) ∂u∂t=v, ∂v∂t+α∂2u∂x2+βF(u)=0, u(x,0)=f(x), v(x,0)=g(x), u(0,t)=u(L,t), v(0,t)=v(L,t).

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 Crank-Nicolson method. Consequently the obtained numerical scheme with well-controlled precision is used for the finite-dimensional representation of infinite-dimensional 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:

 u(x,t)=a0(t)+N∑k=1ak(t)cos(2πLkx)+N∑k=1bk(t)sin(2πLkx),v(x,t)=c0(t)+N∑k=1ck(t)cos(2πLkx)+N∑k=1dk(t)sin(2πLkx).

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) da0dt=c0, daldt=cl, dbldt=dl

follow, where . In the same manner, after substituting them to the 2nd equation of (2), multiplying and , and integrating by for , we obtain

 (4) Ldc0dt+β∫L0F(u)dx=0,L2dcldt−α2π2Ll2al+β∫L0F(u)cos(2πLlx)dx=0,L2ddldt−α2π2Ll2bl+β∫L0F(u)sin(2πLlx)dx=0.

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 right-hand side of (4). The operator-conversion 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 two-step treatment substantially reduces the calculation costs arising from the nonlinearity. In fact, using the Crank-Nicolson type numerical scheme, those integrals are approximated by

 (5) ∫L0F(u)cos(2πlxL)dx≃LJJ−1∑j=0cos(2πlxjL),∫L0F(u)sin(2πlxL)dx≃LJJ−1∑j=0sin(2πlxjL),∫L0F(u)dx≃LJJ−1∑j=0F(uj).

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 right-hand 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 left-hand sides and right-hand 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) da0dt=c0,daldt=cl,dbldt=dl,Ldc0dt+βLJJ−1∑j=0F(uj)=0,L2dcldt−2π2l2αalL+LβJJ−1∑j=0F(uj)cos(2πlxjL)=0,L2ddldt−2π2l2αblL+LβJJ−1∑j=0F(uj)sin(2πlxjL)=0,

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

 an+1l−anlΔt=θcn+1l+(1−θ)cnl,

where is a real number satisfying . This is equivalent to

 an+1l=anl+Δt[(1−θ)cnl+θcn+1l],

and similarly

 bn+1l=bnl+Δt[(1−θ)dnl+θdn+1l]

with a notation and . In the second, using notations , and , the fifth equation of (6) become

 (7) L2(cn+1l−cnl)Δt−2απ2l2Lθan+1l−2απ2l2L(1−θ)anl  +βL2θ^Fn+1l+βL2(1−θ)^Fnl=0

which is equivalent to

 (8) cn+1l=cnl+α(2πlL)2Δt[(1−θ)anl+θan+1l]−βΔt[(1−θ)^Fnl+θ^Fn+1l],

and in the same way

 (9) dn+1l=dnl+α(2πlL)2Δt[(1−θ)bnl+θbn+1l]−βΔt[(1−θ)^Gnl+θ^Gn+1l].

Consequently the original problem is fully discretized as

 (10) an+1l=anl+Δt(1−θ)cnl+Δtθcn+1l,\parbn+1l=bnl+Δt(1−θ)dnl+Δtθdn+1l,\parcn+1l=cnl+(1−θ)Δt[α(2πlL)2anl−β^Fnl]+θΔt[α(2πlL)2an+1l−β^Fn+1l],\pardn+1l=dnl+(1−θ)Δt[α(2πlL)2bnl−β^Gnl]+θΔt[α(2πlL)2bn+1l−β^Gn+1l],\parcn+10=cn0−β(1−θ)Δt^Fn0−βθΔt^Fn+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 right-hand 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) an+1,ν+1l=anl+Δt(1−θ)cnl+Δtθcn+1,νl,\parbn+1,ν+1l=bnl+Δt(1−θ)dnl+Δtθdn+1,νl,\parcn+1,ν+1l=cnl+(1−θ)Δt[α(2πlL)2anl−β^Fnl]+θΔt[α(2πlL)2an+1,νl−β^Fn+1,νl],\pardn+1,ν+1l=dnl+(1−θ)Δt[α(2πlL)2bnl−β^Gnl]+θΔt[α(2πlL)2bn+1,νl−β^Gn+1,νl],\parcn+1,ν+10=cn0−β(1−θ)Δt^Fn0−βθΔt^Fn+1,ν0.

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) ^Fn+1,ν0=1JJ−1∑j=0F(un+1,νj),^Fn+1,νl=2JJ−1∑j=0F(un+1,νj)cos(2πLlxj),^Gn+1,νl=2JJ−1∑j=0F(un+1,νj)sin(2πLlxj).

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

 F(u)=u,α=−1,β=1,Ω=[0,L].

The corresponding problem is written by

 (13) ∂v∂t−∂2u∂x2+u=0, ∂u∂t=v, u(x,0)=0,u(0,t)=u(L,t), v(x,0)=cos(2πLx),v(0,t)=v(L,t).

An exact solution of this problem is given by

 (14) u(x,t)=1ωsin(ωt)cos(2πLx),v(x,t)=cos(ωt)cos(2πLx),ω=√1+(2πL)2,

and the corresponding numerical solution is shown in Figs. 1 and 2. In the actual numerical calculations, we further assume

 θ=12,J≥2N+1,L=8.

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

 Error=maxxjerr(xj),

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 th-order polynomial of , the left- and right-hand 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

 F(u)=sinu,α=−1,β=1,Ω=[0,L].

The master equation is known as the Sine-Gordon equation. Using Jacobi’s elliptic functions sn, cn, dn [3], the corresponding problem is written by

 (15) ∂v∂t−∂2u∂x2+sinu=0,∂u∂t=v, v(x,0)=−√2 cn(x,12) dn(x,12)√1−14sn(x,12), v(0,t)=v(L,t), u(x,0)=2sin−1[12sn(x,12)], u(0,t)=u(L,t).

An exact solution of this problem is written by

 (16) u(x,t)=   2sin−1[12 sn(x−√2t,12)], v(x,t)=   −√2 cn(x−√2t,12) dn(x−√2t,12)√1−14 sn2(x−√2t,12).

Using the perfect elliptic integral calculation, can be written by

 (17) L=4F(π2,12)=4∫π/201√1−(12)2sin2θ dθ=6.743001419250385098⋯.

The corresponding numerical solution is shown in Figs. 7 and 8. In the actual numerical calculations, we further assume

 θ=12,J≥2N+1,L=8.

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 9-digit 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 finite-dimensional representation of infinite-dimensional dynamical systems, see Ref. [4]. For the preceding works treating solitons in sub-atomic physics, see Refs.[5-7]. The scheme used in [5-7] is based on the finite-difference 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)”, Springer-Verlag, 2000.
• [3] M. Ohmiya, “Classical analysis of nonlinear waves” (in Japanese), Morikita Publishers, 2008.
• [4] Y. Iwata, and Y. Takei, “Finite-dimensional representation of infinite-dimensional dynamical systems”, to appear in the proceedings of SNA+MC 2020.
• [5] Y. Iwata, “Energy-dependent 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 time-reversal symmetry in many-nucleus systems”, New J. Phys. 21 (2019) 043010.
• [7] Y. Iwata, “Solitons in nuclear time-dependent density functional theory”, Front. Phys. 8:154, 2020.