About one method of numeral decision of differential equalizations in partials using geometric interpolants

Abstract

The article presents a new vision of the process of approximating the solution of differential equations based on the construction of geometric objects of multidimensional space incident to nodal points, called geometric interpolants, which have pre-defined differential characteristics corresponding to the original differential equation. The incidence condition for a geometric interpolant to nodal points is provided by a special way of constructing a tree of a geometric model obtained on the basis of the moving simplex method and using special arcs of algebraic curves obtained on the basis of Bernstein polynomials. A fundamental computational algorithm for solving differential equations based on geometric interpolants of multidimensional space is developed. It includes the choice and analytical description of the geometric interpolant, its coordinate-wise calculation and differentiation, the substitution of the values of the parameters of the nodal points and the solution of the system of linear algebraic equations. The proposed method is used as an example of solving the inhomogeneous heat equation with a linear Laplacian, for approximation of which a 16-point 2parameter interpolant is used. The accuracy of the approximation was estimated using scientific visualization by superimposing the obtained surface on the surface of the reference solution obtained on the basis of the variable separation method. As a result, an almost complete coincidence of the approximation solution with the reference one was established.

1. Introduction

Traditionally, one of the possible results of the numerical solution of differential equations (DE) is a certain geometric model, the visualization of which allows you to visually evaluate the result. Thus, for most abstract solutions, there is a geometric interpretation. For example, the solution to an ordinary differential equation is a line, and the solution to the inhomogeneous heat equation of the rod is the surface compartment. Those, the result of solving the differential equation is a geometric object. Change the causal relationship to the inverse. Then it turns out that in order to solve the differential equation it is necessary to simulate some geometric object that has the required differential characteristics. A similar approach was implemented in [ 1, 2 ]. Of course, DE have a wide variety of varieties, and not for every differential equation there is an exact solution. Therefore, for the numerical solution of the differential equation it is enough that the required differential characteristics are provided at some discrete points (network nodes) that belong to the simulated geometric object. In this case, the intermediate values of the resulting solution will be determined using multidimensional interpolation. Then, to approximate the solution of the DE, it is convenient to immediately use one of the geometric interpolants.

A bit about geometric interpolant

A geometrical interpolant is a parameterized geometrical object passing through predetermined points, whose coordinates correspond to the initial experimentalstatistical information, or possessing the necessary, predetermined, properties. In accordance with the geometric theory of multidimensional interpolation [ 3-5 ], the geometric interpolant is formed by analytically describing the tree of the geometric model.

So, for a one-dimensional geometric interpolant (1parameter interpolant) the tree of the geometric model is just one line (Fig. 1), passing through the predetermined points.

Fig. 1. 1-parameter interpolant

In the BN-calculus [ 6-8 ], such an interpolant can be represented as the following point equation of a oneparameter set M of points:

n M = ∑ M i pi (u ), (1)

i=1 where M – the current point of the arc of the curve of the line passing through the predetermined points; Mi – the starting points through which the arc of the curve should pass; pi(u) – function of the parameter u; u – is the current parameter, which varies from 0 to 1; n – the number of starting points of the arc of the curve line; i – serial number of the starting point.

Moreover, the condition is that the one-parameter set belongs to the space of the selected dimension: n ∑ pi (u ) = 1 . This condition is satisfied using special i=1 algebraic curves obtained on the basis of the Bernstein polynomial [ 9 ]. The fulfillment of this condition is mandatory for all subsequent interpolants and is not given in the article below, since it is calculated in a similar way.

The point equation (1) is a symbolic notation. Having performed the coordinate-wise calculation for twodimensional space, we obtain a system of the same type parametric equations:  n xM = ∑ xMi pi (u );  yM = ∑i=n1 yMi pi (u ).

i=1

Similarly, any point equation for a space of any dimension can be represented as a system of parametric equations. Moreover, the presented system of parametric equations is an analytical description of the projections of the arc of a plane curve on the axis of the global coordinate system.

A two-dimensional geometric interpolant represents a two-parameter set of points – the surface of 3dimensional space passing through predetermined (Fig. 2).

The computational algorithm for determining a 2parameter geometric interpolant can be represented as the following sequence of point equations, which include m of 1-parametric interpolants at the stage of tree formation of the geometric model (Fig. 2):  n M1 = ∑ M1 j p1 j (u );  j=1 ...............................  n M i = ∑ M ij pij (u );  j=1 (2) ...............................  n M m = ∑ M mj pmj (u );  j=1 M = ∑im=1 M i qi (v), where qi(v) – function of the parameter v.

