The ADI (alternate directions implicit) method is widely used for the numerical solution of multidimensional parabolic PDE (partial differential equations). . Although the method is known for a long time and is well described in the text-books, its practical realizations sometimes appear to be inaccurate . The inaccuracy arises every time when one neglects the correct account of the so-called intermediate boundary conditions. This neglect can become the cause of instabilities even when the used ADI-scheme is known to be unconditionally stable in the frame of von Neumann spectral analysis technique . The procedures of a correct account of the intermediate boundary conditions (for the Peaceman-Rachford ADI-scheme), are described in [4, 5, 6].
Below, we are going to consider the correct account of the intermediate boundary conditions for the Douglas-Rachford ADI-scheme [3, 7]:
where -is the Laplace operator; - is the unknown function, defined on the following mesh:
i, j, sare the nodes numbers of the spatial mesh along the Y, X and Z directions respectively. - are the total numbers of nodes of the spatial mesh along the X, Y and Z directions respectively. - are the steps values along the directions X, Y and Z respectively. - is the value of time sub-step : the full time step is equal to .
Let the following boundary conditions are set on the edges of a rectangular modeling region (we will not precise the type of the boundary conditions):
Let us now re-write the Eqs. (1), (2) and (3) in such a form, that the unknowns , and , stand in l.h.s.:
and determine which quantities are not known explicitly. For this purpose, let us consider the sets, on which the Eqs (13) – (15) are defined (i.e., the ranges of definition).
As one can see, Eq. (13) contains operators and . For this reason, one cannot calculate the expressions , , and , retaining the approximation order ) that is inherent for the Douglas-Rachford scheme. Thus, the range of definition of Eq. (13) is (see. Figure 1):
Fig. 1 – Schematic illustration of the modeling region. The internal parallelepiped corresponds to the region defined by Eq. (16), in which Eq. (13) is solved. Red lines are the contours of planes, in the points of which the obtained do not possess physical sense (see. discussion below)..
In the r.h.s. of Eq. (14) and (15) we have finite difference operators that refer only to those spatial directions, along which the quasi-one-dimensional problem is solved. Thus, the range of definition of Eqs. (14) and (15) is the whole modeling region:
Solving the Eq. (13), one finds (, , the number 1 over means that the quantity if found from the first equation of the Douglas-Rachford scheme) with the indices i, j, s from Eq. (16).
Here one must ask the following question – what one should consider as a boundary conditions on the edges with and on the intermediate time step ? Namely the quantities with and are necessary for the calculation of using the Eq. (14) – see. Figure 2.
Fig. 2 – Schematic illustration of the modeling region. Red lines are the contours of planes, in the points of which one must know the correct values of boundary conditions to perform calculations using Eq. (14).
At the same time, (please, attention here!) The values at and at , which are calculated using Eq. (13) – are not correct boundary conditions on the time sub-step !
This follows from the fact that Eq. (13) provides ,in which only the time evolution along the Y – axis is taken into account. And this quantity - , - does not contain correct information about the boundary conditions on the time sub-step .
Nevertheless, we have everything we need for the determination of correct boundary conditions at and at (i.e., we implicitly have all the necessary information).
Let us find using Eq. (14):
adn , (18)
аnd using Eq. (15):
Using Eq. (19) one can find the values of at
because, according to Eq. (4), the correct boundary conditions are known on every time steps (not intermediate sub-steps).
Thus, knowing the boundary values of , using Eq. (18) one can calculate the correct boundary values of at and at . And the latter, in their turn, can be used to determine the values of вthe internal points of the modeling region. Then, using Eq. (15), one finds for the all internal points, thus passing on the next time layer (see Eq. (4)). And then everything repeats for every time layer.
The procedure described above is valid for the case when the heat equation is linear and the boundary conditions are also linear. The situation becomes much more complicated when one solves non-linear equations, e.g. – heat equation with temperature dependent thermal conductivity and heat capacity. Solving a non-linear equation is a complicated problem itself, and when one has to consider in addition non-linear equations arising from the correct account of boundary conditions on the intermediate time layers in an ADI-scheme, the problem becomes almost unsolvable.
Let us now consider such a problem a bit closer.
Re-write the Douglas-Rachford scheme (1) – (3) in the following form, the range of definition being defined as before by Eq. (4):
where - is the temperature dependent heat capacity.
One can apply to Eqs. (21) – (23) the described above procedure of determination of intermediate boundary conditions. In this case, a complicated functional dependence of the heat capacity on temperature leads to appearance of the following terms in Eqs. (21) – (23):
1. In Eq. (22):
2. 2. In Eq. (23):
where - are the indices of the nodes on the boundaries of modeling region.
Thus, in order to obtain and , one has to solve equations containing terms of the form (24) – (25). This problem can, in principle, be solved numerically if one possesses a good initial approximation. But the solution can require significant time resources and further analysis of obtained roots (there can be several roots among which one would need to make an appropriate choice).
In conclusion, we would note that the problem of a correct account of intermediate boundary conditions – is the price one has to pay for the time efficiency of an ADI-scheme.
For many cases it has been shown that the appropriate values of the input parameters and consistent set of initial and boundary conditions lead the error (associated with the incorrect account of intermediate boundary conditions) that does not exceed degree C. Nevertheless, the risk of instability – increase of the error’s magnitude when passing from one time layer to another, is still present. Such an unstable behavior can potentially arise at some configurations of the boundary conditions and some values of the input parameters. In this regard, the general analysis of numerical stability of an applied finite difference scheme acquires the crucial importance.
The author acknowledges Prof. T.I. Lakoba for helpful discussions.
1. Ch. Gao , Y. Wang “A general formulation of Peaceman and Rachford ADI method for the N-dimensional heat diffusion equation.” Int. Comm. Heat Mass Transfer, Vol. 23, No. 6, pp. 845 – 854 (1996)
2. W.Y. Yang, W. Cao, T.-S. Chung, J. Morris “Applied numerical methods using MATLAB”. Wiley-Interscience (2005).
3. Peaceman D.W. Fundamentals of Numerical Reservoir Simulation. - Amsterdam : Elsevier, 1977.
4. Strikwerda, J. (2004). Finite difference schemes and partial differential equations. Madison: SIAM.
5. Lakoba, T. (2012). Lecture notes on Numerical Differential Equations. Vermont, USA: University of Vermont. – internet resource: http://www.cems.uvm.edu/~lakobati/.
6. A.R. Gourlay, A.R. Mitchell Intermediate boundary conditions for split operator methods in three dimensions BIT Numerical Mathematics. v. 7 (1). (1967), pp. 31 – 38.
7. Douglas J. Alternating direction methods for three space variables // Numerische Mathematik. – v. 4 (1). (1962), pp. 41-63 : Т. 4.