Linear and Nonlinear Splitting Schemes Conserving Total Energy and
Mass in the Shallow Water Model
YURI N. SKIBA
Instituto de Ciencias de la Atmósfera y Cambio Climático,
Universidad Nacional Autónoma de México ,
Av. Universidad #3000, CU/UNAM, Coyoacán, Mexico City, 04510,
MEXICO
Abstract: - Two linear and one nonlinear implicit unconditionally stable finite-difference schemes of the second-
order approximation in all variables are given for a shallow-water model including the rotation and topography
of the earth. The schemes are based on splitting the model equation into two one-dimensional subsystems. Each
of the subsystems conserves the mass and total energy in both differential and discrete (in time and space)
forms. One of the linear schemes contains a smoothing procedure not violating the conservation laws and
suppressing spurious oscillations caused by the application of central-difference approximations of spatial
derivatives. The unique solvability of the linear schemes and convergence of iterations used to find their
solutions are proved.
Key-Words: - Shallow-water model, splitting method, linear and nonlinear numerical schemes, conservative
laws, unique solvability of the linear schemes, convergence of the iterative process.
1 Introduction
The shallow-water equations describe a thin layer of
fluid of constant density in hydrostatic balance,
bounded from below by the bottom topography and
from above by a free surface. They exhibit a rich
variety of features, because they have various
conservation laws. Shallow-water equations can be
used to model Rossby and Kelvin waves in the
atmosphere, rivers, lakes, and oceans as well as
gravity waves in a smaller domain (e.g. surface
waves in a bath), [1], [2], [3], [4], [5], [6].
The equations are derived from depth-integrating
the Navier-Stokes equations, in the case where the
horizontal length scale is much greater than the
vertical length scale. Under this condition,
conservation of mass implies that the vertical
velocity scale of the fluid is small compared to the
horizontal velocity scale. It can be shown from the
momentum equation that vertical pressure gradients
are nearly hydrostatic, and that horizontal pressure
gradients are due to the displacement of the pressure
surface, implying that the horizontal velocity field is
constant throughout the depth of the fluid. Vertically
integrating allows the vertical velocity to be
removed from the equations. The shallow-water
equations are thus derived.
The classic system of the shallow-water
equations taking account of the Coriolis force and
topography can be written as
+
+
0
u u u h
u v fv g
t x y x
v v v h
u v fu g
t x y y
HuH vH
t x y
(1)
In, [7], [8], [9], where
( , , )u x y t
and
( , , )v x y t
are
the components of the velocity vector,
0
f f y

is
the Coriolis parameter,
( , , )h x y t
is the free surface
height,
( , , ) ( , , ) ( , )
T
H x y t h x y t h x y
(2)
and
( , )
T
h x y
is the topography height (Figure 1).
System (1) is considered with initial conditions
0
( , ,0) ( , )u x y u x y
,
,
0
( , ,0) ( , )H x y H x y
(3)
in a channel
( , ) : 0 ; 0D x y x X y Y
with periodic conditions in x:
Received: December 11, 2022. Revised: October 18, 2023. Accepted: November 13, 2023. Published: December 14, 2023.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2023.18.18
Yuri N. Skiba
E-ISSN: 2224-347X
193
Volume 18, 2023
( , , ) (0, , )u X y t u y t
,
( , , ) (0, , )v X y t v y t
( , , ) (0, , )H X y t H y t
(4)
In addition, at the boundaries the y-component of the
velocity vector must vanish:
( ,0, ) ( , , ) 0v x t v x Y t
(5)
Fig. 1: Shallow-water model.
Numerical simulation of the problem of
atmospheric dynamics, [9], [10], [11], [12], [13].
[14], shows that the best difference schemes are
those that retain the most important properties of the
original differential model. Full conservatism is one
of such properties, [15]. As is known, [8], system
(1)-(5) has several conservation laws, in particular,
it preserves the mass
D
M h dD
and the total
(kinetic plus potential) energy
( ) ( ) ( )E t K t P t
:
( ) 0
dMt
dt
(6)
( ) 0
dEt
dt
(7)
where
22
1
2
D
K u v H dD




,
22
1
2T
D
P g h h dD




are the kinetic and potential energies of the system,
respectively, and
dD dxdy
. Note that
2 2 2
T
DD
dd
g h h dD g h dD
dt dt