To describe a 2-parameter interpolant (Fig. 2), a 3dimensional Cartesian coordinate system is used (although the proposed equations are also valid for an affine coordinate system). In addition, such a geometric interpolant can exist in a space of higher dimensions. In this case, the point equation will remain unchanged, but when performing the coordinate-wise calculation of the parametric equations of the system, there will be more, and their number will directly depend on the dimension of the space in which the simulated geometric object is located.

Similarly, a three-parameter interpolant is defined by a 3-parameter set of points – a hypersurface of 4dimensional space passing through predetermined points (Fig. 3).

The computational algorithm for determining the 3parameter interpolant will include m of 2-parametric interpolants forming an even more extended tree of the geometric model (Fig. 3):  l M ij = ∑ M ijk pijk (u );  k =1 ............................... M i = ∑n M ij qij (v); (3)  j=1 ...............................  m M = ∑ M i ri ( w),  i=1 where ri(w) – function of the parameter w.

Summarizing this approach, we can obtain a geometric interpolant of n dimension corresponding to nparameter set of points or hypersurfaces (n+1)-th space passing through the predetermined points. At the same time, the belonging of the nodal interpolation points to the simulated geometric interpolant is ensured by the passage of all points through the guide lines (onedimensional interpolants) at each stage of the formation of the model tree: Fig. 1 → Fig. 2 → Fig. 3.

It should be noted that the geometric theory of multidimensional interpolation was developed and is effectively used to model and optimize multifactor processes and phenomena based on any experimental statistical information [ 10-12 ]. However, in the context of the above studies, it is used for a different purpose, namely, for solving DE.

To solve the equations of mathematical physics [ 13 ], the choice of a geometric interpolant depends primarily on the dimension of the Laplacian. So, for the numerical solution of the inhomogeneous heat equation with linear Laplacian ∂2U/∂x2 served as a two-parameter interpolant U=f(x,t) [ 2 ]. Then with a flat Laplacian ∂2U/∂x2 + ∂2U/∂y2 the use of a three-parameter interpolant is necessary U=f(x,y,t), and with spatial ∂2U/∂x2 + ∂2U/∂y2 + ∂2U/∂z2 four-parameter interpolant U=f(x,y,z,t).

3. Analytical description of geometric interpolants

For the analytical description of geometric interpolants, the point equations of algebraic curves arcs passing through the predetermined points obtained on the basis of Bernstein polynomials [ 9 ] are used. The need to determine such curves lies in the fact that when modeling multi-factor processes for each separate problem, it is necessary to solve systems of linear algebraic equations (SLAE) in determining the desired equation. To obtain a universal approach to modeling multifactor processes [35], it was necessary to obtain such equations of arcs of algebraic curves into which you can substitute any values of the points coordinates (both fixed and variable), and immediately obtain the desired result. For this, the SLAE solution process was laid down directly at the stage of curve modeling. As a result, we obtained the point equations of algebraic curves arcs passing through predetermined points, which are the main tool of the geometric theory of multidimensional interpolation and approximation.

It should be noted a very important distinguishing feature of the obtained equations. For point equations, the belonging of a geometric object to a space of a specific dimension is determined by the sum of functions of a parameter (condition to equation (1)), which must be equal to 1. The using of Bernstein polynomials made it possible to ensure that this condition is met regardless of the dimension of the space of the global coordinate system. The functions of the parameter are determined by the Newton binomial, which is expanded for the parameter and its complement to 1. By this, it provides the condition that the arc of the curve belongs to a specific space, regardless of its dimension. In other words, the obtained parametric equations of the arc of the curve can be used for a space of any dimension and, accordingly, for solving differential equations with a Laplacian of any dimension

Another important feature of the obtained equations of the curve arc is the uniform distribution of the parameter values, which was originally laid down in the method for determining the curve arc passing through predetermined points. Moreover, for each specific coordinate axis having a uniform distribution of the coordinates of the source points, a linear relationship between the natural value of the factor belonging to the projection axis and the current parameter is valid. This significantly reduces the amount of necessary calculations when approximating the solution of the differential equation, allowing us to consider them on a regular multidimensional network of points. Moreover, the method is universal in nature and without making any changes, it can be fully used for both regular and irregular network of points.

