On Some Problems of the Numerical Implementation of Nonlinear Systems on Example of KPI Equation

Analytical and numerical peculiarities of solving nonlinear problems are considered on examples of wave equations like KdVB and Kadomtsev-Petviashvili-I equation (KPI). KPI is represented in integro-differential form. Main attention is paid to the problem of asymptotical behavior of solution and appearance of nonphysical artefacts. The numerical solution is carried out by the finite difference method. For a correct representation of the boundary condition along the y axis in numerical simulation a method is proposed for introducing small artificial convection into the original equation in the indicated direction. Along with the introduction of artificial convection, the procedure of trimming of the integral on the bands adjacent to the upper and lower boundaries of the calculated region is used. The results obtained by numerical testing, showed sufficient accuracy and validity of this procedure.


Introduction
The study of nonlinear waves has long been in the focus of researchers [1]. The greatest efforts are associated with the study of waves in solids [2] and in viscous flows [3]. Of particular interest is the study of tsunamis, where the importance of describing two-dimensional wave motions is manifested in its entirety [4]. Although mathematical description of nonlinear waves is well established [5], there are still a lot of questions in realizing numerical equivalent of such a model. In particular, a lot of questions have to be answered with non-integrable numerical analogue of completely integrable system [6]. A good review of the theory and applications of such a waves is given in [7], but there are still a lot of technical questions have to be settled and we try to fill in some of the gaps in our research.
1. Non-physical sources appear. The outcome is the distorted solution and possible loss of stability of the wave.
2. At large times, the solution reaches the fictitious boundary of the computational domain and at the same time: (a) A reflected wave arises; (b) The bogus sources do not go to zero.
So, the correction of the discrete variant of the problem is necessary. In many cases, the infinite region is replaced by a calculated domain with finite dimensions, which makes it necessary to specify the corresponding discrete boundary conditions.

Method to keep artefacts under control
In the exact, analytical formulation of the problem, various forms of the equation representation, both in differential and integral forms are identical, but when passing to a discrete, finite-dimensional space, this equivalence is violated. The method in which the conservation laws are best fulfilled must be chosen.
To overcome problems noted above some steps must be taken in discrete representation of the problem: 1. Sources and sinks that appear in calculation when approaching the fictitious boundary should be zeroed out using a smooth finite function.
2. At the fictitious boundary, boundary conditions are set to exclude reflections.
3. If conditions 1, 2 are not consistent with the transfer part of the equation, then compensating terms should be added to it.

Method demonstration
For testing the proposed method, we carried out test calculations for the two-dimensional Kadomtsev-Petviashvili-I equation (KPI): Equation (1) for the function u(x, y, t) is usually called the equation KPI and is considered in the area with t ≥ 0, x, y ∈ (−∞, ∞), β, η > 0. In general case, the solution of such an equation is sought numerically.
The handling of equation (1) in the finite-difference approach is not effective due to the presence of the mixed derivative u tx .
In the numerical integration of the KPI equation, instead of the original equation (1), its integrodifferential analog is considered, The computational area (plane) of the equation (2) is replaced by a finite computational area [x min , x max ] × [y min , y max ] × [0, T ]. In the computational domain, a uniform difference grid is specified with respect to time t and spatial coordinates x and y: , y min = y 1 , y max = y L , t n = n∆t, n = 0, 1, 2, . . . , T/∆t−1, with ∆x, ∆y -being the spatial coordinate steps, ∆t -being the time step, j 0 -being the index, corresponding to the values of x = 0, u[( j − j 0 )∆x, (k − L/2)∆y, n∆t] = u n j,k . The solution of equation (2) in the half-space t ≥ 0 is sought for the initial distribution u(x, y, 0) = q(x, y).
The Gaussian distribution is considered as the initial distribution: A series of previous papers [8][9][10][11][12][13][14] was devoted to the numerical simulation of non-linear KdVB and KPI type equations. In these papers, a difference scheme was chosen, its testing was carried out, and numerical results were obtained. Various methods for numerical solution of KPI were discussed in [11].
We briefly recall the main features of the difference scheme for the KPI equation [10,11,13]. Derivatives are approximated by the central-difference operators. The difference scheme for equation (2) can be written by approximating the derivatives with the traditional method, for example: with ∆u n+1 j,k = u n+1 j,k − u n j,k and after linearization of (u 2 ) n+1 j,k ) = (u 2 ) n j,k + 2u n j,k ∆u n+1 j,k + O(∆t 2 ). The integral S (x, y, t) is calculated by the trapezoid method, the integrand derivative u yy is approximated by second order central differences: The difference scheme approximating equations (2) can be represented as: The coefficients a j , b j , c j , d j , e j , f j in (4) are determined from the above written relations for the derivatives.
Difference scheme (4) is an implicit linearized finite-difference scheme. The order of approximation of the difference scheme in the computational domain is O(∆t, ∆x 2 , ∆y 2 ). In some necessary cases, a flow correction procedure is used (FCT) [8]. Quasi-one-dimensional system (4) is solved by the five-point sweep method.