[10], proposed a spatial discretization of the
shallow-water equations that precisely conserves the
energy and potential enstrophy in semi-discrete
(discrete in space) system. Their algorithm is
important because the simultaneous conservation of
energy and potential enstrophy is known to prevent
the spurious cascade of energy to high
wavenumbers.
However, the numerical schemes proposed in,
[9], [10], [11], [16], conserve the energy and
enstrophy only when the time derivatives are taken
in the differential form. The discretization of these
derivatives leads to the violation of these laws in a
completely discrete system. In this paper, three
implicit schemes are proposed that ensure the
conservation of mass and energy in a completely
discrete (in space and time) system (1)-(5). One of
the schemes is nonlinear, and the other two are
linear. Since the integral in (7) can be considered as
the square of the norm of the solution, the
constructed schemes are unconditionally stable. The
linear schemes have a second-order approximation
in all spatial and temporal variables.
It is well known that implementing multi-
dimensional implicit schemes requires a large
amount of computation. To simplify the calculations
at each small-time interval, a geometric division of
the original two-dimensional problem into two one-
dimensional ones was used. Every split differential
system satisfies laws (6) and (7). The corresponding
discrete one-dimensional systems and therefore the
entire finite-difference scheme also conserve these
two laws in discrete form.
2 Modified Shallow-Water Equation
Assuming that in the domain D the height of the
atmospheric layer
( , , )H x y t
is always positive, we
introduce new dependent variables, [14].
ZH
,
U Zu
,
V Zv
(8)
Then according to the formulas
U u Z
Zu
t t t

and
V v Z
Zv
t t t

system (1) can be rewritten as
11
+ ( ) ( )
22
U U U
uU u vU v
t x x y y
h
fV gZ x
,
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2023.18.18
Yuri N. Skiba
E-ISSN: 2224-347X
194
Volume 18, 2023
11
+ ( ) ( )
22
V V V
uV u vV v
t x x y y
h
fU gZ y
,
0
HZU ZV
t x y
(9)
Multiplying the equations of system (9) by U, V,
and gh, respectively, integrating the obtaining
equations over D, and summing the results, we
arrive at the energy conservation law (7)
2 2 2
0
22
D
d U V g h
H dD
dt







(10)
Law (6) is also valid for system (9). The main
advantage of system (9) over system (1) is “the
divergent form” of its equations. Indeed, for
example, the term
() U
uU u
xx


(11)
when multiplied by U, takes the divergent form,
which after integration over x and using periodic
conditions (4) gives
2
( ) 0
D
uU dD
x
(12)
The same is true for the remaining terms of the
system (9) enclosed in brackets. In addition, the sum
of the terms
( / )gZ h x
and
/ZU x
previously
multiplied by U and gh, respectively, also leads to
the divergent form
( ) /gZUh x
. To explain the
usefulness of transforming system (1) to the
divergent form (9), consider the central-difference
approximation
1 1 1 1 1 1
22
i i i i i i
i
u U u U U U
u
xx


(13)
of the divergent form (11), where the fixed index j
of the grid functions
( , )
ij i j
u u x y
and
( , )
ij i j
U U x y
is omitted for simplicity. Then it is
easy to see that multiplying (13) by
i
U
with
subsequent summation over all internal nodes of the
grid from
1i
to
iI
leads to
1 1 1
1
1()
2
I
i i i i i i
iu UU uU U
x
1 1 1
1
1( ) 0
2
I
i i i i i i
iu U U uU U
x
due to periodic conditions
1 1 0 1 1 0
, , ,
I I I I
u u u u U U U U

(14)
approximating (4). Thus, the finite-difference form
(13) preserves the property (12) of the divergent
form (11). This property is of great importance for
constructing a scheme that conserves the energy of
the system.
3 Splitting of the Model Operator
As mentioned above, our goal is to construct a
numerical scheme that preserves laws (6) and (7) in
the discrete form (both in time and space) and has a
second-order approximation in all independent
variables. For this reason, only implicit schemes
should be used. However, it is known that in the
case of a multi-dimensional problem the use of
implicit schemes leads to numerical algorithms, the
implementation of which is too expensive.
Therefore, to construct an accessible numerical
algorithm, we use a geometric splitting of the
original two-dimensional system into two one-
dimensional subsystems that admit a simple
solution. To illustrate the main idea of the splitting
method, [17], [18], [19], [20], consider the linear
homogeneous problem:
1 2 1
0 , ( )A A t g
t