In this way, point equations of arcs of curves of 2–10 order, passing through 3–11 points, respectively, were obtained. For example: 1. The point equation of an arc of a curve of the 2nd order passing through 3 predetermined points:

M = M1u (1 − 2u ) + 4uuM 2 + M 3u ( 2u −1) , (4) where u = 1 − u - parameter addition u to 1. 2. The point equation of an arc of a third-order curve passing through 4 predetermined points: M = M1 (u 3 − 2, 5u 2u + uu2 ) + M 2 (9u 2u − 4, 5uu2 ) + M 3 ( −4, 5u 2u + 9uu2 ) + M 4 (u 2u − 2, 5uu2 + u3 ). (5)

General approach to the approximation of the solution of differential equations

The main idea of the proposed approximation method is that at the nodes of the selected interpolation network of points the condition of the original differential equation is satisfied. For its implementation, the following fundamental computational algorithm was formed: 1. Depending on the source differential equation, form a network of points of the required dimension and density, which will be the basis for creating the tree of the geometric model. 2. Select arcs of approximating curves for an analytical description of a geometric interpolant, thereby forming a computational sub algorithm. 3. Perform coordinate-wise calculation and, in the case of using a regular network of points, go from the parametric equations system of the geometric interpolant analytical description to its equation in an explicit form. 4. Enter the coordinates of the points corresponding to the initial and boundary conditions. 5. To differentiate the obtained equations and substitute them in the original DE. 6. Substitute parameter values at the nodal points, thereby forming a local system of linear algebraic equations (SLAE). 7. In the case of using piecewise approximation, we repeat the first 6 points of the computational algorithm several times, thus accumulating local SLAEs to form a global SLAE. 8. We solve the obtained SLAE and determine the necessary values at the nodal points of the interpolant. After that, we substitute the result of calculations in the approximation equation from the 5th point. 9. We analyze the result and check its reliability. In the case of insufficiently accurate results, we increase the number of nodal points of the geometric interpolant.

Of course, each engineering task is separate in nature and has its own characteristics, but this will not affect the fundamental approach to solving differential equation. For example, with a large number of nodal points, it is possible to use composite approximating curves that will form composite geometric objects in multidimensional space. If necessary, they can be docked with the required smoothness order [ 14-19 ]. And with an increase in the order of the differential equation, an increase in the order of the approximating curve is necessary. Moreover, in order to obtain the correct result of solving the differential equation, it is necessary that the order of the approximating curve be greater than the order of the original differential equation.

It should be noted that the result of the implementation of the proposed computational algorithm will be a general control solution. It can have an infinite number of particular solutions. A specific solution is distinguished from a variety of particular solutions using initial and boundary conditions, which, unlike most methods for solving differential equation, must be laid down in the form of input data at the stage of creation and analytical description of the geometric interpolant, thereby forming point 4 of the computational algorithm. In other words, the geometric display of the initial and boundary conditions are also some geometric objects: points, lines, surfaces, etc. Thus, the desired geometric interpolant must be a carrier of geometric objects corresponding to the initial and boundary conditions.

5. An example of approximation of the solution of the heat equation using a two-parameter interpolant

Consider the use of the proposed method on the example of solving the following inhomogeneous heat equation: ∂U =a2 ∂2U

∂x2 + 2x +1, ∂t 0 < x < 1, t > 0, U (0, t ) =1,

U (1, t ) = 2, U ( x, 0) = x +1.

To approximate the solution of equation (6), we use a 16-point 2-parameter interpolant [ 2 ]. Using the point equation of the arc of a third-order curve passing through 4 predetermined points, we obtain the following computational algorithm for determining a 16-point 2parameter interpolant, which is determined using the point equation (5) by the following sequence of point equations: indicate the corresponding coordinates of the points along the axis U: UM1111=UM1112=UM1113=UM1114=1,

Thus, it remains to determine the values of the geometric interpolant at 6 points: M1122, M1123, M1124, M1132, M1133 and M1134 that correspond to the following values of the parameters of the nodal points of the interpolant: =32; xM1124 =1; =23; xM1134 =1.

(1− e−(π na)2 t ) sin (π nx).

