UDC 519.248:
GRID APPROXIMATION METHODS FOR COMPLEX EVENT SYSTEMS
O.Yu. Vorobyov, O.Yu. Tarasova, A.N. Ovsyannikova
A newly introduced method of grid approximation of an unknown distribution of a set of random events is presented.
Introduction
In view of a large number events in real statistical systems, it becomes difficult to determine the states in which the system may find itself. One way to overcome difficulties of this kind is to find a grid distribution of the state of the system, defined on the introduced grid and close to the desired unknown distribution. Matrix marketing models prompted the idea of methods for approximating distributions, similar to grid methods for solving differential equations.
Formulation of the problem
The essence of the grid method is as follows: instead of the initial space of elementary events, its grid analogue is introduced. This grid model is described by probabilities that are defined only on grid events. Unknown distributions, i.e. the laws in accordance with which the space of elementary events evolves are replaced by the corresponding grid counterparts. As a result, the original problem is replaced, or, as they say, is approximated by a system of grid distributions - a grid scheme. In other words, the approximated set of events is “meshed” in order to try to imagine how the “unknown” events from the approximated set behave within the “known cells”, the role of which will be given to the events-terraces that form the grid.
Basic concepts of eventology and probability theory
Eventology is a new direction that has arisen within the framework of probability theory and studies the distribution of sets of events, the structure of dependencies of sets of events.
Definition 1.1. A probability space is a triple (C1, 7% P), where O is the space of elementary events, T7 is the algebra of events, and P is the probability defined on the elements of the algebra P - random events x, y, ... eP.
Definition 1.2. finite set selected events AeF selected from the probability space algebra (O, P, P) and consisting of ^Y = |x| events is called the set of random events.
Definition 1.3. The set of random events X generates various sets of so-called terrace events, among which there are terrace events for lcX in the form of an intersection: ler(X)= P|x P)*c.
Definition 1.4. Two random events x, y e X (x φ 0, y φ 0) are called nested if only two relations are possible between them
that is, one of these events is nested within the other: X C y OR V C L".
Event grid, or eventological grid (E-grid)
The initial stage in the construction of a grid scheme is the replacement of the initial space of elementary events with some grid of events that form its partition.
Definition 2.1. An eventological grid (E-grid) is a set of 5 s / 7 non-intersecting random events selected from the algebra P of the probability space (O, P, P) and forming a partition of the space of elementary events O ..
E-grid 5 (grid set of random events 5, grid eventological distribution of set 5) is defined if and only if:
1) a set of BarR non-intersecting random events is chosen, forming a partition
2) a set of probabilities is given<7(5) = Р(«), 5 е
For a "one-dimensional" problem, the simplest example of an E-grid is the equiprobable division of the space of elementary events into N equiprobable events, the probability of which is equal to WE (equiprobable E-grid). Equally probable events of the E-grid are called E-terraces of the grid (E-network nodes), and their probabilities are also called E-grid steps. The set of E-terraces forms a set of events where E-grid distributions are defined.
Definition 2.2. The eventological grid of the n-th order 8" of the space of elementary events O is the Minkowski intersection of the partitions A1,..., An with P:
Despite the apparent simplicity, the question of choosing an E-grid deserves attention. On the one hand, it is desirable to take a large number of terrace events, i.e. use fine, detailed E-grids. By more accurately conveying the region of change of the E-argument, we intuitively expect to better approximate the desired E-distribution by grid distributions. On the other hand, practical considerations, and primarily the limited speed and memory capacity of computers, forces us to turn to E-grids with a relatively small number of E-terraces. Non-equiprobable E-grids often serve as a solution to this problem. If there is information about the E-distribution, for example, the “location” in the space of elementary events of some of its features is known, for the “resolution” of which a fine E-grid is needed, then it is possible, without increasing the total number of terraces, to thicken the grid in the “neighborhood” of these features , and make the grid sparse in the “smooth” distribution area.
Definition 2.3.. An approximating 5”-grid mapping is called a mapping cp: B" -» 2x, which on the E-cell
Definition 2.4. - an approximation of the set of events X is the set of co-
E-grid elements - view terrace events
They are called E-cells of the E-grid b1”. It's clear that
E-cells form a partition of the space of elementary events O. Thus, the E-grid has the form:
<Р(-ІСГ{а".а"]^(ХЄХ: *П%>..a«)*0b
which is abbreviated
where every event
Xph \u003d (xf: 1b!}
is called a 5”-approximation of the event xx X . The terrace event will be denoted
u!G*(L")= GK P(^)C’ X^X
Let us give examples of an X-approximating 5 mapping with different degrees of error: a
Y x, y "L X X
G Y Y> 7~ x, y, r x, r X
Y Y > 2 g<т 0
Rice. 1. Eventological grid of the second order V2 = (a.b.c, c!, r)(n)>a, [), y, e, s), approximating the initial E-distribution of the triplet of events A" = (x, y ,d) with zero error (on the left).
X<р -{хР"уф"г^} на э-ячейках (справа)
Rice. 2. Eventological grid of the second order 5 = (a, b, c, c/, e)(n)(a, /?, y, s, e), approximating the initial E-distribution of the triplet of events X = (x, y, d) with some error (on the left). E-grid approximation X "r \u003d (x‘r ^ , x1 ') on E-cells (right)
Types of E-grid approximation
The accuracy of the E-grid approximation of the E-distributions is determined by the structure of the local dependence of the approximated set of events, which is characterized by the ratio of the dependency structures of the two sets of events involved in the E-grid approximation: the E-grid 51" and the approximated set X .
Let us consider three types of E-grid approximation, each of which is valid under one of the three assumptions about the local structure of the dependence of the approximated set of events X: nested, independent, and least intersecting.
Rice. 3. Eventological grid of the second order S = (a, b, c, d, e)(n)(a, P, y, S, e), approximating the initial E-distribution of the triplet of events X - (x, y, z) with the largest error (left).
E-grid approximation X9 = (Q,Q,Q) on E-cells (right)
Definition 3.1. The set of events X has a locally nested structure with respect to
body of the E-grid Sn, if for any E-cell ter, - C] a" e Sn and any event
xnter, "=<,
(a "a) lteV..vv
those. or ter | a, is contained in l, or does not intersect with it.
Two properties follow from this definition:
U (*ntev..v>J=teV..vv
2) events xn ter, „ _, where xe X , „ are equiprobable.
(a ...a ) a ... a 1 1
Definition 3.2. The set of events X has a locally least intersecting structure with respect to the E-grid Sn if for any E-cell ter , „ eSn and any set
set of events X , „ with X:
I) X (*nte^.aJ=teVv)’
2) the events xn ter, „ 1, where r e X , „ are equiprobable.
"(a ... a) a ... a 1
Definition 3.3. The set of events X has a locally independent structure with respect to the E-grid Sn if for any E-cell ter, "i e S" and any set of events
i „I „nil 1 L ()
(a...a) / (a...a)
2) events хгМег, „, where ххХ, „ are equally probable,
7 (a...a) a...a 1
3) the events x n leg(a, an, where x ∈ X^ an, are independent in the aggregate.
Theorem 3.1. Eventological distribution of the 5th -approximation of the set of events X
p "p (X) ^ RCegf (X)), XaX
takes one of the following forms:
1) for a locally nested structure of the set of events X:
p*(x)= X p(4er(v"...*")>'
ХШ*г(а>а„))
2) for the locally least overlapping structure of the set of events X:
XZ
3) for a locally independent structure of the set of events X:
p1p(X) = £ (P(x n ter(a, (P(ter(a, - P(x n ter!a, qH)))1A"11cX,
X<^<р{\ег, „)
where (p(xn 1er(t)1 wi,)) = 1 -(1 -P(1er(a, a n])\ln. Proof. It is obvious that
X (ter^ (X) n ter
ter. n eS V (0x ap) -
Substituting (1) into (2), we obtain
ter-(T)= IV P<^)“ =
XЄ(p\ ter i ^ ^ a1...a"
Y (xnteral v) n X Y ntera\..a")
xeg>\ ter , „ xєХс XH(p\ ter 1 i
V V a ■ a) J V ° a) y
Using the obtained expression for Xxx9 (X), we obtain
a, ^ = n „in-,) A 1 * cn%> .. "",)
heh x<Ехс
or in other notation
ler"(DT)nx(a, t = T) (*nem(a,
where X = ter ,
\ x n ter, „ I is the complement of x to ter, | e S'1, and xc = Q \ .
o 1 v (a...a [ / (a...a! V
(a ". an)! "............" V
addition X to P.
Thus, the expression under the sum sign in (3) can be written as
and formula (3)
fl^ntev....")) P (*ntev v/
p(X)= £ p(ter^(X)nter , a)=
I P P (^nter!a>..o",) P (JCnteV..V>>C"
Formula (4) is the eventological distribution of the Sn -approximation of the set of events X , which has an arbitrary local structure. Using formula (4), we can write formulas for the eventological distribution S" -approximation of the set of events X, which has a locally nested, locally least overlapping, and locally independent structure.
For example, consider a set X that has a locally least overlapping structure. By the definition of a locally least intersecting structure, the events d:nter , „ , where xf/?(ter. , u) do not intersect and have the same probability. So
\0 ... O (\C1 ... (I)
the eventological distribution of S"-approximations of the set of events X with a locally least intersecting structure can be written as:
As an approximate set X, consider the set of strategies offered by the Arthur D. Little company. The distribution of the set of strategies is unknown, so we will use the grid approximation method to find this distribution.
In the classical presentation of the marketing model AOL / LC, the principle of constructing the AEL matrix suggests the type of eventological grid. The intersection of two sets of events - four events of the life cycle of the industry and five events of its competitive position form
eventological grid of the second order. V2 = (a, b, c, c1)(u)(a, /3, y, 3, e). Each cell of the grid corresponds to a set of strategies from a set of basic strategies offered by AOB specialists as a guide to action. Having made assumptions about the local dependence structure of the approximated set of strategies A" (for example, based on the analysis of these strategies
marketing specialists), we obtain the eventological distribution of the 52-approximation of the set of events X .
The following results were obtained in the work:
1. The concepts of the eventological theory of grid methods are defined: E-grid, X-approximating $n grid mapping, 8P-approximation of the set of events.
2. Three types of local dependency structures are defined.
3. A theorem on the eventological distribution of the 5"-approximation of the set of events X is formulated and proved for various local dependency structures of this set.
1. Efremov, B.C. Classical models of strategic analysis and planning: ADL / LC / B.C. Efremov // Management in Russia and abroad. M.: Finpress. - 1998. - No. 1 (http://www.cfm.ru/press/rnanagment/1998-l/09.shtml).
2. Vorobyov, O.Yu. Introduction to eventology / O.Yu. Vorobyov - Krasnoyarsk: KrasGU, INM SO RAN, 2005.-512 p.
3. http://www.r-events.narod.ru
Xs<р(\£Х - „) X ,
(a -a > a...a
where the summation is over all ter , „<= 5”, содержащим X с X .
Application of the results
Conclusion
Literature
Consider the differential problem in operator form (2.1) and the operator form of finite difference scheme (2.3).
Let us introduce the norm of the grid function using the expression
Definition 3. The finite-difference scheme (2.3) approximates the differential problem on the exact solution if any norm of the difference (not necessarily in the form (2.9)) tends to zero at:
Definition 4. Finite-difference scheme (2.3) approximates the differential problem on the exact solution with order p in time and order q in space if any norm of the difference satisfies the equality
Thus, if a finite-difference scheme approximates a differential problem, then we are talking about the closeness of the differential and finite-difference operators at the grid nodes.
From the definition of the order of approximation, it is clear that the higher the order of approximation, the better the finite-difference scheme approaches the differential problem. This does not mean that the difference scheme solution can be as close to the solution of the differential problem, since the difference scheme can be conditionally stable or completely unstable.
To find the order of approximation, we use the Taylor series expansion apparatus for exact (unknown, but differentiable) solutions of the differential problem at grid nodes (we emphasize: the values of the grid function are discrete, therefore, not differentiable and therefore cannot be expanded into Taylor series).
In accordance with the definition of the order of approximation, we analyze the order of approximation of the finite-difference scheme (2.6), for which we write this scheme on the exact solution :
(2.12)
Let us expand the values in Taylor series in the variable x in the vicinity of the node up to the fourth derivative inclusive, and the value in Taylor series in the variable t in the vicinity of the node up to the second derivative inclusive, we obtain
(2.15)
Substituting (2.13)-(2.15) into (2.12), we find
In this way,
i.e., the explicit scheme (2.6) for the heat equation has the first order of approximation in time and the second in the spatial variable. Similarly, the same order of approximation can also be obtained for the implicit scheme (2.8).
Definition 5. The solution obtained using the finite difference scheme (2.3) converges to the exact solution U if any norm of the difference tends to zero as the grid characteristics tend to zero:
Definition 6. Finite-difference scheme (2.3) has the pth order of convergence (order of accuracy) in time and the qth order of convergence in space if any norm of the difference satisfies the equality
Thus, the order of convergence (order of accuracy) characterizes the closeness of the finite-difference and exact (unknown) solutions.
2.1.3 Studying the stability of finite difference schemes
Let the input data in the finite-difference scheme (2.3) be perturbed and take values. Then the grid function will also receive a perturbation and take the value
Definition 7. A finite-difference scheme (2.3) is stable with respect to the input data if there is such a bounded constant independent of the grid characteristics of the input data that the inequality
Thus, the concept of stability is interpreted as follows: a finite difference scheme is stable if, for small perturbations of the input data (initial boundary conditions and right-hand sides), the finite difference scheme (2.3) ensures small perturbations of the grid function , i.e., the solution using a finite difference scheme is under the control of the input data.
If the input data includes only initial conditions, or only boundary conditions, or only right-hand sides, then one speaks of stability with respect to initial conditions, boundary conditions, or right-hand sides, respectively.
Definition 8. Finite-difference scheme (2.3) is called absolutely stable if inequality (2.17) is satisfied for any ratio of steps u.
Definition 9. Finite-difference scheme (2.3) that is unstable for any ratio of steps is called absolutely unstable.
Definition 10. A finite-difference scheme (2.3) is said to be conditionally stable if inequality (2.17) is satisfied for the grid characteristics u, which are subject to certain restrictions.
Since stability is one of the main characteristics of finite difference schemes, this section discusses various methods for studying the stability of finite difference schemes with respect to initial conditions. The most common methods for studying stability are the following:
Method of harmonic analysis (Fourier);
The principle of maximum;
Spectral method;
energy method.
Each of these methods has advantages and disadvantages.
Method of harmonic analysis. It is known from mathematical physics that the solution of initial-boundary value problems is represented as the following series:
where are the eigenvalues, and are the eigenfunctions obtained from the solution of the corresponding Sturm-Liouville problem, i.e. the solution can be represented as a superposition of individual harmonics, each of which is the product of a function of time t and a function of the spatial variable x, the latter in terms of modulus is limited from above by one for any values of the variable x.
At the same time, the time function , called the amplitude part of the harmonic, is not limited in any way, and, in all likelihood, it is the amplitude part of the harmonics that is the source of the growth of the function uncontrolled by the input data and, therefore, the source of instability.
Thus, if the finite-difference scheme is stable, then the ratio of the amplitude part of the harmonic in the upper time layer to the amplitude part in the lower time layer must be less than unity in absolute value.
If we expand the value of the grid function into a Fourier series in terms of eigenfunctions:
where the amplitude part can be represented as a product
Un is the dimensional and constant factor of the amplitude part, and k is the exponent (corresponding to the number of the time layer) of the time-dependent factor, then, by substituting (2.18) into the finite-difference scheme, we can estimate the ratio of the amplitude parts on neighboring time layers by modulo.
However, since the summation operation is linear and the eigenfunctions are orthogonal for different summation indices, it suffices to substitute one harmonic of expansion (2.18) instead of the grid values into the finite-difference scheme (remove the index n from the amplitude part), i.e.
Thus, if the finite difference scheme is stable with respect to the initial data, then
i.e., condition (2.21) is a necessary condition for stability.
Investigation of stability by the method of harmonic analysis of explicit and implicit schemes for the heat equation. We substitute expressions (2.20) into the explicit finite-spaced scheme (2.6) for the heat equation, we obtain
Here we use the formula following from the Euler formula:
and the formula , moreover, since and.
In accordance with (2.22), we obtain the expression
,
,
or, taking into account (2.21), the inequality
From this we obtain the following two inequalities:
from which the right one is always satisfied, and the famous Courant stability condition follows from the left one, or the more stringent condition for
It follows from (2.23) that the explicit scheme for the heat equation is conditionally stable with condition (2.23) imposed on the grid characteristics and h.
Let us now substitute the harmonics (2.20) into the implicit finite-difference scheme (2.8) for the heat equation, we obtain
,
always, since a and the square of the sine are greater than zero.
Consequently, the implicit scheme for the heat equation is absolutely stable, since no restrictions were imposed on the grid characteristics h to satisfy the inequality.
The complex is called the Courant number for the heat equation.
Maximum principle. In mathematical physics, the principle is known, according to which the solution of the initial-boundary value problem inside the computational domain cannot exceed the values of the desired function on the space-time boundary. This principle underlies the method for studying the stability of finite-difference schemes, called the maximum principle.
To use it, consider the explicit finite difference scheme (2.6) for the heat equation in the form
and introduce the norm of the grid function in the form
Then from (2.24) we obtain
then from (2.25) we have the inequality
whence, continuing the chain of inequalities up to the initial condition, we obtain
where is the initial condition from (1.18).
Inequalities (2.27) in computational mathematics are called the maximum principle. It is a sufficient condition for the stability of the explicit scheme (2.24) for the heat equation.
Thus, if the Courant condition (2.26) is satisfied, then chain (2.27) shows that the value of the grid function on any time layer does not exceed the initial condition in the norm, i.e., the scheme under consideration is stable with respect to the initial condition, and condition (2.26) is now not only necessary in accordance with the method of harmonic analysis, but also sufficient.
Spectral method for studying stability. Consider the grid functions and,, on two time layers and represent the finite difference scheme in the following operator form:
where S is the transition operator from layer tk to layer tk+1. Such an operator cannot be constructed for any finite difference scheme
(for example, the sweep method cannot be represented in the form (6.64)). For explicit finite difference schemes (2.6), the operator S is represented by the following transition matrix:
.
We compose the norm operation from the left and right parts of equality (2.28) and use the property of the norm: the norm of the product of operators does not exceed the product of norms, we obtain
If an inequality of the form
then conditions (2.29) and (2.30) imply the maximum principle
Thus, if the scheme is stable, then the norm of the transition operator S does not exceed 1, and hence condition (2.30) is a necessary condition for the stability of finite difference schemes.
Energy method for studying the stability of finite difference schemes. As can be seen from the previous sections, the harmonic analysis method and the spectral method are necessary conditions for the stability of finite difference schemes, and the maximum principle is a sufficient condition for stability. In this paragraph, one of the most powerful and widespread methods is considered - the energy method developed in the works of A.A. Samarsky and based on the concepts of energy space with energy norm, energy identity (inequality) and the maximum principle. It will be shown below that the conditions used in the energy method are sufficient conditions for the stability of finite difference schemes.
To understand the energy method, consider its application to study the stability of finite difference schemes in the numerical solution of the following first initial-boundary value problem for the heat equation with homogeneous boundary conditions:
T > 0; (2.31)
X = 0, t > 0; (2.32)
X = 1, t > 0; (2.33)
T = 0; (2.34)
On grid (2.2), we will approximate this problem using explicit (2.6) and implicit (2.8) finite-difference schemes written in vector-operator form as follows:
where the finite difference operator approximates the differential operator with respect to the space variable x, i.e.
Energy space. Let us introduce the energy space HA of grid functions , which is a Hilbert space in which the scalar product for two elements is defined:
,
(2.37)
and hence with the norm
.
(2.38)
As you know, a Hilbert space is a complete normed space in which the scalar product is defined. Here, completeness is defined in the sense that if a sequence of grid functions converges to its limit at (in this case, to the solution of a differential problem), then it is fundamental, i.e. the Cauchy condition is satisfied
Indeed, if the finite difference scheme approximates the differential problem and is stable, then by the equivalence theorem the solution using the finite difference scheme converges to the solution of the differential problem when the grid is refined.
For two grid functions u (two elements of the Hilbert space HA) on different grids with steps hn and hm, the concept of completeness means that when the grid is refined, i.e., when (or) the sequences converge to the same limit, i.e., (2.39).
In what follows, we will need the following concepts characterizing finite difference operators: conjugacy, self-adjointness, and positive definiteness.
Definition 11. A finite difference operator A* is called adjoint to the operator A if the equality
For example, if the operator A is a symmetric matrix with real elements (A = AT), then A is the adjoint operator (this can be verified directly).
Definition. A finite difference operator A is called self-adjoint if the equality
Definition 12. A finite difference operator A is called positive definite or positive semidefinite
on the Hilbert space of grid functions, if
It can be shown that the difference operator is self-adjoint, i.e.,
In order to determine the eigenfunctions and eigenvalues of the finite difference operator A, we first consider the problem for the eigenvalues and eigenfunctions of the operator :
Native functions must satisfy the following conditions:
be orthogonal on the segment at, i.e., for example;
satisfy homogeneous boundary equations (2.32), (2.33);
their number must match the number of eigenvalues.
These conditions are satisfied by the functions:
To find the eigenvalues, we substitute (2.44) into (2.43), we get
It can be seen from (2.45) that all eigenvalues of the operator A are negative, while the eigenvalues of the operator are positive, i.e., the operator is positive definite.
In studying the stability of the explicit finite-difference scheme (2.35) by the energy method, we use the following identities:
We multiply the scalar explicit scheme (2.35) by the vector , we obtain
or, after substituting here the identities (2.46),
Due to the self-adjointness of the operator A and the commutativity of the scalar product, the last two terms cancel, after which we obtain the following energy identity:
If the operator
then from (2.47) we obtain the following energy inequality:
from which the maximum principle immediately follows
which is a sufficient condition for stability.
If we now calculate any norm from inequality (2.48), for example, a norm that is equal to the maximum modulo eigenvalue of the operator A, then using expression (2.45) we obtain
Thus, the Courant stability condition (2.33) of the explicit finite-difference scheme for the heat equation, derived using the method of harmonic analysis, is also a sufficient condition.
Let us now study the stability of the implicit finite-difference scheme (2.36) by the energy method, for which we add to identities (2.46) the identity
(2.51)
We multiply the scheme (2.36) scalarly by the vector , we obtain
Thus, for the implicit scheme, the energy identity has the form
Here the first term is always positive definite, so the energy inequality has the form
whence follows the maximum principle
Thus, the implicit scheme (2.36) is unconditionally stable, since the operator is always positive definite.
When constructing grid equations, two approaches are possible: 1) consideration of the GPZ within the framework of a formal problem of mathematical physics and, as a consequence, the use of methods and algorithms developed in other scientific areas operating with similar systems of equations for solving; 2) formulation of a numerical problem based on concepts and regularities arising directly from the physical consideration of the process under study.
The first approach is more universal, it allows one to draw a direct analogy between physical processes of different nature (for example, diffusion, thermal "and filtration), if they are described by isomorphic systems of equations, which at first served as an incentive for the widespread introduction of HFFM. However, later it became it is clear that the isomorphism of differential equations is not a sufficient condition for the effective use of the same numerical method in solving applied problems of different physical nature Insufficient consideration of the physics of the process led, for example, to the construction of numerical schemes based on the simplest non-conservative difference approximations (the Lax method) , which for filtration problems in inhomogeneous media entailed a significant loss of accuracy or instability of the solution.
In contrast, the second approach - from the physical formulation of the problem, bypassing the differential equation, directly to the numerical scheme - is based on a quite obvious requirement: on the difference grid, both -balance relations (scheme conservatism) and the main patterns of groundwater movement should be approximated. Note, however, that the opposite
the statement of these two approaches is rather conditional: the balance relations can also be satisfied on the basis of differential equations. However, this requirement is often not satisfied, since finite-difference approximations of differential equations formally do not require this. As an example of the second approach to the construction of difference schemes that can be implemented on an AVM, one can cite the Liebmann method, which, due to its accessibility and physical clarity, has found the widest application in the practice of solving the GFZ. This approach was first implemented by G. N. Kamensky, who proposed using difference methods to solve filtering problems.
In essence, the constructions mentioned are special cases of the integro-interpolation method, which allows you to most fully link the physical and mathematical formulations of a specific numerical problem and obtain, on this basis, a priori estimates of stability and convergence, taking into account not only the type of boundary value problem, but also its specific features. This approach will be considered below using the example of constructing a difference scheme for a planned problem of geofiltration in a heterogeneous reservoir.
Under the difference scheme, we mean a set of difference equations that approximately describe the filtration processes and additional conditions (boundary and initial) that characterize the behavior of the desired head function H(x, y, f) on the internal and external boundaries, as well as its distribution within the filtration area at the initial time. To construct difference equations describing the filtration process, according to the integro-interpolation method, the following is required:
1) replace the region of continuous variation of the argument with a discrete one: WA=W(xh y(, tn), where W is the continuous region of variation of the argument; x, y are coordinates; t is time; i, j, n are the numbers of points of the discrete domain fV& ;
2) write down for the constructed area Wg the balance identities relating the changes in the flow rate of the filtration flow within the element with area Pu) = &x(Ay) with the intensity of the change in capacitive reserves in it, as well as with the flow rates of additional sources (sinks) attributed to this element ( integral side of the method);
3) to express all the components of the balance identity in terms of the pressure values at the nodal points of the area, the parameters of the modeled system and the intervals of the space-time breakdown of the difference grid (the interpolation side of the method).
Currently, when solving problems of mathematical physics, the following main types of grids are used: 1) rectangular, uniform and non-uniform; 2) triangular and polygonal, uniform and uneven; 3) orthogonal curvilinear. The advantages of grids of the second and third types are a more accurate approximation of the outer and inner boundaries with a complex
geometry, as well as the possibility of greater detail of the structure of the filtration flow within individual subareas. However, in this variant, the construction of the grid is a rather complicated and non-universal procedure. It is also important that the use of such grids often does not lead to an increase in the real accuracy of modeling at all due to the lack of hydrogeological information necessary for their construction.
Proceeding from this, it seems inexpedient to use grids of the second and third types in solving the GPZ, especially since it is difficult to estimate the practical error of the solution, which depends on the accuracy of the grid construction (the so-called “positional” errors). As practice shows, increasing the accuracy of constructing such grids requires a significant complication of the computational algorithm, which makes the programs created on their basis even less suitable for mass use. Finally, it is important that the use of the “fictitious regions” method when approximating curvilinear boundaries on a rectangular grid makes it possible to construct numerical algorithms that are almost similar in accuracy to those used in the FEM.
Meanwhile, to build grids of the first type, a minimum of information is required, which greatly simplifies communication with the program. Now two main modifications of such grids are used, which differ in the position of the nodal points relative to the faces of the calculated blocks: in the first version, the nodal points coincide with the centers of the blocks (block mesh), and in the second - with the intersection points of the faces (nodal mesh). Block grids have undoubted advantages over nodal grids when using homogeneous difference schemes based on the idea of through counting (see Sections 2.3 and 2.4); the use of nodal grids is preferable when the difference scheme is not regular, b. also in solving mass transfer problems requiring explicit identification of singularities. Note that when referring to the method of "fictitious regions" for approximation on the boundaries of conditions of the first or second kind, these two modifications turn out to be equivalent here. _
So, to construct a difference grid y,), approximate
which forms a continuous filtration region W\x, y), the simulated field is covered with a non-uniform rectangular grid with steps along the OX and 07 axes equal to Axf and Ayj, respectively (Fig. 1).
Consider the balance of flow rates of the filtration flow within an elementary block with coordinates M(xl+1, yj+1)1 (see Fig. 1) at time tL
1 The choice of point indices /+1 and j+1 (instead of i and j) is due to the need to approximate the boundary conditions using fictitious boundary blocks.
Rice. 1. Scheme of breakdown of the filtration area when constructing a balance identity:
1 anchor point; 2 - axial line of the calculation block; J - combed block border; 4 - the boundary of the block, relative to which the balance identity is built. The direction of the flow is given by the spherules
where L is an index that determines the choice of tL for a particular balance element; n = 1, 2, ..., K\ Ar -ve
the name of the time interval; K is the number of time steps. In general, expenses
filtration flow through the faces of the block can be recorded for different Time Points belonging to the selected interval. By virtue of the law of conservation of the filtration flow mass, the total increment of these costs is equal to the total intensity of water inflow into the block under consideration due to the depletion of capacitive reserves QE(x, y, t); overflows O.P(x, Y, t); inflow from imperfect reservoirs QB(x, y, t) or Qu(x, y, t) (if their sizes are commensurate with the grid spacing) or QP(x, y, t) [if their sizes are much smaller than it]; infiltration recharge QM(x, y, z) and well operation QC(x, y, /).
Taking into account all the balance elements, and also assuming at first that t1 = t2 = ti = t4 = t, for the considered hydrodynamic scheme, we can write the balance identity, which is valid for any fixed moment t:
Q A 1/2, j+1 - 3/2, j+ 1 + 6 i3+l,j+ 1/2 ~Qi\ 1 ,j+ 3/2 -
The nonlinearity of the process associated with the dependence of conductivity on the pressure (depth) of the flow can often be taken into account by linearization, carried out by introducing a constant calculated conductivity over a given time interval.
The value of the discharge QE, due to the change in capacitive reserves in the element of the aquifer with area F, is determined by the expression
where (1° and u* are the gravitational and elastic capacities (water loss) of the reservoir; H is the level of the free surface; R is the average head in the reservoir.
When specifying the elastic capacity, it is assumed that rock deformations are conditionally instantaneous elastic in nature and are linearly related to pressure (head) changes. In reality, the rock is a heterogeneous system and its stress-strain state within each element can be significantly inhomogeneous. To take this circumstance into account, a heterogeneous block model or a scheme of a medium with a double capacity is used, in which the rock is represented as consisting of a quasi-homogeneous system of weakly permeable blocks, uniformly separated It is assumed that only the flow in the channels directly reacts to changes in the hydrodynamic situation, and the reaction of blocks slows down due to their resistance.Then, the flow rate QE includes the capacities of fractures and blocks, multiplied by the rate of change of pressure in the corresponding elements of the heterogeneous medium.
As is known, the gravitational capacity can change for two reasons: due to the lithological variability of the cover layers, within which the free surface passes, and due to the influence of the capillary zone. With a significant role of the second mechanism, especially in problems with areal supply, it is recommended to switch to the joint flow model in a saturated-unsaturated medium.
A significant role in the formation of groundwater flows belongs to the areal recharge with a flow rate QH=eF (where c is the intensity of the areal recharge, reflecting the total effect of infiltration and evaporation and depending on the depth of the groundwater level). When solving predictive regional problems for reclaimed territories, the QH value includes the module of drainage flow u, which varies depending on the pressure
The flow rate QII through the separating reservoir with the crossflow coefficient AP=% = kr1tr (where kr is the filtration coefficient; tr is the reservoir thickness) is determined, according to the simplest and most common design scheme of hard overflow, by the expression
where H+ is the head in the adjacent aquifer from which the overflow occurs. Note that the value of the flow coefficient can vary significantly over the area due to manifestations of the structural and lithological variability of rocks. In a more general flow model, the elastic regime in separating layers should be taken into account, the manifestations of which, apparently, can be significantly complicated by the influence of the heterogeneity of low-permeable rocks, as well as violations of the linear filtration law.
We limited ourselves to the six-term form of writing the right-hand side of identity (2.1), which takes into account only infiltration recharge, overflow from neighboring layers and reservoirs at constant pressures in them and a hard filtration regime in the separating layer, as well as the presence of wells. This is due to three reasons: 1) a wide range of problems can be reduced to the considered geofiltration scheme, which makes it sufficiently universal for practical use; 2) the analysis of the possibilities of efficient solution of the systems of equations obtained on the basis of the integral identity (2.1) shows that any additional increase in its components inevitably leads to a sharp increase in the cost of the volume of the operational filling device (RAM) and the counting time (computer resources). The decrease in the universality of identity (2.1), i.e., the exclusion of individual terms from its right-hand side, practically does not lead to a decrease in the used computer resources; 3) the accepted form of representation of additional sources of supply allows direct generalization to the case of multilayer aquifers.
Thus, already at the first stage of developing the difference model, the authors tried to avoid excessive universalism and limit the range of problems to be solved. This aspect of computational schematization (CSN) will further allow us to build an efficient elementary computational module, which can then be generalized to multilayer systems as necessary, using additional computer resources for this.
1. Examples of difference approximations.
Various methods for the approximate replacement of one-dimensional differential equations by difference equations have been studied previously. We recall examples of difference approximations and introduce the necessary notation. We will consider a uniform grid with a step h, i.e. set of points
wh=(x i =ih, i=0,± 1, ± 2,…}.
Let u(x) is a sufficiently smooth function defined on the segment . DenoteDifference relations
are called respectively the right, left and central difference derivatives of the function u(x) at the point x i, i.e. at a fixed x i and for h®0 (and thus for i®¥) the limit of these relations is u'(x i). Expanding according to the Taylor formula, we obtain
u x,i – u’(x i) = 0.5hu’’(x i) + O(h 2),
u x,i - u’(x i) = -0.5hu’’(x i) + O(h 2),
u x,i - u'(x i) = O(h 2),
This shows that the left and right difference derivatives approximate u'(x) with first order h, and the central difference derivative is of the second order. It is easy to show that the second difference derivative
approximates u''(x i) with second order h, and the decomposition
Consider the differential expression
with variable ratio k(x). We replace expression (1) with the difference relation
where a=a(x) is a function defined on the grid w h . Find the conditions that the function must satisfy a(x) in order for the relationship (au x) x,i approximated (ku')' at the point x i with second order h. Substituting into (2) expansions
where u i' = u'(x i), we get
On the other side, Lu = (ku')' = ku'' + k'u',
those.
From here it is clear that L h u–Lu = O(h 2) if conditions are met
Conditions (3) are called sufficient conditions of the second order of approximation. When deriving them, it was assumed that the function u(x) has a continuous fourth derivative and k(x) is a differentiable function. It is easy to show that conditions (3) are satisfied, for example, by the following functions:
Note that if we put a i = k(x i), then we get only the first order of approximation.
As the next example, consider the difference approximation of the Laplace operator
Let us introduce a rectangular grid on the plane (x 1 , x 2 ) with a step h 1 in the x 1 direction and with a step h 2 in the x 2 direction, i.e. set of points
wh= ((x i 1 , x j 2) | x i 1 = ih 1 , x j 2 = jh 2 ; i, j = 0,± 1, ± 2,…},and denote
It follows from the previous considerations that the difference expression
approximates the differential expression (4) with the second order, i.e. L h u ij – Lu(x i 1 , x j 2) = O(h 2 1) + O(h 2 2). Moreover, for functions u(x 1 , x 2) with continuous sixth derivatives, the expansion
Difference expression (5) is called five-point difference Laplace operator, since it contains the values of the function u(x 1 , x 2) at five grid points, namely at the points (x 1 i , x 2 j), (x 1 i ± 1 , x 2 j), (x 1 i , x 2 j ± 1). The indicated set of points is called the difference operator template. Difference approximations of the Laplace operator are also possible on templates containing a greater number of points.
2. Study of approximation and convergence
2.1. Approximation of a differential equation. Previously, we considered the boundary value problem
(k(x) u’(x))’ – q(x) u(x) + f(x) = 0, 0< x < l, (1)
– k(0) u’(0) +bu(0) =m 1 , u(l) =m 2 , (2)
k(x)³ c1 > 0,b³ 0,
for which the difference scheme was constructed by the integro-interpolation method
Denote by Lu(x) the left side of equation (1) and through L h y i is the left side of equation (3), i.e.
Let u(x) is a fairly smooth function and u(x i) is its value at the point x i grids
wh= (x i = ih, i = 0, 1, …,N, hN = l) (7)
They say that difference operator L h approximates the differential operator L at the point x=x i if the difference L hui– L hu(x i) tends to zero as h®0. In this case, it is also said that the difference equation (3) approximates the differential equation (1).
To establish the presence of an approximation, it suffices to expand by the Taylor formula at the point x=x i values ui ± 1 = u(x i± h), included in the difference expression L hui. Much of this work was done in the previous chapter, where it was shown that under the conditions
(8)the relation
We write the differential problem in operator form
LU = f, where L is one of the differential operators
U(x,y)- the desired function that satisfies the differential problem; f- input data (i.e., initial and boundary conditions, right-hand sides, etc.). The operator form describes the differential problem at the grid nodes, and the operator form describes the finite difference scheme on the exact solution U(x,t), i.e. in the finite difference scheme, the exact (unknown) values of the required function are substituted for the grid values of the grid function. The operator form of the finite difference scheme has the form
Let us introduce the norm of the grid function, for example, using the expression
Definition. The finite difference scheme approximates the differential problem on the exact solution, if any norm of the difference? tends to zero at
Definition. The finite difference scheme approximates the differential problem on the exact solution with the order p in time and order q in a space variable, if any norm of the difference satisfies the equality
Thus, if a finite-difference scheme approximates a differential problem, then we are talking about the closeness of the differential and finite-difference operators at the grid nodes.
Sustainability
Let the input data in the finite-difference scheme be perturbed and become. Then the grid function will also receive a perturbation and become.
Definition. The finite difference scheme is stable with respect to the input data if there is such a bounded constant K, independent of grid characteristics and input data, which satisfies inequality (4)
Thus, the concept of stability is interpreted as follows: a finite-difference scheme is stable if, for small perturbations of the input data (initial-boundary conditions and right-hand sides), the finite-difference scheme ensures small perturbations of the grid function, i.e. the solution using a finite difference scheme is under the control of the input data.
If the input data includes only initial conditions, or only boundary conditions, or only right-hand sides, then one speaks of stability, respectively, in terms of initial conditions, in terms of boundary conditions, or in terms of right-hand sides.
Definition. The finite difference scheme is absolutely (unconditionally) stable if inequality (4) is satisfied for any values of the grid characteristics φ and h, i.e. no restrictions are imposed on the grid steps.
Definition. The finite-difference scheme is conditionally stable if inequality (4) is satisfied for the grid characteristics φ and h, which are subject to certain restrictions.
Explicit schemas
Consider the first initial-boundary value problem for the wave equation. On the space-time grid, we will approximate the differential equation of one of the following finite-difference schemes:
with the template in figure 1 and
with template in figure 2
In this case, scheme (1) is explicit. With its help, the solution is determined immediately, since the values of the grid functions on the lower time layers must be known. According to the template for this scheme, the order of approximation is two, both in the spatial and in the temporal variable. In this case, the explicit finite-difference scheme (1) for the wave equation is conditionally stable with the condition imposed on the grid characteristics f and h.
Implicit Schemas
Scheme (2) is an implicit scheme and is absolutely stable. It can be reduced to a SLAE with a tridiagonal matrix solved by the sweep method.
To determine, you can use the simplest approximation of the second initial condition.
From where we obtain the following expression for the desired values:
Implicit schemas are usually stable.