(15)
in the domain D and on a sufficiently small-time
interval
12
( , )tt
. Here
( , , )t x y
is the solution with
the initial value of
( , )g x y
, and the operators
1
A
and
2
A
are assumed to be positive semidefinite:
0 ( 1,2)
i
D
A dD i
(16)
According to the splitting method, the first split
subsystem
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2023.18.18
Yuri N. Skiba
E-ISSN: 2224-347X
195
Volume 18, 2023
1
1 1 1 1
0 , ( )A t g
t

(17)
is solved within
12
( , )tt
, and its final solution
12
()t
is used as the initial condition for solving the second
subsystem
2
2 2 2 1 1 2
0 , ( ) ( )A t t
t

(18)
on the same interval
12
( , )tt
. Then the final solution
22
()t
of subsystem (18) will approximate the
solution
2
()t
of the original problem (15):
2 2 2
( ) ( )tt
. In particular, if
1
A
and
2
A
are the
one-dimensional operators (in the x- and y-
directions, respectively), then the two-dimensional
problem (15) reduces to the sequential solution of
one-dimensional problems (17) and (18), which
greatly simplifies the calculations. The smaller the
interval
12
( , )tt
, the closer
22
()t
to
2
()t
. Thus,
by choosing a sufficiently small interval
12
( , )tt
, the
error of the splitting procedure can be made
arbitrarily small. The same approach can also be
applied to nonlinear problems.
Let us split system (9) on a small-time interval
12
( , )tt
in accordance with algorithms (17), (18). The
subsystem
11
+ ( )
22
U U h
uU u fV gZ
t x x x



(19)
11
+ ( ) 0
22
VV
uV u fU
t x x



(20)
0
HZU
tx



(21)
is first solved with the initial conditions
1
( , , )U x y t
,
1
( , , )V x y t
and
1
( , , )H x y t
. Then its solution
2
( , , )U x y t
,
2
( , , )V x y t
and
2
( , , )H x y t
at
2
tt
is
used as the initial condition at
1
tt
to solve the
subsystem
11
+ ( ) 0
22
UU
vU v fV
t y y



(22)
11
+ ( )
22
V V h
vV v fU gZ
t y y y



(23)
0
HZV
ty



(24)
on the same time interval
12
( , )tt
. The solution of this
subsystem obtained at time
2
tt
, approximates the
solution of the original non-split system (9).
Obviously, the mass conservation law (6) is satisfied
for each of the split systems (19)-(20) and (22)-(24).
Now let’s multiply the equations (19) and (22) by U,
(20) and (23) by V, and (21) and (24) by gh, sum the
results, and integrate the final equation over D.
Using the divergent form of terms and boundary
conditions (4) and (5), it is easy to show that the
conservation law (7) is also valid for each of split
systems (19)-(21) and (22)-(24).
4 Linear and Nonlinear Implicit
Schemes Conserving the Mass and
Total Energy
Due to a certain symmetry in the structure of
subsystems (19)-(21) and (22)-(24), the numerical
scheme constructed for the first subsystem can be
easily modified and applied to the second. Let us
introduce net functions at the grid nodes
( , , )
i j n
x y t
,
where
( , )
ij
xy
is the internal point of the domain D
(
1iI
,
1jJ
),
n
t
is a discrete moment in
time (
1n
), and denote
1nn
tt

,
1ii
x x x
,
1jj
y y y
,
()
jj
f f y
,
( , , )
n
ij i j n
R R x y t
where
R
can be one of the functions
, , , , ,u v h Z U V
or
H
, and
,
x
and
y
are small steps of regular
grids. Moreover, we will use the notation
1
1()
2
nn
ij ij ij
R R R

(25)
At each interval
1
( , )
nn
tt
and for each fixed
j
y
,
system (19)-(21) is solved using the implicit scheme
1
1 1 1 1 1 1
1[]
2 2 2
nn
i i i i i i i i
i
U U u U u U U U
u
xx


11
1
22
ii
j i i
hh
f V gZ x

(26)
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2023.18.18
Yuri N. Skiba
E-ISSN: 2224-347X
196
Volume 18, 2023
1
1 1 1 1 1 1
1[]
2 2 2
nn
i i i i i i i i
i
V V u V u V V V
u
xx


10
2ji
fU
(27)
1
1 1 1 1 0
2
nn
i i i i i i
H H Z U Z U
x


(28)
with periodic boundary conditions. The functions U,
V, and h are defined here by formula (25), and the
fixed index j is omitted for simplicity of notation.
Further, for each fixed
i
x
, the net functions
11
,
nn
ij ij
UV