For a visual comparison of the obtained results, we will visualize the obtained surfaces and superimpose them on each other (Fig. 4). In this case, the green solution shows the reference solution obtained by the method of separation of variables. As can be seen from Figure 4, with the help of the 16point geometric interpolant, it was possible to achieve an almost complete degree of coincidence with the reference solution. Moreover, further use of the obtained polynomial equation for engineering calculations is more preferable in comparison with the equation obtained by the method of separation of variables. It should be noted that, if necessary, the number of nodal points of the approximating network can be practically any and can always be increased to achieve the required accuracy of the solution.

6. A generalization of the proposed solution of the inhomogeneous heat equation to a multidimensional space

Let us consider a generalization of the proposed solution of the inhomogeneous heat equation for a higher-dimensional Laplacian. In this case, the computational algorithm does not have fundamental differences. Only increases the dimension of the geometric interpolant and the number of equations of coordinate calculation. Based on this, we consider not a particular, but a general solution of the heat equation, given in a general form for a three-dimensional Laplacian: ∂U ∂t =a2  ∂2U ∂2U ∂2U 

 ∂x2 + ∂y2 + ∂z2  + f ( x, y, z ).

As a geometric interpolant, we choose a 4-parameter hypersurface belonging to a 5-dimensional space. As an example, let us take a curve of the third order passing through 4 forward given points (5) as an approximating arc. Then the computational algorithm for determining the 4-parameter interpolant takes the following form: M1 = M11 ( w3 − 2, 5w2w + ww2 ) + M12 (9w2w − 4, 5ww2 ) + M13 (−4, 5w2w + 9ww2 ) + M14 ( w2w − 2, 5ww2 + w3 ) , M 2 = M 21 ( w3 − 2, 5w2w + ww2 ) + M 22 (9w2w − 4, 5ww2 ) + M 23 (−4, 5w2w + 9ww2 ) + M 24 ( w2w − 2, 5ww2 + w3 ) ,  M 3 = M 31 ( w3 − 2, 5w2w + ww2 ) + M 32 (9w2w − 4, 5ww2 ) + M 33 (−4, 5w2w + 9ww2 ) + M 34 ( w2w − 2, 5ww2 + w3 ) , M 4 = M 41 ( w3 − 2, 5w2w + ww2 ) + M 42 (9w2w − 4, 5ww2 ) + M 43 (−4, 5w2w + 9ww2 ) + M 44 ( w2w − 2, 5ww2 + w3 ) , M = M1 (ϕ 3 − 2, 5ϕ 2ϕ +ϕϕ 2 ) + M 2 (9ϕ 2ϕ − 4, 5ϕϕ 2 ) + M 3 (−4, 5ϕ 2ϕ + 9ϕϕ 2 ) + M 4 (ϕ 2ϕ − 2, 5ϕϕ 2 +ϕ 3 ) ,  where w = 1− w and ϕ = 1−ϕ .

It should be noted that the sequence (10) did not include 16 2-parameter interpolants Mij, which also must be determined by analogy with the sequence (7). Thus, the desired geometric interpolant will pass through 256 nodal points. Accordingly, for the inhomogeneous heat equation with a flat Laplacian ∂2U/∂x2 + ∂2U/∂y2 the number of nodal points will be 64.

We perform the coordinate-wise calculation of the sequence of equations (10) for a 5-dimensional space. To do this, we adopt a Cartesian coordinate system with axes: x, y, z, t and U. Given the special properties of the (9) (10) where ax, ay, az, at, bx, by, bz, bt – parameters that are determined depending on the initial and boundary conditions for solving the differential equation.

Further, taking into account the linear dependence of the first 4 equations of system (11), we proceed to the equation given explicitly U=f(x,y,z,t). We differentiate it in accordance with equation (9) and, substituting the parameter values at the nodal points of the interpolant one by one, we compose a SLAE, solving which we obtain the desired numerical solution of the inhomogeneous heat equation.

Similarly, other arcs of curves passing through forward given points of a higher order or obtained in some other way can be used to approximate the solution of the differential equations. It is also possible to create mixed geometric interpolants, including arcs of curves of various orders. Thus, the number of nodal points can be any at each separate stage of the formation of the geometric interpolant and depends primarily on the initial and boundary conditions of the differential equation.

7. Conclusion