The boundary conditions
The correct representation of the differential boundary conditions in finite-difference form is often a difficult task. When the transformation of the original equation (1) to its integro-differential form (2) is done in the continuous domain of the differential equation solution (infinite-number space), then this transformation is unique. Unfortunately, similar transformations with finite-difference schemes in finite-dimensional spaces, with few exceptions, are not identical.
In the transition from the continuous domain of a differential equation to the grid domain of finite difference equations, the infinite region is replaced by a region with finite dimensions. This is quite natural procedure.
Solutions to the KPI equation (1) are decaying proportionally to the square of the distance O[1/(x 2 + y 2 )], as a result, in practical calculations, it is impossible to select a sufficiently large area size for the use of homogeneous Dirichlet boundary conditions. But this immediately raises the question of the definition of appropriate boundary conditions. In such cases, at the boundaries of the calculated area [x min , x max ] × [y min , y max ] traditionally there are used so-called "free flow conditions": u x = u xx = 0 along the boundary lines x min and x max , and u y = 0 along the boundary lines y min and y max . In a differential form, the differential conditions look like: with artificial grid points being introduced in the system (5): We meet here a peculiar feature of the application of difference boundary conditions (5). In fact, these conditions are written for the original Kadomtsev-Petviashvili differential equation (1). When applying it to the corresponding integro-differential analog (2), we do not get the same result for the difference form (5). When calculating the system of difference equations (4) that approximate equation (2) with boundary conditions (5), the flow conditions along the y axis are not correctly realized, but lead to the reflection condition. Fig. 1 shows visible effects of the interaction of propagating waves with the boundary y = y min and y = y max . For this picture and for all following ones, the main numerical values of the parameters are the same. The mesh is 800 × 600 (x × y), ∆x = ∆y = 0.1, ∆t = 5 · 10 −5 , β = η = 1. In our case it is y min = −31.9, y max = 48.

Artificial convection
The original KPI equation (1) in the y direction is of elliptic type. When converting it to the integrodifferential form (5), it becomes the KdV (Korteweg de Vries) equation with the source S (x, y, t) in the right hand side. When calculating such an equation by a finite-difference method in a bounded rectangular region, the integral S introduces a significant error in the solution. Because of the finite region of calculation, the integral S at the boundaries of the region does not tend to zero. In order to get rid of the distortion of the numerical solution, due to the reasons described above, the following procedure is proposed.
1. The sources that appear during the calculation, when approaching the fictitious boundary, are zeroed with the help of a smooth trimming (finite) function.
2. On the fictitious (numerical) boundary, flow conditions are set that exclude reflections.
3. If the conditions 1, 2 are not consistent with the transport part of the equation, then compensating terms should be added.
To be able to correctly set the boundary condition of the flow type u y = 0 at y = y min , y = y max in the left part of equation (2) a small quantity is introduced, which can be called artificial convection over y: −ε · sign(u y ) u y , where the parameter ε > 0 is small enough (ε 1), such as not to introduce a significant error in the solution of the corresponding equation (2). The function sign(u y ) determines the sign of the direction of propagation of the perturbation along the y axis. In the transition to the grid space, it is necessary, if possible, to correctly approximate the function sign(u y ).
In the domain of the equation (2), the integral In our case, when x ≤ x max , this condition is not met, which leads to numerical errors. It is proposed to introduce the trim band δ of the integral on which the integral converges to zero. This trim strip is adjacent to the vertical borders at the bottom. y = y min and top y = y max .
Instead of integral S (x, y, t) in equation (2) we will use the expression: where h − (y), y ∈ [y min , y min + δ] h + (y), y ∈ [y max − δ, y max ] The functions h − (y) and h + (y) fulfil the symmetry condition: h − (y min + y * ) = h + (y max − y * ) at y * ∈ [0, δ]. The choice of the width δ of such band and of the corresponding function h + (y) is discussed below.

