- Email: [email protected]
Violation of conservation laws when solving difference equations by iteration methods
131
Solving difference eqlcarions b_~,iteration methods 8. DAUTOV. R. Z.. and LAPIN. A. V.. Mesh schemes of arbitrary order of accuracy for quasi-linear elliptic equations. Ix I‘XOI: Ser. marem, No. 10 (209), 24-37, 1980. 9. DAUTOV, R. Z., On superconvergence of mesh schemes for solving the third boundary value problem fol quasi-linear elliptic equations, in: ,Vumerical merhods of the mechanics of continuous media (Chisl. meted) mekhan. sploshnoi sredyj. Vol. 2. No. 1.62-81. Nauka. Novosibirsk, 1980.
10. MAKAROV, V. L., and SAMARSKII. A. A.. Application of exact difference schemes to the estimation of the convergence of the method of straight lines. Zh. v,$hisl Mof. mu?. Filir, 20, No. 2. 371-387. 1980.
11. ODEN. J. T., and REDDY. J. H.. An Introduction to the Mathematical Theor), of Finite Elements, H’ileyInterscience. New York. 1976. 12. TIKHONOV. A. N.. and SAMARSKII, A. A.. On homogeneous difference schemes of high order of accurac!. DokL Akad. Nouk SSSR, 131, No. 3,514-517. 1960. 13. TIKHONOV, A. N.. and SAMARSKII. A. A., Homogeneous difference schemes on non-uniform meshes, Zk v_i&isL Mot. mat. Fii., 2, No. 5, 812-832, 1962. 14. TIKHONOV. A. N.. and SAMARSKII, A. A., On homogeneous difference schemes of high order of accurac! on non-uniform meshes, Zh v~chjsl. Mar. mat. Fiz., 3, No. 1, 99-108, 1963. 15. KORNEEV. V. G.. Mesh operators, energ)-wise equivalent to the operators generated by piecewise Hermitian extensions, Zh y:chisL Mot. mat. Fii., 19, NO. 2,402-416, 1979. 16. KORNEEV, V. G. Finite element method: convergence to piecewise polynomial interpolation of the e\act solurion and convergence on certain sets of points, Part I. Dep. at VINITI, No. 4754-81, DEP? 1981; Part II, No. 4755-81 DEP; Part III, No. 4756-81 DEP. 17. KORNEEV. V. G., On exact mesh schemes. Zh v.TchisLMa?. mot. Fiz., 22. No. 3. 646-654,
1982.
18. BRA\lBLE. J. H., and HILBERT. S. R., Estimation of linear functionals on Sobole\ spaces with apptication to Fourier transforms and \pline interpolation, SIAM, J. IVumer. Anoljxsis. 7. 113- 124. 1970. 19. KORNEEV. V. G.. Schemes oj the finite-element method of high orders of occurac?, (Skhem!, metoda konechn! kh elementov v!,sokikh por!,adkov tochnosti). Izd-vo Leningr. un-ta. Leningrad. 1977. 20. SEARLE. F.. Method offinite
elemenrs for ellipric problems (Russian translation. Mir. \los
1980).
21. KORNEEV. V. G.. On the construction of variational difference schemes of high orders 01‘ accurac!, I’.csm Leningr. un-ta, 4, No. 19. 28-40. 1970.
U.S.S.R. Comput. Moths. Moth Pl7~s. Vol. 22. No. 5. pp. 131-139. Printed in Great Britain
0041-5553/82507.50+.00 Gl964. Perfamon Press Ltd.
1982.
VIOLATION OF CONSERVATION LAWS WHEN SOLVING DIFFERENCE EQUATIONS BY ITERATION METHODS* M. Yu. SHASHKOV Moscow (Received
75 Sepfember 1980; revised 16 Februac,
IT IS SHOWN that the use of iteration to conservative
difference
free of this disadvantage
methods for solving difference
systems may result in a considerable is introduced.
Results of calculations
lZil @chisL Mot. mat. Fi:., 22, 5, 1149-1156.
1982.
198 I )
equations
energy imbalance. are given.
which correspond A class of methods
M Yu Shashkov
132
Introduction When solving difference equations
it is advisable for several reasons to use iteration methods
111. When constructing
a difference scheme one endeavours
to represent correctly, the basic
properties of the physical process described by the differential basic conservation
laws in difference form is an important
Suppose a conservative scheme is constructed difference equations obtained
are solved by an iteration
equations.
requirement
The satisfaction
imposed on the scheme [2].
by some method, and the corresponding
method. The mesh function
(or several functions)
satisfies the difference scheme with a certain accuracy. Consequently,
is conservative,
the solution obtained may not satisfy the conservation
much the magnitude
of the
of the resulting energy. imbalances
even if the scheme
laws. The question of how
depends on the accuracy of the iteration
process. and how much on the parameters of the difference scheme itself (e.g. on the time and space steps). therefore becomes important. imbalance
in multidimensional
It is particularly
important
to know how to assess the
problems where highly accurate calculations
are impossible.
A way out may be to set up iteration processes for which an imbalance does not occur for any itetation
accuracy. Another possibility
enable any iteration
is to form difference schemes of a special kind, which
process in which the imbalance
These possibilities
does not occur to be employed.
are discussed in the present paper using an example of the second boundary
value problem for the one-dimensional
equation
of heat conduction.
Below it is shown that in using the non-conservative
Seidel method there is a considerable
imbalance in the amount of heat. and this imbalance may increase with time. A theoretical assessment of the magnitude of the imbalance is given, and it is compared with the actual)!,occurring value obtained when solving a model problem. The conditions for the conservation of a two-layer iteration process are obtained, and an example of a conservative process is given. A flow scheme for the equation of heat conduction which eliminates conservative properties of iteration processes is described.
1. Example of a non-conservative Consider the second boundary
the problems connected
iteration
with the non-
process
value problem for the one-dimensional
equation
of heat
conduction:
au
PU
at
dU
dl
I
-q(t). IPO
II
are given functions. where cc, 3, u. instant of time t is given by
(1.1)
OCSCI,
ax"
-=-(
au
asI*-i=zC:(t),
(1.2)
1,,,=u,(s). The amount of heat Q contained
(1.3) in the system at the
Solvingdifferewe equations by iteration methods
Q(t)=
133
ju(x,t)dz.
(1.1)
0
It follows from Eq. (I .I) and the boundary
conditions
(1.2) that
iIQ!at=$-q, This equation
(1.5)
expresses the fact that the
because of the heat inflow through the boundary The conservative
difference
of the region.
system for which the difference analogue of Eq. (1 S) is satisfied
has the form
U,-Tii --
Ui+,--2U,iUi-* ‘h’
St
uz-uj
---p= h
h
=
u,--vu*
2At
UN-UN--I +_ h
%
h
In the above, Ui are the values of the mesh function misunderstanding,
the mesh and continuous
i=2,3, . . . , K- 1,
0,
2
Cl .6)
u&i.y At
=
(1.7)
$.
u (below: where this does not lead to a
functions
are designated identically)
in the nodes of
a uniform mesh on a unit section with spacing lz = I I(.% - I), At is the interval with respect to time. and Gi is the temperature
on the previous la),er with respect to time. The difference
(1.6). (I .7) is of the second order of approximation
with respect to h and the first with respect to
At. Because system (1.6), (I .7) is fully implicit, it is unconditionally second order of accuracy along the spacing of a three-dimensional The difference
scheme
stable and therefore
is of
mesh [Z]
analogue of the amount of heat Q has of course to be determined
as follows
(1.8)
Using Eqs. (1.6) and (1.7). we obtain
Q-d,$-q At
(I.91
.
Relationship (1.9) is a difference analogue of Eq. (1.5). Thus, the difference and during its use no additional inputs or outflows of heat are present,
scheme is conservative.
To solve the system of equations (1.6) and (1.7) we use Seidel’s method (see, for example. [l] ). In this case, the calculation formulae of the method are:
134
M. Yu. Shashkov
(‘J
u2
(rtl)
-u,
h
--
(‘+ll
(s+l)
,S+l,
u .x
h u,,’
--U.S.-~ T--
h
Ui
=cF7
*t
2
h
(I+11
ul(s+l’-&
Ui+l
At
_2Ui(s+J)+Uy =
h2
-
(1.10)
zc‘.
=
The imtial approximation iteration
0:
-ii,
At
2
( .),
-iii
Ui co) for the values Ui can be specified arbitrarily; process (I. IO) converges [ 11.
in any case the
Let us write Eqs. (I. 10) in a form analogous to Eqs. (1.6) and (1.7): (St!’
(*+I)
U?
h
-U1
_-
h Ui
(StJJ -la!,
&‘+‘)
1
-ii,
At
2
=
(‘+L.)._2Ui(~+~liU~~f:l)
”
=
At
After carrying out a certain number of iterations convergence,
the process is terminated.
and
ui‘+‘)
(1.1 I)
h2
h2
-
U-kl)
(‘I
u,+t -u,+1
Ui+l
necessary to satisfy the chosen criterion of is taken as the value of the temperature
the next layer. with respect to time. Thus. the mesh function
u,=u,(‘+‘),
obtained
on
by the iteration
process satisfies Eqs. (1.1 1) where the right-hand sides contain terms that are additional
with
respect to the initial system of Eqs. (I .6). (1.7) and which can be regarded as inputs and outflows of heat. Let the condition
for the iteration
process to converge have the form
where E is the prescribed accuracy’. Let us estimate the imbalance Using Eqs. (I .I I) and the definition
It follows that the imbalance
caused by using the Seidel method.
of Q, in the difference case we have
of the amount of heat is determined
by the relationship (1.13)
Solvingdifference equations b_xs iteration methods
Using the convergence
condition
135
we arrive at the estimate
lAQl$ The expression on the right-hand
(N-1)hAt
(1.14)
= $,
side gives the upper limit of the imbalance occurring at one
step with respect to time, The results of calculations the end of the paper the results obtained
by Seidel’s method are given in Section 4. At
are discussed and compared with the estimate (I .14).
2. Condition of conservatism of the two-layer iteration process We write Eqs. (1.6) and (1.7) in the operator form:
5
(2.1)
+Au=/,
where
(Au) I=
i=l, i=2,3,. . , , N-i,
-2 (u,-u,) /h’, -(~d+i-Zui+~i-g)/hZ,
I
i=N,
2(UN--U.M)lh2,
i-2,3,.
If we introduce
,. ,N-1,
a scalar product in the space of mesh functions
(2.2,
by the equation (2.3)
it becomes obvious that the amount of heat Q is given by the equation
where / is the mesh function, The two-layer iteration
which equals unity at all nodes, Ii = I. process for Eq. (2.1) has the form _‘iw
u(a+l)
B T
where B is the given operator,
u”)__ii + Au(')=f,
+ -
(2.4)
At
and r is the iteration
parameter.
We reduce Eq. (2.4) to a form analogous to (1 .ll j, i.e. we shall clarify what kind of inputs and outflows occur as a result of using the iteration as
process. It is clear that Eq. (2.4) can be written
136
M. Yu. Shashkov
(2.5)
where E is the unit matrix. The requirement
that an iteration process should be conservative means that the inputs and
outflows appearing on the right-hand side of (2.5) should compensate
each other and not contribute
to the chan_pe of the amount of heat in the whole system. By the definition
of Q and by Eq. (2.5), the requirement
of conservatism
is equivalent
to the
condition
(u(~+‘)-u(~)), 1) =o,
( [EjAt+A-B/z] which. because of the arbitrariness
of
u(‘), u(‘+‘),
is in turn equivalent
(2.6) to the condition
[ElAt+A-Blz]‘Z=O, or
[ElAt+A*-B’/T]l=O.
(2.7)
Since in our case A =A * and in addition AI = 0, Eq. (2.7) becomes,
[E/At-B./T For example. the method of simple iteration.
(2.8)
]I=O.
for which B = E and 7 = Af. satisfies the above
equation.
3. Flow scheme for the heat conduction The equation
of heat conduction
equation
(I .I) is written in the flow form:
adat-
dz a (WI,
(3.1)
w=adax, where
w
is the heat flow. By Eq. (3.2), the boundary
w 1.=o=rp, Equation (3.1) expresses the law of conservatjon
w+-i=
(3.2) conditions
11‘.
(I .2) take the form (3.3)
of heat. It is desirable to organize the procedure
for the numerical solution of system (3.1), (3.2) so that the difference analogue (3.1) occurs. The flow system proposed in [3-51 satisfies this requirement.
137
Solvingdifference equations by iteration methods
A difference scheme is constructed
Ui-Xj
for the system of equations
Wi+t-Wi
i-$2,.
-=
h
At
Here the temperature the surrounding
(3.4)
. . , N-i,
’
LZf-Ui-1
wi =
(3.1), (3.2). It has the form
h
(3.5)
. . , N-l.
i-2,3,.
’
Ui refers to the centres of the cells, and the flows
wi
and
~i+~
refer to
nodes.
Let us describe the procedure
by which a transition
is made to the next step with respect to
ui using Eq. (3.4) and substitute
time. We express the temperature
Wi=At
wl+t-2W,+Wi-i
+
h2
it into Eq. (3.5). Then we have
ii{-iii_i
(3.6)
h
’
or
!!Lwi+,_ ( l+zg)w,+;w,_,=_ -y-’ Thus, an equation with the boundary
conditions
(3.6) for determining
.
(3.7)
the flow on the next
layer with respect to time is obtained. After an approximate value of the flow is worked out using any iteration process, the approximate temperature on the next layer with respect to time is easily found from (3.4). We shall show that in such an organization the method of solving Eq. (3.7). an imbalance
of the computing
process, which is independent
of
in the amount of heat does not occur. In fact, from
Eq. (3.4) it follows that for’the amount of heat Q which in this case can of course be determined
b>
the equality, N-i
Q = x uil~, z-1 the relationship accurately
(Q-Q)/A~=w~-uJ~
satisfied and redefinition
(Q-@/At-q-q,
holds. Because the boundary
conditions
for the flow are
of the values w,, w.% does not take place, we have
which it was required to prove.
4. Results of calculations We give below the results of calculations boundary
and initial conditions
by the Seidel method for a model problem. The
for Eq. (1.1) are given as follows:
cp=o,q=o,
uo(5)
=i+coszz.
The accurate solution of problem (1.1) (4.1) is
(4.1)
138
M. Yu. Shashkov
Obviously, the amount of heat in this case for any moment of time is
0
The first two calculations
were performed
to clarify the dependence
of the imbalance AQ on
time. Gzlculatior~ 1. The initial approximation graph of Ag as a function
is
u:‘)=2,
At=0.04, h='l,:, ~=0.01.The
oft is the upper line in Fig. 1,
4Q
t
The choice of the initial approximation
is explained by the desire to simulate a case where the
solution varies strongly with time, and therefore there is no good initial approximation Gdcularion 2.
calculation.
u,“‘=ii,;
The result is illustrated
It can be seen from Calculation
the quantities
for ui.
Ar, h and E are the same as in the previous
by the lower curve in Fig. I. 1 that the imbalance builds up in time and reaches 40% of the
true value of the amount of heat. The imbalance occurring during the transition to the next step with respect to time in this calculation is around 7 X 10-2 at each step. The estimate of the upper limit of imbalance,
in accordance with formula (1.14) is 7.84 X 10-2.
Thus, we may say that situations exist where the estimate (1.14) is fairly close, and this estimate should be considered in deciding the standard of accuracy of the iteration process. In Calculation 2 the imbalance is considerably in time it can reach a considerable value.
less, but it also has a tendency
to build up and
Three more calculations were made. Their purpose was to establish the real nature of the dependence of the imbalance on Ar, E and A’. Here the imbalance was calculated for one step with respect to time. In all the calculations which follow ui(O) E 2.
Soh’ing difference equations by iteration methods
139
Calculation 3. The dependence
of AQ on Ar for e = 0.01 and h = l/14 was studied.
Calcularion 4. The dependence
of AQ on E for At = 0.04 and h = l/14 was studied
Calculations
3 and 4 showed that the dependence
of AQ on both At and E has a linear form,
which agrees with the estimate (1.14).
Calculation 5. Figure 2 shows AQ as a function of the number of points N = 1 t l/h for E = 0.01 and At = 0.04. It can be seen that the imbalance increases as the number of nodes in the mesh increases? and this relation is quadratic,
which agrees with the estimate (1 .I 4).
Conclusions When solving difference equations considerable
heat imbalance,
bp non-conservative
iteration
methods the occurrence
of a
which may increase with time, is possible. This should be taken into
account in selecting the quantity
E which governs the accuracy of an interation
process. Calculations
showed that the estimate (1.14) was indeed reached in a number of cases, For example, in order that the law of change in the amount of heat for the difference
scheme
(1.6). (I .7) for the instant of time rk be satisfied with accuracy h2, i.e. with the accuracy of the difference
scheme. the quantity
E should be chosen from the condition
or
h’
&-
%(/?--I)
*
The author thanks A. A. Samarskii for his interest. The author is also grateful to A. P. Favorskii, V. F. Tishkin and P. nl. Babishchev who took part in a discussion of the results obtained
and made a number of valuable suggestions. Translated by W. Chrzczonowicz
REFERENCES 1. SAMARSKII, A. A. and NIKOLAEV, E. S., Methods ofsolving mesh equations (Meted), resheniya setochnykh uravnenii), Nauka, Moscow. 1978. 2. SAMARSKII, A. A., Theory of difference schemes (Teoriya raznostnykh
skhem), Nauka,
Moscow,
1977.
3. TISHKIN, V. F., FAVORSKII, A. P. and SHASHKOV, M. Yu., An algorithm for the numerical solution of the second boundary value problem for the equation of heat conduction on a non-rectangular mesh. Preprint No. 113, IPM AN SSSR. 4. TISHKIN. V. F., FAVORSKII, A. P. and SHASHKOV, M. Yu., Variational-difference systems for the heatconduction equation on irregular meshes, DokL Akad. Nauk SSSR 246, No. 6, 1342-1345, 1978. 5. KORSHIYA, T. K. et al, A variational approach to the construction of difference schemes on curvilinear meshes for the heat-conduction equation. Zh vjbhisl Mat. mat. Fiz., 20, 2.401-421, 1980.
Solving difference eqlcarions b_~,iteration methods 8. DAUTOV. R. Z.. and LAPIN. A. V.. Mesh schemes of arbitrary order of accuracy for quasi-linear elliptic equations. Ix I‘XOI: Ser. marem, No. 10 (209), 24-37, 1980. 9. DAUTOV, R. Z., On superconvergence of mesh schemes for solving the third boundary value problem fol quasi-linear elliptic equations, in: ,Vumerical merhods of the mechanics of continuous media (Chisl. meted) mekhan. sploshnoi sredyj. Vol. 2. No. 1.62-81. Nauka. Novosibirsk, 1980.
10. MAKAROV, V. L., and SAMARSKII. A. A.. Application of exact difference schemes to the estimation of the convergence of the method of straight lines. Zh. v,$hisl Mof. mu?. Filir, 20, No. 2. 371-387. 1980.
11. ODEN. J. T., and REDDY. J. H.. An Introduction to the Mathematical Theor), of Finite Elements, H’ileyInterscience. New York. 1976. 12. TIKHONOV. A. N.. and SAMARSKII, A. A.. On homogeneous difference schemes of high order of accurac!. DokL Akad. Nouk SSSR, 131, No. 3,514-517. 1960. 13. TIKHONOV, A. N.. and SAMARSKII. A. A., Homogeneous difference schemes on non-uniform meshes, Zk v_i&isL Mot. mat. Fii., 2, No. 5, 812-832, 1962. 14. TIKHONOV. A. N.. and SAMARSKII, A. A., On homogeneous difference schemes of high order of accurac! on non-uniform meshes, Zh v~chjsl. Mar. mat. Fiz., 3, No. 1, 99-108, 1963. 15. KORNEEV. V. G.. Mesh operators, energ)-wise equivalent to the operators generated by piecewise Hermitian extensions, Zh y:chisL Mot. mat. Fii., 19, NO. 2,402-416, 1979. 16. KORNEEV, V. G. Finite element method: convergence to piecewise polynomial interpolation of the e\act solurion and convergence on certain sets of points, Part I. Dep. at VINITI, No. 4754-81, DEP? 1981; Part II, No. 4755-81 DEP; Part III, No. 4756-81 DEP. 17. KORNEEV. V. G., On exact mesh schemes. Zh v.TchisLMa?. mot. Fiz., 22. No. 3. 646-654,
1982.
18. BRA\lBLE. J. H., and HILBERT. S. R., Estimation of linear functionals on Sobole\ spaces with apptication to Fourier transforms and \pline interpolation, SIAM, J. IVumer. Anoljxsis. 7. 113- 124. 1970. 19. KORNEEV. V. G.. Schemes oj the finite-element method of high orders of occurac?, (Skhem!, metoda konechn! kh elementov v!,sokikh por!,adkov tochnosti). Izd-vo Leningr. un-ta. Leningrad. 1977. 20. SEARLE. F.. Method offinite
elemenrs for ellipric problems (Russian translation. Mir. \los
1980).
21. KORNEEV. V. G.. On the construction of variational difference schemes of high orders 01‘ accurac!, I’.csm Leningr. un-ta, 4, No. 19. 28-40. 1970.
U.S.S.R. Comput. Moths. Moth Pl7~s. Vol. 22. No. 5. pp. 131-139. Printed in Great Britain
0041-5553/82507.50+.00 Gl964. Perfamon Press Ltd.
1982.
VIOLATION OF CONSERVATION LAWS WHEN SOLVING DIFFERENCE EQUATIONS BY ITERATION METHODS* M. Yu. SHASHKOV Moscow (Received
75 Sepfember 1980; revised 16 Februac,
IT IS SHOWN that the use of iteration to conservative
difference
free of this disadvantage
methods for solving difference
systems may result in a considerable is introduced.
Results of calculations
lZil @chisL Mot. mat. Fi:., 22, 5, 1149-1156.
1982.
198 I )
equations
energy imbalance. are given.
which correspond A class of methods
M Yu Shashkov
132
Introduction When solving difference equations
it is advisable for several reasons to use iteration methods
111. When constructing
a difference scheme one endeavours
to represent correctly, the basic
properties of the physical process described by the differential basic conservation
laws in difference form is an important
Suppose a conservative scheme is constructed difference equations obtained
are solved by an iteration
equations.
requirement
The satisfaction
imposed on the scheme [2].
by some method, and the corresponding
method. The mesh function
(or several functions)
satisfies the difference scheme with a certain accuracy. Consequently,
is conservative,
the solution obtained may not satisfy the conservation
much the magnitude
of the
of the resulting energy. imbalances
even if the scheme
laws. The question of how
depends on the accuracy of the iteration
process. and how much on the parameters of the difference scheme itself (e.g. on the time and space steps). therefore becomes important. imbalance
in multidimensional
It is particularly
important
to know how to assess the
problems where highly accurate calculations
are impossible.
A way out may be to set up iteration processes for which an imbalance does not occur for any itetation
accuracy. Another possibility
enable any iteration
is to form difference schemes of a special kind, which
process in which the imbalance
These possibilities
does not occur to be employed.
are discussed in the present paper using an example of the second boundary
value problem for the one-dimensional
equation
of heat conduction.
Below it is shown that in using the non-conservative
Seidel method there is a considerable
imbalance in the amount of heat. and this imbalance may increase with time. A theoretical assessment of the magnitude of the imbalance is given, and it is compared with the actual)!,occurring value obtained when solving a model problem. The conditions for the conservation of a two-layer iteration process are obtained, and an example of a conservative process is given. A flow scheme for the equation of heat conduction which eliminates conservative properties of iteration processes is described.
1. Example of a non-conservative Consider the second boundary
the problems connected
iteration
with the non-
process
value problem for the one-dimensional
equation
of heat
conduction:
au
PU
at
dU
dl
I
-q(t). IPO
II
are given functions. where cc, 3, u. instant of time t is given by
(1.1)
OCSCI,
ax"
-=-(
au
asI*-i=zC:(t),
(1.2)
1,,,=u,(s). The amount of heat Q contained
(1.3) in the system at the
Solvingdifferewe equations by iteration methods
Q(t)=
133
ju(x,t)dz.
(1.1)
0
It follows from Eq. (I .I) and the boundary
conditions
(1.2) that
iIQ!at=$-q, This equation
(1.5)
expresses the fact that the
because of the heat inflow through the boundary The conservative
difference
of the region.
system for which the difference analogue of Eq. (1 S) is satisfied
has the form
U,-Tii --
Ui+,--2U,iUi-* ‘h’
St
uz-uj
---p= h
h
=
u,--vu*
2At
UN-UN--I +_ h
%
h
In the above, Ui are the values of the mesh function misunderstanding,
the mesh and continuous
i=2,3, . . . , K- 1,
0,
2
Cl .6)
u&i.y At
=
(1.7)
$.
u (below: where this does not lead to a
functions
are designated identically)
in the nodes of
a uniform mesh on a unit section with spacing lz = I I(.% - I), At is the interval with respect to time. and Gi is the temperature
on the previous la),er with respect to time. The difference
(1.6). (I .7) is of the second order of approximation
with respect to h and the first with respect to
At. Because system (1.6), (I .7) is fully implicit, it is unconditionally second order of accuracy along the spacing of a three-dimensional The difference
scheme
stable and therefore
is of
mesh [Z]
analogue of the amount of heat Q has of course to be determined
as follows
(1.8)
Using Eqs. (1.6) and (1.7). we obtain
Q-d,$-q At
(I.91
.
Relationship (1.9) is a difference analogue of Eq. (1.5). Thus, the difference and during its use no additional inputs or outflows of heat are present,
scheme is conservative.
To solve the system of equations (1.6) and (1.7) we use Seidel’s method (see, for example. [l] ). In this case, the calculation formulae of the method are:
134
M. Yu. Shashkov
(‘J
u2
(rtl)
-u,
h
--
(‘+ll
(s+l)
,S+l,
u .x
h u,,’
--U.S.-~ T--
h
Ui
=cF7
*t
2
h
(I+11
ul(s+l’-&
Ui+l
At
_2Ui(s+J)+Uy =
h2
-
(1.10)
zc‘.
=
The imtial approximation iteration
0:
-ii,
At
2
( .),
-iii
Ui co) for the values Ui can be specified arbitrarily; process (I. IO) converges [ 11.
in any case the
Let us write Eqs. (I. 10) in a form analogous to Eqs. (1.6) and (1.7): (St!’
(*+I)
U?
h
-U1
_-
h Ui
(StJJ -la!,
&‘+‘)
1
-ii,
At
2
=
(‘+L.)._2Ui(~+~liU~~f:l)
”
=
At
After carrying out a certain number of iterations convergence,
the process is terminated.
and
ui‘+‘)
(1.1 I)
h2
h2
-
U-kl)
(‘I
u,+t -u,+1
Ui+l
necessary to satisfy the chosen criterion of is taken as the value of the temperature
the next layer. with respect to time. Thus. the mesh function
u,=u,(‘+‘),
obtained
on
by the iteration
process satisfies Eqs. (1.1 1) where the right-hand sides contain terms that are additional
with
respect to the initial system of Eqs. (I .6). (1.7) and which can be regarded as inputs and outflows of heat. Let the condition
for the iteration
process to converge have the form
where E is the prescribed accuracy’. Let us estimate the imbalance Using Eqs. (I .I I) and the definition
It follows that the imbalance
caused by using the Seidel method.
of Q, in the difference case we have
of the amount of heat is determined
by the relationship (1.13)
Solvingdifference equations b_xs iteration methods
Using the convergence
condition
135
we arrive at the estimate
lAQl$ The expression on the right-hand
(N-1)hAt
(1.14)
= $,
side gives the upper limit of the imbalance occurring at one
step with respect to time, The results of calculations the end of the paper the results obtained
by Seidel’s method are given in Section 4. At
are discussed and compared with the estimate (I .14).
2. Condition of conservatism of the two-layer iteration process We write Eqs. (1.6) and (1.7) in the operator form:
5
(2.1)
+Au=/,
where
(Au) I=
i=l, i=2,3,. . , , N-i,
-2 (u,-u,) /h’, -(~d+i-Zui+~i-g)/hZ,
I
i=N,
2(UN--U.M)lh2,
i-2,3,.
If we introduce
,. ,N-1,
a scalar product in the space of mesh functions
(2.2,
by the equation (2.3)
it becomes obvious that the amount of heat Q is given by the equation
where / is the mesh function, The two-layer iteration
which equals unity at all nodes, Ii = I. process for Eq. (2.1) has the form _‘iw
u(a+l)
B T
where B is the given operator,
u”)__ii + Au(')=f,
+ -
(2.4)
At
and r is the iteration
parameter.
We reduce Eq. (2.4) to a form analogous to (1 .ll j, i.e. we shall clarify what kind of inputs and outflows occur as a result of using the iteration as
process. It is clear that Eq. (2.4) can be written
136
M. Yu. Shashkov
(2.5)
where E is the unit matrix. The requirement
that an iteration process should be conservative means that the inputs and
outflows appearing on the right-hand side of (2.5) should compensate
each other and not contribute
to the chan_pe of the amount of heat in the whole system. By the definition
of Q and by Eq. (2.5), the requirement
of conservatism
is equivalent
to the
condition
(u(~+‘)-u(~)), 1) =o,
( [EjAt+A-B/z] which. because of the arbitrariness
of
u(‘), u(‘+‘),
is in turn equivalent
(2.6) to the condition
[ElAt+A-Blz]‘Z=O, or
[ElAt+A*-B’/T]l=O.
(2.7)
Since in our case A =A * and in addition AI = 0, Eq. (2.7) becomes,
[E/At-B./T For example. the method of simple iteration.
(2.8)
]I=O.
for which B = E and 7 = Af. satisfies the above
equation.
3. Flow scheme for the heat conduction The equation
of heat conduction
equation
(I .I) is written in the flow form:
adat-
dz a (WI,
(3.1)
w=adax, where
w
is the heat flow. By Eq. (3.2), the boundary
w 1.=o=rp, Equation (3.1) expresses the law of conservatjon
w+-i=
(3.2) conditions
11‘.
(I .2) take the form (3.3)
of heat. It is desirable to organize the procedure
for the numerical solution of system (3.1), (3.2) so that the difference analogue (3.1) occurs. The flow system proposed in [3-51 satisfies this requirement.
137
Solvingdifference equations by iteration methods
A difference scheme is constructed
Ui-Xj
for the system of equations
Wi+t-Wi
i-$2,.
-=
h
At
Here the temperature the surrounding
(3.4)
. . , N-i,
’
LZf-Ui-1
wi =
(3.1), (3.2). It has the form
h
(3.5)
. . , N-l.
i-2,3,.
’
Ui refers to the centres of the cells, and the flows
wi
and
~i+~
refer to
nodes.
Let us describe the procedure
by which a transition
is made to the next step with respect to
ui using Eq. (3.4) and substitute
time. We express the temperature
Wi=At
wl+t-2W,+Wi-i
+
h2
it into Eq. (3.5). Then we have
ii{-iii_i
(3.6)
h
’
or
!!Lwi+,_ ( l+zg)w,+;w,_,=_ -y-’ Thus, an equation with the boundary
conditions
(3.6) for determining
.
(3.7)
the flow on the next
layer with respect to time is obtained. After an approximate value of the flow is worked out using any iteration process, the approximate temperature on the next layer with respect to time is easily found from (3.4). We shall show that in such an organization the method of solving Eq. (3.7). an imbalance
of the computing
process, which is independent
of
in the amount of heat does not occur. In fact, from
Eq. (3.4) it follows that for’the amount of heat Q which in this case can of course be determined
b>
the equality, N-i
Q = x uil~, z-1 the relationship accurately
(Q-Q)/A~=w~-uJ~
satisfied and redefinition
(Q-@/At-q-q,
holds. Because the boundary
conditions
for the flow are
of the values w,, w.% does not take place, we have
which it was required to prove.
4. Results of calculations We give below the results of calculations boundary
and initial conditions
by the Seidel method for a model problem. The
for Eq. (1.1) are given as follows:
cp=o,q=o,
uo(5)
=i+coszz.
The accurate solution of problem (1.1) (4.1) is
(4.1)
138
M. Yu. Shashkov
Obviously, the amount of heat in this case for any moment of time is
0
The first two calculations
were performed
to clarify the dependence
of the imbalance AQ on
time. Gzlculatior~ 1. The initial approximation graph of Ag as a function
is
u:‘)=2,
At=0.04, h='l,:, ~=0.01.The
oft is the upper line in Fig. 1,
4Q
t
The choice of the initial approximation
is explained by the desire to simulate a case where the
solution varies strongly with time, and therefore there is no good initial approximation Gdcularion 2.
calculation.
u,“‘=ii,;
The result is illustrated
It can be seen from Calculation
the quantities
for ui.
Ar, h and E are the same as in the previous
by the lower curve in Fig. I. 1 that the imbalance builds up in time and reaches 40% of the
true value of the amount of heat. The imbalance occurring during the transition to the next step with respect to time in this calculation is around 7 X 10-2 at each step. The estimate of the upper limit of imbalance,
in accordance with formula (1.14) is 7.84 X 10-2.
Thus, we may say that situations exist where the estimate (1.14) is fairly close, and this estimate should be considered in deciding the standard of accuracy of the iteration process. In Calculation 2 the imbalance is considerably in time it can reach a considerable value.
less, but it also has a tendency
to build up and
Three more calculations were made. Their purpose was to establish the real nature of the dependence of the imbalance on Ar, E and A’. Here the imbalance was calculated for one step with respect to time. In all the calculations which follow ui(O) E 2.
Soh’ing difference equations by iteration methods
139
Calculation 3. The dependence
of AQ on Ar for e = 0.01 and h = l/14 was studied.
Calcularion 4. The dependence
of AQ on E for At = 0.04 and h = l/14 was studied
Calculations
3 and 4 showed that the dependence
of AQ on both At and E has a linear form,
which agrees with the estimate (1.14).
Calculation 5. Figure 2 shows AQ as a function of the number of points N = 1 t l/h for E = 0.01 and At = 0.04. It can be seen that the imbalance increases as the number of nodes in the mesh increases? and this relation is quadratic,
which agrees with the estimate (1 .I 4).
Conclusions When solving difference equations considerable
heat imbalance,
bp non-conservative
iteration
methods the occurrence
of a
which may increase with time, is possible. This should be taken into
account in selecting the quantity
E which governs the accuracy of an interation
process. Calculations
showed that the estimate (1.14) was indeed reached in a number of cases, For example, in order that the law of change in the amount of heat for the difference
scheme
(1.6). (I .7) for the instant of time rk be satisfied with accuracy h2, i.e. with the accuracy of the difference
scheme. the quantity
E should be chosen from the condition
or
h’
&-
%(/?--I)
*
The author thanks A. A. Samarskii for his interest. The author is also grateful to A. P. Favorskii, V. F. Tishkin and P. nl. Babishchev who took part in a discussion of the results obtained
and made a number of valuable suggestions. Translated by W. Chrzczonowicz
REFERENCES 1. SAMARSKII, A. A. and NIKOLAEV, E. S., Methods ofsolving mesh equations (Meted), resheniya setochnykh uravnenii), Nauka, Moscow. 1978. 2. SAMARSKII, A. A., Theory of difference schemes (Teoriya raznostnykh
skhem), Nauka,
Moscow,
1977.
3. TISHKIN, V. F., FAVORSKII, A. P. and SHASHKOV, M. Yu., An algorithm for the numerical solution of the second boundary value problem for the equation of heat conduction on a non-rectangular mesh. Preprint No. 113, IPM AN SSSR. 4. TISHKIN. V. F., FAVORSKII, A. P. and SHASHKOV, M. Yu., Variational-difference systems for the heatconduction equation on irregular meshes, DokL Akad. Nauk SSSR 246, No. 6, 1342-1345, 1978. 5. KORSHIYA, T. K. et al, A variational approach to the construction of difference schemes on curvilinear meshes for the heat-conduction equation. Zh vjbhisl Mat. mat. Fiz., 20, 2.401-421, 1980.