and
1n
ij
H
, obtained from (26)-(28), are
used as initial conditions
n
ij
U
,
n
ij
V
and
n
ij
H
to solve
system (22)-(24) on the same time interval
1
( , )
nn
tt
according to the implicit scheme
1
1 1 1 1 1 1
1[]
2 2 2
nn
j j j j j j j j
j
U U v U v U U U
v
yy


10
2jj
fV
(29)
1
1 1 1 1 1 1
1[]
2 2 2
nn
j j j j j j j j
j
V V v V v V V V
v
yy


11
1
22
jj
j j j
hh
f U gZ y

(30)
1
1 1 1 1 0
2
nn
j j j j j j
H H Z U Z U
y


(31)
(the fixed index i is also omitted here) with the
boundary conditions
0 1 1
11
( ) 0 , ( ) 0
22
JJ
V V V V
0 1 1
11
( ) 0 , ( ) 0
22
JJ
v v v v
0 1 0 1 1 1
, , ,
J J J J
Z Z h h Z Z h h

Then this procedure is repeated on the interval
12
( , )
nn
tt

, while the functions
1n
ij
U
,
1n
ij
V
and
1n
ij
H
found on the previous interval are taken as
initial conditions when solving system (26)-(28).
Note that, depending on the choice of the
functions
,
ij ij
uv
and
ij
Z
, at least three different
schemes can be constructed:
1) If
,,
ij ij ij ij ij ij
u u v v Z Z
(32)
where
and,
ij ij ij
u v Z
are defined by (25), both
schemes (26)-(28) and (29)-(31) are nonlinear.
2) If
,,
n n n
ij ij ij ij ij ij
u u v v Z Z
(33)
where
, , and
n n n
ij ij ij
u v Z
are the solutions obtained on
the previous time interval
1
( , )
nn
tt
, schemes (26)-
(28) and (29)-(31) are linear.
3) If
,,
( ) , ( )
n n n n
ij km km ij km km
k m k m
u ij u v ij v



,
,
()
nn
ij km km
km
Z ij Z
(34)
where
, , ,
( ) 1, ( ) 1, ( ) 1
n n n
km km km
k m k m k m
ij ij ij
the schemes (26)-(28) and (29)-(31) are also linear.
Here
and, ,
n n n
ij ij ij
u v Z
are the same as in (33), and
( ) , ( )
nn
km km
ij ij

and
()
n
km ij
are certain
coefficients of the interpolation used to smooth the
functions
,
ij ij
uv
and
ij
Z
. The latter procedure can
be effective for suppressing spurious oscillations
caused by the use of central difference
approximations.
Schemes (26)-(28) and (29)-(31) conserve mass
and total energy regardless of the choice of (32),
(33), or (34) for the functions
,
ij ij
uv
and
ij
Z
.
Indeed, if we multiply (28) by
xy

and sum over
all internal grid nodes
( , )
ij
xy
(
1iI
,
1jJ
).
Then:
1
,,
nn
ij ij
i j i j
x y H x y H

(35)
since
1, 1, 1, 1,
[ ] 0
I j I j j j
j
Z U Z U


WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2023.18.18
Yuri N. Skiba
E-ISSN: 2224-347X
197
Volume 18, 2023
due to the periodic conditions
1, 1,I j j
ZZ
and
1, 1,I j j
UU
(see (14)). In complete analogy with
(28), equation (31) leads to (35) due to the
conditions
,0 , 1 0
i i J
VV

(see (5)). Next, if we
multiply (26), (27) and (28) by
i
x yU

,
i
x yV

and
i
x y gh

, respectively, and sum
the resulting equations over all internal grid nodes
( , )
ij
xy
of the domain D then the use of (2), (25)
and the boundary conditions leads to the
conservation of the total energy
()
nn
E E t
in the
finite-difference system
2 2 2
1 1 1 1
,
1
2
n n n n
ij ij ij
ij
E x y U V g h
2 2 2
,
1
2
n n n n
ij ij ij
ij
x y U V g h E
(36)
Relation (36) is also fulfilled for system (29)-
(31). Defining the solution norm as
1/2
nn
E
we get
1nn
This means that the difference scheme (26)-(31)
is unconditionally stable.
Note that the linear scheme (26)-(31) has the
second order of accuracy in
x
and
y
, but the
first order of accuracy in
, since the matrices
approximating the dynamic operators
1
A
and
2
A
of
the splitting algorithm (17), (18) do not commute:
1 2 2 1
A A A A
However, a second-order approximation in
can be achieved using the symmetrized algorithm,
[18], when the original 2D problem is solved on
each double interval
11
( , )
nn
tt