A method for the numerical solution of differential equations using a geometric interpolant is proposed. Moreover, it can easily be generalized to multidimensional space and therefore can be used to solve differential equations with a large number of variables, by analogy with the geometric modeling [2021] of multifactor processes and phenomena [ 3-5 ]. The proposed method is considered as an example of solving the inhomogeneous heat equation using a 16-point twoparameter interpolant. In this case, a generalization of the proposed solution of the heat equation to multidimensional space is made. In a similar way, the number of nodal points of a geometric interpolant can be increased, which allows you to geometrically simulate the solution of differential equations with any predetermined accuracy. For this, not only arcs of curves passing through predetermined points can be used, but also contours of the required smoothness order. Also, the proposed approximation method can be effectively generalized not only in the direction of increasing the dimensionality of space, but also in the direction of increasing the order of the initial differential equation, which is the prospect of further research.

The geometric theory of multidimensional interpolation can also be used to solve other engineering problems of modeling and visualization multi-factor processes and phenomena [ 20-30 ].

Acknowledgments

This work was completed and published with financial support from the Russian Foundation for Basic Research, grant 19-07-01024.

References

[1] Konopatskiy E.V. The solution of differential equations by geometric modeling methods . Proceedings of the 28th International Conference on Computer Graphics and Machine Vision GraphiCon 2018. September 24-27 , 2018 . Tomsk: TPU, 2018 . pp. 322 - 325 .

[2] Konopatskiy E.V. Modeling an approximating 16- point compartment of the response surface as applied to solving the inhomogeneous heat equation . Geometry and Graphics . - Moscow: Infra-M, 2019 . Vol. 7 . No. 2 . pp. 38 - 45 . DOI: 10 .12737/article_5d2c1a551a22c5. 12136357 .

[3] Konopatskiy E.V. An approach to the construction of geometric models of multifactorial processes of multidimensional interpolation . Software Engineering . Moscow: 2019 . Vol. 10 . No. 2 . pp. 77 - 86 .

[4] Konopatskiy E.V. The principles of building computer models of multifactor processes by the method of multidimensional interpolation . Proceedings of the II International Scientific and Practical Conference: “Software Engineering: Methods and Technologies for the Development of Information and Computing Systems” (November 14-15 , 2018 ). Donetsk: DonNTU , 2018 . pp. 277 - 287 .

[5] Konopatskiy E.V. Geometric modeling and optimization of multifactor processes by the method of multidimensional interpolation . Proceedings of the International scientific conference on physical and technical informatics CPT2018 (May 28-31 , 2018 ). Moscow-Protvino, 2018 . pp. 299 - 306 .

[6] Balyuba I.G. Constructive geometry of manifolds in point calculus: dis . ... Dr. tech. Sciences: 05.01 .01. I.G. Balyuba. Makeevka , 1995 . 227 p.

[7] Balyuba I.G. , Naidysh V.M. Point calculus . Melitopol: Melitopol State Pedagogical University named after B. Khmelnitsky , 2015 . 236 p.

[8] Bumaga A.I. , Konopatskiy E.V. , Krysko A.A. , Chernisheva O.A. Introduction to the mathematical apparatus of BN-calculus. Materials of the VII International Scientific and Practical Internet Conference "Problems of the quality of graphic preparation of students in a technical university: traditions and innovations." Perm: PNIPU , 2017 . Issue. 4. pp. 76 - 82 .

[9] Konopatskiy E.V. Modeling of arcs of curves passing through predetermined points . Bulletin of computer and information technology. Moscow: 2019. No. 2 . pp. 30 - 36 . DOI: 10 .14489/vkit. 2019 . 02 .pp. 030 - 036 .

[10] Konopatskiy E.V. , Bezditnyi A.A. Geometric modeling and optimization of multidimensional data in Radischev integrated drawing . IoP conference series: Journal of Physics: Conf. Series 1260 ( 2019 ) 072006 . - DOI: 10.1088/ 1742 -6596/1260/7/072006.

[11] Konopatskiy E.V. , Bumaga A.I. , Krysko A.A. , Voronova O.S. Geometric modeling and optimization of physico-mechanical properties of tar polymer. Information technologies in design and production . Moscow: Compass, 2019 . No. 1 ( 173 ). pp. 20 - 24 .

[12] Voronova O.S. , Konopatskiy E.V. Geometric modeling of the parameters of the physical state of water and water vapor . Bulletin of cybernetics . Surgut: Surgut State University, 2019 . No. 1 ( 33 ) 2019 . pp. 29 - 38 .

