DEVELOPMENT OF AN ALGORITHM FOR CALCULATING STABLE SOLUTIONS OF THE SAINT- VENANT EQUATION USING AN UPWIND IMPLICIT DIFFERENCE SCHEME

R a k h m a t i l l o A l o e v Doctor of Math-Physics Sciences, Professor* A b d u m a u v l e n B e r d y s h e v Corresponding author Doctor of Math-Physics Sciences, Professor, Head of Department** Сhief Researcher*** E-mail: berdyshev@mail.ru A z i z a A k b a r o v a Postgraduate Student* Z h a r a s b e k B a i s h e m i r o v PhD, Associate Professor** Сhief Researcher*** *Department Computational Mathematics and Information Systems National University of Uzbekistan Universitet str., 4, Tashkent, Republic of Uzbekistan, 100174 **Department of Mathematics and Mathematical Modeling Institute of Mathematics, Physics and Informatics Abai Kazakh National Pedagogical University Dostyk аve., 13, Almaty, Republic of Kazakhstan, 050010 ***Institute of Information and Computational Technologies Pushkin str., 125, Almaty, Republic of Kazakhstan, 050010 The problem of numerical determination of Lyapunov-stable (exponential stability) solutions of the Saint-Venant equations system has remained open until now. The authors of this paper previously proposed an implicit upwind difference splitting scheme, but its practical applicability was not indicated there. In this paper, the problem is solved successfully, namely, an algorithm for calculating Lyapunov-stable solutions of the Saint-Venant equations system is developed and implemented using an upwind implicit difference splitting scheme on the example of the Big Almaty Canal (hereinafter BAC). As a result of the proposed algorithm application, it was established that: 1) we were able to perform a computational calculation of the numerical determination problem of the water level and velocity on a part of the BAC (10,000 meters) located in the Almaty region; 2) the numerical values of the water level height and horizontal velocity are consistent with the actual measurements of the parameters of the water flow in the BAC; 3) the proposed computational algorithm is stable; 4) the numerical stationary solution of the system of Saint-Venant equations on the example of the BAC is Lyapunov-stable (exponentially stable); 5) the obtained results (according to the BAC) show the efficiency of the developed algorithm based on an implicit upwind difference scheme according to the calculated time. Since we managed to increase the values of the difference grid time step up to 0.8 for calculating the numerical solution according to the proposed implicit scheme


Introduction
It is known that the behavior of water in rivers, lakes, oceans, as well as in small reservoirs is described by the Saint-Venant equations (also called the shallow water equation in the literature). The Saint-Venant equations are a system of hyperbolic partial differential equations that describe flows under the surface of a liquid.
The one-dimensional theory of non-stationary flows without discontinuities, i.e., flows that are not accompanied by the formation of discontinuous waves, can be considered the most developed today (including the problem of natural flooding).
The shallow water approximation (shallow water equations) has great theoretical and practical importance. Characteristic of the processes under study, an analytical solution is almost impossible, therefore, solutions to practical problems are found using numerical methods. The application of the Saint-Venant equation to the solution of various practical problems has a rich history. The variety of problems associated with the study of fluid flows with a free surface has given rise to a significant number of numerical methods that take into account the features of a particular problem under consideration. Among them, we can distinguish such methods as the finite difference method, the particle method, the finite element method, the finite volume method and others, more specific, based, for example, on the splitting of steps, the Godunov method, etc. The approaches used cannot be considered perfect. The main difficulties in constructing a numerical approximation of the Saint-Venant equation are associated with obtaining a stable difference solution to the problem. The question of the stability of the Lyapunov numerical solution (exponential stability) of difference schemes is particularly acute. This question is practically not investigated for hyperbolic systems. However, this issue is practically related to the management and monitoring of water resources, especially for open canals. Other difficulties arise when approximating the terms of the lower terms. For their approximation, for example, the methods of splitting by the lowest terms are used. All this leads to the complication of difference algorithms and often deprives them of uniformity. Thus, the task of improving and developing effective algorithms for mathematical modeling of flows in the Saint-Venant approximation is urgent.
In recent years, the development of numerical models of the dynamics of surface waters in the shallow water approximation has been actively developed.
Great progress has been achieved in studying the general properties and patterns of shallow water modeling in systems with open canals and reservoirs. However, it is important to study specific hydrological objects. Modeling of river systems and reservoirs seems to be the most difficult. As positive examples of constructing multidimensional models for such objects, we point to the lower course of the Bureya River behind the Bureyskaya HPS [1], the Danube Delta [2], the Kuibyshev reservoir [3], the Don estuary area [4], the middle course of the river Don [5], the Medveditsa River [6].
Thus, the task of developing algorithms for mathematical modeling of flows in the Saint-Venant approximation is urgent.