. Moreover, first the
scheme (26)-(28) and then the scheme (29)-(31) are
solved on
1
( , )
nn
tt
. After this, on the interval
1
( , )
nn
tt
the order changes and first the scheme
(29)-(31) is solved, and then the scheme (26)-(28).
Thus, within the framework of the symmetric
algorithm, scheme (29)-(31) can be immediately
resolved on the double interval
11
( , )
nn
tt

. Note that
when using the symmetric algorithm on the double
time interval
11
( , )
nn
tt

, conditions (33) must be
replaced by solutions obtained on the previous
double time interval
31
( , )
nn
tt

:
1 1 1
,,
n n n
ij ij ij ij ij ij
u u v v Z Z
5 Unique Solvability of the Linear
Schemes
Schemes (26)-(28) and (29)-(31) are solved using an
iterative method. Let us now show that in cases (33)
and (34) when these schemes are linear, their
solutions are unique. Indeed, each of these schemes
can be written as a vector equation
10
nn

(37)
where
1
1
2()
nn
, [21], the components of
vector
n
are values of net functions
,
nn
ij ij
UV
and
n
ij
hg
, and the time-dependent matrix
()
n
t
is
skew-symmetric due to the energy conservation law
(36):
0
T
(38)
Here the superscript T means transposition.
Therefore, there is a complete orthogonal system of
eigenvectors
k
of matrix
:
k k k

,
T
k m km
(39)
where
km
is the Kronecker symbol and the
eigenvalue
k
is pure imaginary for each k. Let’s
rewrite (37) as
Bw F
(40)
where
1n
w

is the solution to be found,
F
is a
known vector depending on
n
t
,
1
()
2
BE
and E is the unit matrix. By virtue of (38), the
matrix B is positive definite:
0
TB
.
Moreover, it has the same eigenvectors
k
as
the matrix
with eigenvalues
1
2
1
kk

:
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2023.18.18
Yuri N. Skiba
E-ISSN: 2224-347X
198
Volume 18, 2023
k k k
B
. Since every
k
is purely
imaginary,
1/2
2
2
1
4
11
kk
,
the matrix B is nonsingular, and problem (40) (and
therefore the linear scheme (26)-(31)) has a unique
solution
w
at each time step
1
( , )
nn
tt
.
6 Convergence of the Iterative Process
We now derive the condition under which the
iterative process:
1()
p p p
w w Bw F
(41)
used to solve (40) converges. Here
is the
relaxation parameter. The error vector
pp
ww

satisfies the equation:
1pp
T

(42)
where
T E B

is the transition matrix from
p
to
1p
. Using Fourier series:
()
p m m
m
p
,
() T
m m p
p
,
eigenvectors (39) and equation (42) we obtain:
( 1) (1 ) ( )
m m m
pp
[(1 ) ] ( )
2mm
p
Therefore, the Fourier coefficients
()
mp
tend to
zero as
p
if and only if:
max (1 ) 1
2m
m
or
22
1
4
2
1

(43)
where
max m
m