[13] Chuprov I.F. , Kaneva E.A. , Mordvinov A.A. The equations of mathematical physics with applications to the problems of oil production and pipeline transport of gas: Textbook . Ukhta: USTU, 2004 . 128 p.

[14] Konopatskiy E.V. Principles of modeling multifactor processes with a large amount of input data. Information technology in design and production . Moscow: Compass, 2018 . No. 4 ( 172 ). pp. 20 - 25 .

[15] Figueira G. , Almada-Lobo B . Hybrid simulationoptimization methods: A taxonomy and discussion . Simulation Modelling Practice and Theory . 2014 . Vol. 46 . pp. 118 - 134 .

[16] Foster I. , Zhao Y. , Raicu I. , Lu S. Cloud Computing and Grid Computing 360-Degree Compared // eprint arXiv: 0901 .0131, 2008 http://arxiv.org/ftp/arxiv/papers/0901/0901.0131.pdf

[17] Gallopoulos E.N. Houstis and J.R. Rice , «Computer as Thinker/Doer: Problem-Solving Environments for Computational Science» IEEE Computational Science &amp; Eng., Vol. 1 , No. 2, Summer 1994 . pp. 11 - 23 .

[18] Szalay A . Extreme data-intensive scientific computing . Computing in Science & Engineering , 2011 . Vol. 13 . No 6. pp. 34 - 41 .

[19] Zadeh L. Fuzzy logic, neural networks and soft computing . Соmmutation on the ASM-1994 . Vol. 37 . No 3. рр. 77 - 84 .

[20] Voloshinov D.V. Constructive geometric modeling: theory, practice, automation . Saarbrucken: Lambert Academic Publishing, 2010 . 355 p.

[21] Voloshinov D.V. , Solomonov K.N. Constructive geometric modeling as a perspective for teaching graphic disciplines . Geometry and graphics , 2013 . Vol. 1 . No. 2 . pp. 10 - 13 .

[22] Gerget O.M. Bionic models for identification of biological systems . Journal of Physics: Conference Series , 2017 . Vol. 803 . No. 1 . pp. 12 - 24 .

[23] Zakharova A.A. , Korostelyov D.A. , Fedonin O.N. Visualization Algorithms for Multi-criteria Alternatives Filtering ( 2019 ). Scientific Visualization 11 .4: 66 - 80 , DOI: 10.26583/sv.11.4.06.

[24] Zakharova , A. , Shklyar , A. Basic principles of data visual models construction, by the example of interactive systems for 3D visualization ( 2014 ) Scientific Visualization, 6 ( 2 ), pp. 62 - 73 .

[25] Zakharova , A.A. , Nebaba , S.G. , Zavyalov , D.A. Algorithms and Software Increasing the Efficiency of Processing Multidimensional Heterogeneous Data. Programming and Computer Software , 2019 . 45 ( 4 ). pp. 196 - 201 . doi: 10 .1134/S0361768819040108.

[26] Filinskikh A.D. , Merzlyakov I.N. Assessment of geometric models based on the structure of their parameters . Information measuring and control systems , 2015 . Vol. 13 . No. 3 . pp. 69 - 74 .

[27] Bondarev A.E , Galaktionov V.A. Parametric Optimizing Analysis of Unsteady Structures and Visualization of Multidimensional Data . International Journal of Modeling, Simulation and Scientific Computing , 2013 . Vol. 04 , No. 01 , 13 p., DOI 10 .1142/S1793962313410043.

[28] Bondarev A.E. On visualization problems in a generalized computational experiment ( 2019 ). Scientific Visualization 11 .2: 156 - 162 . DOI: 10 .26583/sv.11.2.12.

[29] Bondarev A.E. , Galaktionov V.A. Construction of a generalized computational experiment and visual analysis of multidimensional data . CEUR Workshop Proceedings , Vol. 2485 , 2019 , pp. 117 - 121 . DOI: 10 .30987/graphicon-2019-2- 117 -121.

[30] Alekseev A.K. , Bondarev A.E. Estimation of the Distance between True and Numerical Solutions . Computational mathematics and mathematical physics , 2019 . Vol. 59 . No. 6 . pp. 857 - 863 . DOI: 10 .1134/S0965542519060034.