Literature review and problem statement
Currently, various methods for solving the Saint-Venant equation have been developed, such as the method of characteristics, the method of grids, the Lyapunov quadratic function method, the variational method, etc. In [7], for solving one-dimensional equations of unsteady motion in open canals, an implicit difference grid scheme was developed that allows calculations to be performed with a large time step. It was especially important for calculating floods in large rivers when the calculation process took quite a long time. However, there is no description of the algorithm and justification of the finite difference method itself. [8] is devoted to the exponential stability of non-linear Saint-Venant equations in differential form. An explicit quadratic Lyapunov function is constructed, and the local exponential stability is proved.
In [9], the theory of a symmetric hyperbolic system is described. Particularly, in the case of two spatial and one temporal variable, the existence theorem of a dissipative mixed problem is stated. Numerical calculations of simple models are given as examples.
The monograph [10] is devoted to the study of mixed problems for one-dimensional hyperbolic systems in canonical form. Lyapunov stability is established in various functional spaces, in particular, many practical models are considered. The problems of numerical solutions of mixed problems are not considered.
In [11], for the telegraph equation, the discrete Lyapunov function is constructed and its decreasing is proved. Using such an approach to the Saint-Venant equation is due to the difficulties that require additional research.
For two-dimensional hyperbolic equations with dissipative boundary conditions, the exponential stability of the solution is established by the Lyapunov function method in [12].
In [13], algebraic conditions for exponential stability of the solution of mixed problems of linear Saint-Venant equations are obtained. The question of the numerical solution is not considered.
For one-dimensional quasi-linear hyperbolic systems, problems with dissipative boundary conditions that guarantee exponential stability of classical solutions are considered in [14]. There are no studies on numerical calculations.
In [15], using the second kind Volterra transformation and reversible Fredholm transformation, optimal management problems for general linear hyperbolic systems are investigated.
Note that in the works [8,10,[12][13][14][15] we study the issues related to the theoretical aspects of the solvability and stability of mixed problems for hyperbolic systems and the issues of constructing numerical solutions and the stability of difference schemes are not considered. During the numerical calculation of mixed problems for hyperbolic systems, the reason for the above is that the dimension of the linear algebraic equations system increases with an increasing dimension of the considered area. This leads to an unreasonably large amount of computing and requires the involvement of high-performance computing equipment.
In [16], the use of various difference grids (rectilinear and curved) for the numerical calculation of linear partial differential equations is shown. The paper does not consider the problems of the Saint-Venant equation, this is apparently due to the nonlinearity of the equation.
In [17], the discretization of equations describing the unsteady flow of a viscous incompressible fluid is considered based on the finite difference method and the splitting scheme by physical factors on a rectangular non-uniform grid with a staggered arrangement of nodes. However, there is no justification for the convergence and stability of the difference scheme.
In [18], the authors propose a class of difference schemes for hyperbolic systems of equations that have several forms of notation. The stability of the proposed difference schemes is investigated using the technique of energy integrals. However, their application to the study of exponential stability is a rather difficult task.
The work [19] is devoted to the study of initial boundary value problems for a class of three-dimensional quasi-linear hyperbolic systems. An a priori estimate of the problem solved using the energy integral method is obtained. The problems of the numerical solution and its stability have not been studied.
In [20], a linear initial boundary value problem of the dynamics of fluid-saturated porous media, described by three elastic parameters in a reversible hydrodynamic approximation, is solved numerically. The issue of computational model adequacy remains open.
In [21], a problem for the Saint-Venant-Exner equation (SVE), which describes the dynamics of water in a canal filled with sediments with arbitrary values of the canal bottom slope, friction, porosity, as well as the interaction of water and sediment under subcritical or supercritical flow regime is considered. However, this is the subject of further research. We consider the case of a canal without sediment.
The book [22] is devoted to the solution methods of high-order algebraic systems, which appear when using the grid method to the problems of mathematical physics. Along with iterative methods, which are most widely used in computational practice in solving these problems, direct methods are also described. Here we use the matrix sweep method about the Saint-Venant equations.
The development of numerical models of surface water dynamics in the shallow water approximation has been actively developed in recent years by the efforts of researchers. Great progress has been achieved in studying the general properties and patterns of shallow water modeling in systems with open canals and reservoirs. However, it is important to study specific hydrological objects. Modeling of river systems and reservoirs seems to be the most difficult. As positive examples of constructing multidimensional models for such objects, we point to the lower course of the Bureya River behind the Bureyskaya HPS [1], the Danube Delta [2], the Kuibyshev reservoir [3], the Don estuary area [4], the middle course of the river Don [5], the Medviditsa River [6].
In [1], the results of numerical modeling of the spillway of the Bureyskaya HPS at various discoveries spillway spans are given.
A modified shallow water model was used as a mathematical model. The results of numerical modeling and experimental studies of currents at the spillway dam and in the downstream of the Bureyskaya HPS are in good agreement with the experimental data. However, the work does not provide numerical results.
On the basis of numerical modeling in [2], the influence of the Danube river runoff on the formation of the hydrological structure of waters and the peculiarities of circulation on the northwestern shelf of the Black Sea is investigated. This research uses a three-dimensional sigma-coordinate numerical model of the joint dynamics of a shallow sea. No numerical results are available.
In [3], the spatial variability of bottom sediments and the heterogeneity of the distribution of phosphorus in the water area of the near-dam reach of the Kuibyshev reservoir are considered. The results show the calculations of sedimentation and bottom erosion during the rise and fall of water.
In [4], HEC-RAS was used in hydrological processes modeling in the Don's delta, as well as digital terrain model (DTM), which is created by combining Japan DTM, topographic maps and channel data of the united navigable deep-water system of the European part of Russia. The authors carried out calculations of the flow distribution in the delta channels, periods of water changes in the canals and modeling of floods using a wind surge.
In [5], on the basis of computational experiments in a two-dimensional formulation, the features of the change in the dynamics of the flows of p are considered. Don when implementing various bank protection options. The initial materials for constructing the model were the results of a detailed bathymetric survey of this section of the watercourse, the materials of the geodetic survey of the coastal strip, the results of hydrometric measurements, which included the measurement of flow rates and the determination of the granulometric composition of bottom sediments. To ensure a correct assessment of the distribution of the velocity field and the degree of bank destruction in the considered section of the river Don (500 m above the Basovsky arm channel -500 m below the channel), calculations were performed in a two-dimensional formulation. A two-dimensional (in the horizontal plane) model for the considered section of the river Don is made using a specialized hydrodynamic package SMS v.10.1 of the American company AQUAVEO LLC, developed by order and with the participation of the US Center for Hydraulic Research. The water surface levels and horizontal components of the velocity field are calculated. Numerical results are not shown.
The possibility of solving various kinds of problems related to the dynamics of surface waters in different territories in non-stationary conditions is the use of numerical modeling methods together with parallel technologies on various supercomputers. In [6], the results of modeling the site of clearing the Medveditsa river channel are presented. Numerical integration of the equations of hydrodynamics, in the Saint-Venant approximation, taking into account a wide range of physical factors: realistic terrain, operation of hydraulic structures, evaporation, infiltration, bottom friction, wind regime, rotation of the Earth, etc., forms the basis of a computer model of the dynamics of surface water. Numerical results are not presented in the paper.
This review allows the claim that it is advisable to conduct a study on the numerical calculation of exponentially stable solutions of the Saint-Venant equation.

The aim and objectives of the study
The aim of this study is to develop an effective algorithm for finding a numerical solution using an implicit upwind difference scheme for a linear Saint-Venant equations system in the general case. This allows for a numerical calculation to determine the level and speed of water on the BAC.
To achieve this aim, the following objectives were set: -to formulate an implicit upwind difference scheme; -to apply the matrix sweep method to the Saint-Venant equations to determine its numerical solution; -to conduct a computational experiment for the problem of numerical determination of the water level and speed on the BAC section (10.000 m) located in the Almaty region, the Republic of Kazakhstan.

Materials and methods
According to [8], we consider the following hyperbolic system with variable coefficients and with minor terms: and with boundary conditions at x = 0, L respectively: and with initial data at t = 0: Here: H * (x), V * (x) are known stationary solutions of the Saint-Venant equations system; H(t, x) is the water depth, V(t, x) is the horizontal water velocity, they are unknown functions of two variables (t, x); C(x)∈C 2 ([0, L]) is the slope of the channel depth; g is a constant acceleration of gravity and k is a constant coefficient of friction.
For the numerical solution of the Saint-Venant equation, we propose an implicit upwind difference scheme. In this case of considering intervals on time and spatial variables, a uniform grid is used. The Lyapunov stability theory is used to establish the stability of the numerical solution of the exponential stability theorems proved by the authors of the pre sent work. The method of solving different systems of algebraic equations is the matrix sweep method. The theory of correctness and stability of matrix sweep algorithms is used [22].
In the course of conducting the numerical experiment, the following geological and hydrogeological data were obtained from the Big Almaty Canal under the Committee of Water resources management of the Ministry of Ecology, Geology and Natural Resources of the Republic of Kazakhstan and GoogleEarthPro application.
The numerical experiment was carried out on the Big Almaty Canal named after Dinmukhamed Akhmedovich Kunaev (Fig. 1). The canal is located on the territory of the Republic of Kazakhstan (Almaty region and Almaty) in the Ili River basin.
The canal waters are used for irrigation, large-scale industrial, small-scale economic and recreational needs. The length of the canal is 168 km. The route of the canal crosses the territories of three districts of the Almaty region: Enbekshikazakh, Talgar, Karasay and Almaty, its waters flow by gravity. The starting point of the canal is the Bartogai reservoir on the Chilik River (43°27¢11″ n. l. 78°23¢46″ e.l.), where the BAC begins. Then the canal crosses the Issyk, Talgar, Bolshaya and Malaya Almatinka rivers through aqueducts. According to the project, the final point of the BAC had to be the Kurtin reservoir on the Kurty River (43°19¢44″ n.l. 76°36¢55″ e.l.), but its waters do not always reach the last one due to seepage, evaporation, and water intake for household needs. The canal flow is fast. The banks of the canal are concreted, vegetation breaks through them in places.
For the calculation, a segment at the beginning of the canal with a length of 10 km is taken (Fig. 1). Therefore, the starting point of the calculation has the coordinates: 43°27¢11″ n.l. 78°23¢46″ e.l. and the endpoint has coordinates: 43°30¢20.14″ n.l. 78°19¢30.58″ e.l.Hydrogeological data from the BAC: 1. The length of the experimental section of the BAC is 10 km (CL0 -CL100).
2. The shape of the cross-section is trapezoidal. The terrain data for a section of 10 km from the starting point along the BAC, including the coordinates of the starting and ending points, the height above sea level of the points at every 30 m is obtained from the GoogleEarthPro application.
In the paper, we are given these values at the initial 14 and final 14 points with a step of 30 m ( Table 1). Full information can be obtained from the above program. Photos of some regions have an unprecedentedly high resolution.  (1), an upwind implicit difference splitting scheme is used for the lower terms.
We consider the flow in the river mode, so the values H * , V * must satisfy the inequality gH * >V *2 , only in this case λ 1 (x), λ 2 (x) will take positive values (λ 1 (x) > 0, λ 2 (x) > 0). In addition, the following equality for stationary solutions holds H * (x)V * (x) = Q * . Taking into account these conditions, we choose the following value of the functions H * (x), V * (x): Using the data in Table 1, we can calculate the function of the slope angle C(x) of the BAC. We have the values of the function B(x), measured at 335 points x j . The slope angle values: at each node x j can be calculated using the following difference ratio: Next, we are faced with the task of approximating the function C(x) most accurately from the given values C(x i ).
To do this, we use the cubic spline interpolation method.
For the correct statement of the problem, it is necessary to set the initial data along the entire length of the studied section of the riverbed: For calculating the unsteady flow in the canal of the Big Almaty Canal, only the subcritical flow Fr > 1 is considered, which requires setting one boundary condition at the upper and lower ends of the considered section.
According to Theorem 1 of [8], in the case of L = 10,000, we take the parameters of the boundary conditions from the following interval: Then the numerical solution of the difference initial boundary value problem is exponentially stable in L 2 -norm. In calculations, b 0 = -3, b 1 = 0.0001 satisfy these conditions. Therefore, we calculate the values of the corresponding parameters of the boundary conditions from these values k 0 = -0.2048437158, k 1 = -0.0001011587. Applying the rectangle method to calculate definite integrals, we obtain that r = 0.2048437158, s = -2.8227590592⋅10 -12 . The numerical experiment was carried out on a computer with the following technical characteristics: AMD Ryzen 9 3950X, 64 Gb, NVIDIA GeForce RTX 3090 ОС Windows 10 Pro.

1. Implicit upwind difference scheme
In the domain G = {(t,x): 0 ≤ t ≤ T, 0 ≤ x ≤ L}, we construct a difference grid with ∆t and ∆x steps on direction t and x, respectively. The anchor points of the difference grid (meaning the intersection of straight lines t = t k = k∆t and x = x j = j∆x) are denoted by (t k , x j ). The set of nodal points of the difference grid is denoted by G h , where And the values of the numerical solution at the anchor points are denoted by We select the steps of the difference grid ∆t, ∆x in such a way that the equalities K∆t = T and J∆x = L hold.
To find a numerical solution to a mixed problem (1)-(3) over the difference grid G h , we obtain the following upwind implicit difference scheme: with initial and boundary conditions in the following form: where ρ λ Here, using the discrete Lyapunov function, we obtain the exponential stability of the initial boundary value difference problem (5)-(7).

2. Formulas for the numerical solution of the Saint-Venant equation obtained by the matrix sweep method
Applying the matrix sweep method to (5)-(7), we obtain the following formulas for the numerical solution of equation (5) in the form: The algorithm of the matrix sweep method is presented in Fig. 3.
Using formulas (8) and (9) From the last formula, it can be seen that the height level and the horizontal velocity of water in the canal significantly depend on the stationary solution of the Saint-Venant equation.

3. Calculations of water level and horizontal velocity on the BAC section
Using the geological and hydrogeological data of the BAC, the following values of the height and velocity of the water flow into the BAC are obtained from the formula (11) (Tables 2, 3). In Fig. 4, 5, graphs of the numerical solution are presented. Based on the obtained data, a numerical experiment was conducted for a time of 5 seconds, 10 seconds, 60 seconds, 300 seconds ( Table 4). The results of the numerical experiments of computing the height and horizontal velocity of water flow presented in Tables 2, 3 are coherent with the actual measurements of water flow parameters in the BAC. The results presented in Table 4  .

Discussion of the results of the study of water level and horizontal velocity
Formulation of an implicit upwind difference scheme. As follows from Section 5. 1, we have constructed an implicit upwind difference splitting scheme for the lower terms (the initial boundary difference problem (5)- (7)) for the numerical solution of a mixed problem for the linear system of the Saint-Venant equation (problem (1)- (3)). Why exact ly such a difference scheme is proposed for the numerical calculation of the BAC? What are the advantages of this scheme compared to other difference schemes? First, all explicit difference schemes require the Courant-Friedrichs-Levy condition to be fulfilled. And this imposes strict restrictions on the steps of the difference grid from above. And the implicit scheme (5)- (7) is free from these restrictions. This gives us the opportunity to perform calculations for the BAC, the length of which can be any (depending on the computer performance parameters). In addition, we were able to prove the Lyapunov stability of this difference scheme.
The matrix sweep method for the Saint-Venant equation. It can be seen from Section 5. 2 that it is not an easy task to determine the numerical solution of an implicit upwind difference splitting scheme for the lower terms (5)- (7). First, the difference problem is presented in the form of a two-point scheme, and then, using a specially invented procedure, the two-point scheme is reduced to a threepoint one, in order to apply the standard matrix sweep me thod (8), (9). And then the numerical solution is determined by formulas (10), (11).
Numerical solution. Section 5. 3 presents the results of numerical calculation for the problem of numerical determination of the water level and velocity on the BAC section (10,000 m) located in the Almaty region of the Republic of Kazakhstan. Some initial data (hydrogeological data from the BAC-see above) were kindly provided to us by the administrators of the BAC, we are grateful to them for this. Based on these initial data, we were able to find a numerical solution to the system of Saint-Venant equations, which describes the movement of water to the BAC. Namely: -numerical values of the height of the water level H(t, x) in the BAC (Table 2); -the numerical value of the water velocity V(t, x) in the BAC (Table 3); -a graph of the numerical solution: a) a graph of the water height along the entire length of the considered section; b) a graph of the terrain of the considered area; (Fig. 4); -a graph of the numerical solution of the horizontal water velocity along the entire length of the considered section (Fig. 5).
Based on the data obtained above (the first 4 points), a numerical experiment was conducted for 5 seconds, 10 seconds, 60 seconds, 300 seconds (Table 4). Please note that the numerical values of the water level height and velocity given in Tables 2, 3 reflect only part of the numerical calculation results due to the limited space. However, they are fully represented in the graphs in Fig. 4, 5.
What we have as a result: -the numerical values of the height and horizontal velo city of the water flow presented in Tables 2, 3 fully correspond to the actual measurements of the water flow parameters in the BAC; -the results of the computational experiment presented in Table 4 show that the values of the steps of the difference grid ∆t, ∆x in numerical calculations can be arbitrary in the range 0 < ∆t < 1, 0 < ∆x < 1. Thus, the Courant-Friedrichs-Levy condition does not play any role here. This allows us to perform numerical calculations with large steps. Hence, we can find a numerical solution over the entire length of the BAC; -the results presented in Table 4 for the calculation time (Columns 3 and 8 of Table 4) show the effectiveness of the developed algorithm for the exponentially stable solution of the Saint-Venant equation; -the following limitations for the channel are inherent in this study: the bottom topography is rather smooth; there is no filtration through the bottom of the channel (concrete bottom); there are no additional inflows and outflows.
The subjects of further research: 1. To make calculations (find the numerical values of all the necessary parameters) along the entire length of the BAC.
2. To make calculations (find the numerical values of all the necessary parameters) along the entire length of the BAC