. Thus, the iterative method
(41) is convergent under condition (43).
7 Conclusion
One nonlinear and two linear implicit
unconditionally stable finite-difference schemes of
the second-order approximation in all variables are
given for a shallow-water model that includes the
rotation and topography of the earth. The schemes
are developed using the splitting of the model
equation into two one-dimensional subsystems. Each
subsystem conserves the mass and total energy in
differential and discrete (in time and space) forms.
One of the linear schemes contains a smoothing
procedure not violating the conservation laws and
suppressing spurious oscillations caused by the
application of central-difference approximations of
spatial derivatives. Unique solvability of the linear
schemes is shown and convergence of iterations used
to find their solutions is proved.
Acknowledgement:
This work was partially supported by the National
System of Researchers of Mexico (SNI,
CONACYT) through the grant 14539.
References:
[1] Anderson D., M. Harris, H. Hartle, D.
Nicolsky, E. Pelinovsky, A. Raz and A.
Rybkin, Run-up of long waves in piecewise
sloping U-shaped bays, Pure and Applied
Geophysics, Vol. 174, No. 8, 2017, pp. 3185-
3207.
[2] Battjes, J.A. and R.J. Labeur, Unsteady Flow
in Open Channels, Cambridge University
Press, 2017.
[3] Carrier G. F. and H. Yeh, Tsunami
propagation from a finite source, Computer
Modeling in Engineering & Sciences, Vol. 10,
No. 2, 2005, pp. 113–122.
[4] Didenkulova I. and E. Pelinovsky, Rogue
waves in nonlinear hyperbolic systems
(shallow-water framework), Nonlinearity,
Vol. 24, No. 3, 2011, pp. R1–R18.
[5] Lindborg E. and A.V. Mohanan, A two-
dimensional toy model for geophysical
turbulence, Physics of Fluids, Vol. 29, No.
11, 2017, ID 111114.
[6] Skiba Yu.N. and D.M. Filatov, Numerical
simulation of coastal flows in open multiply-
connected irregular domains. In: Simulation
and Modeling Methodologies, Technologies
and Applications. Advances in Intelligent
Systems and Computing, Springer, Vol. 256,
pp. 71-84, 2014.
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2023.18.18
Yuri N. Skiba
E-ISSN: 2224-347X
199
Volume 18, 2023
[7] Batchelor G.K., An Introduction to Fluid
Dynamics, Cambridge University Press, 1970.
[8] Pedlosky J., Geophysical Fluid Dynamics,
Springer-Verlag, 1979.
[9] Sadourny R., The dynamics of finite-
difference models of the shallow-water
equations, Journal of Atmospheric Sciences,
Vol. 32, No. 4, 1975, pp. 680-689.
[10] Arakawa A. and V.R. Lamb, A potential
enstrophy and energy conserving scheme for
the shallow-water equations, Monthly
Weather Review, Vol.109, No. 1, 1981, pp.
18-36.
[11] Salmon R., Poisson-bracket approach to the
construction of energy- and potential
enstrophy-conserving algorithms for the
shallow-water equations, Journal of
Atmospheric Sciences, Vol. 61, 2004, pp.
2016-2036.
[12] Skiba, Yu.N., Mathematical Problems of the
Dynamics of Incompressible Fluid on a
Rotating Sphere. Springer International
Publishing AG, Switzerland, 2017.
[13] Vreugdenhil C.B., Numerical Methods for
Shallow-Water Flow, Kluwer, 1994.
[14] Weiyan Tan, Shallow-Water Hydrodynamics,
Elsevier, 1992.
[15] Popov Yu.P. and A.A. Samarskii, Complete
conservative difference schemes, USSR
Computational Mathematics and
Mathematical Physics, Vol. 9, No. 4, 1969,
pp. 296-305.
[16] Ringler T.D. and D.A. Randall, A potential
enstrophy and energy conserving numerical
scheme for solution of the shallow-water
equations on a geodesic grid, Monthly
Weather Review, Vol. 130, 2002, pp. 1397-
1410.
[17] Douglas J. and H.H. Rachford, On the
numerical solution of heat conduction
problems in two and three space variables,
Transactions of the American Mathematical
Society, Vol. 82, No. 2, 1956, pp. 421-439.
[18] Marchuk G.I., Methods of Numerical
Mathematics, Springer-Verlag, 1982.
[19] Skiba Yu.N. and D. Filatov, Conservative
schemes based on the separation method, for
the numerical simulation of vortices in the
atmosphere, Interciencia, Vol. 3, No. 1, 2006,
pp. 16-21.
[20] Yanenko N. N., Method of Fractional Steps,
Springer-Verlag, 1971.
[21] Crank J. and P. Nicolson, A practical method
for numerical evaluation of solutions of partial
differential equations of the heat conduction
type, Proceedings of the Cambridge
Philosophical Society, Vol. 43, 1947, pp. 50-
67.
Contribution of Individual Authors to the
Creation of a Scientific Article (Ghostwriting
Policy)
The author contributed in the present research, at all
stages from the formulation of the problem to the
final findings and solution.
Sources of Funding for Research Presented in a
Scientific Article or Scientific Article Itself
This work was partially supported by the National
System of Researchers of Mexico (SNI,
CONACYT) through the grant 14539.
Conflict of Interest
The author has no conflicts of interest to declare.
Creative Commons Attribution License 4.0
(Attribution 4.0 International, CC BY 4.0)
This article is published under the terms of the
Creative Commons Attribution License 4.0
https://creativecommons.org/licenses/by/4.0/deed.en
_US
WSEAS TRANSACTIONS on FLUID MECHANICS
DOI: 10.37394/232013.2023.18.18
Yuri N. Skiba
E-ISSN: 2224-347X
200
Volume 18, 2023