Setting artificial convection
Certain difficulties arise when writing a differential expression for artificial convection (6). First, we don't, as such, have a convection velocity. It is simulated by a constant factor ε, the sign of which is associated with the derivative of the function u in the direction of y. Second, the sign sign(u y ) itself in the grid domain depends on the method of approximation of the derivative u y . It is practically possible to implement three methods of approximation: one-sided forward and backward differences and central differences. Verification calculations helped to choose the optimal representation of the artificial convection: , n = 2, 4, 6, . . . sign(u n j,k − u n j,k−1 ), n = 1, 3, 5, . . . If we do not take into account such alternation and/or approximation by one-sided differences, then we can get the deviation of the motion of the soliton peak from the straight line x = 0. The approximation was carried out by central differences, i.e.
In any case, using the relation (7), such a deviation will be small. This result can be explained by the fact that the values of the u y derivatives at the nodes symmetrically located about the axis of symmetry x = 0, with time due to numerical errors begin to differ from each other, which leads to a drift of the soliton.
For sufficiently large ε, errors can occur that are of the nature of the numerical dispersion and lead to the "collapse" of the solution.

Setting the finite function
When choosing the finite function P(y), and practically the function h + (y), there are two criteria: 1. The function h + (y) should monotonously decrease in the band adjacent to the upper boundary y max (h − (y) should increase in the lower band adjacent to y min ).
2. The derivatives of the functions h − (y) and h + (y) at the boundaries y min and y max , respectively, should be as small as possible.
Other conditions and bandwidth δ must be estimated from numerical tests. When choosing the finite function P(y), various functions h + (y), were considered, four functions were selected out of them (the expressions are given for the top band) withȳ = [y − (y max − δ)]/δ,ȳ ∈ [0, 1] at y ∈ [y max − δ, y max ]:

Numerical tests
An analysis of the results got with various finite functions showed that the most suitable for the initial distribution (3) are two functions: the cubic polynomial and the parabola. (See Fig. 2), ε = 0.01. Practically these two functions give the same results.  The width of the smoothing region δ also has a great influence on the results of calculations. It is important to choose the optimal ratio between the function and the smoothing width. Fig. 3 shows the smoothing results for the case of finite function h + 2 (ȳ) = 0.5[1 − tanh(2ȳ)] with three values of δ: 4, 8, and 12% from the total width | y max − y min |.
It is seen that the results for δ = 8% and for δ = 12% are close between each other. It means, a noticeable improvement of the smoothing does not occur when δ extends some certain value. More important role in this situation plays the form of the finite function. It is obvious that with width δ increasing above a certain value, a noticeable improvement of the smoothing does not occur. More important is the form of the finite function.
From the results shown in Fig. 2 and Fig. 3 (middle graph), it is clear that the choice of a parabola or spline for smoothing is better than the choice of a sigmoid function. It is clear that for each specific task a series of test calculations may be required to choose the finite function and the smoothing bandwidth.

Conclusion
The numerical analysis shows, that proposed approach solves at least the problems connected to complete integrability of the initial model. A calculation method is proposed by introducing a smooth finite function and artificial convection. It works quite effectively at least for nonlinear equations with an infinite region of evolution.
In solving of those problems, the features of different types of discretization are revealed.