This third Note is dedicated to the discussion of application of the 3D Douglas – Rachford ADI scheme to the solution of a non-linear heat equation. We will discuss the Newton – Raphson method and the “method of frozen coefficients”.
Before starting discussion of how to approach solution of a non-linear equation, let us define what will be called a non-linear heat equation.
Below, we will consider a heat equation with temperature dependent thermal conductivity and heat capacity. Thus, the coefficients of the heat equation appear to be temperature dependent:
, (1)
where 
and 
are heat capacity and thermal conductivity respectively. When we talk about non-linearity, we mean it in the above sense – i.e., the coefficients of heat equation depend on temperature (according to [(Lees, 1966)], such equation should be classified as quasi-linear).
Now, let us consider the 3D Douglas – Rachford scheme for the Eq. (1) [(Douglas, 1962)], [(Lees, 1961)]:
, (2)
, (3)
, (4)
where 
- is the unknown function, defined on the following mesh:
(5)
where
, (6)
,(7)
- are the indexes of nodes of the spatial mesh along the 
and 
directions respectively. 
- are the total numbers of nodes of the spatial mesh along the 
and 
directions respectively. 
- are the steps values along the directions 
and 
respectively. 
- is the value of time step. 
and 
- are the intermediate quantities, that are necessary for the calculation of unknown function in the scheme (2) – (4).
First of all, let us discuss the time moments, at which the coefficients 
and 
of the scheme (2) – (4) are taken. These coefficients depend on the unknown function 
and thus, depend implicitly on the time moment at which its values are calculated.
The following three cases can be imagined:
1. If one takes the values 
and 
(i.e., the values at time moment 
) for each entrance of these coefficients in scheme (2) – (4), one obtains the so-called “scheme with frozen coefficients”. Although at each time step the coefficients 
and 
must be re-calculated, both 
and 
are known and the scheme becomes linearized.
2. One could erroneously suggest to consider the following scheme, where the coefficients 
and 
are calculated at “intermediate” time moments – i.e., the intermediate quantities 
and 
are substituted as arguments in and in order to obtain 
, 
, 
,
:
, (8)
, (9)
, (10)
This would be an error because the quantities 
and 
are not “true” values of the unknown function at intermediate time moments. This can be seen from the fact that 
is the solution for Eq. (8) (which is an implicit scheme with respect to Y-direction) and 
is the solution for Eq. (9) (which is an implicit scheme with respect to X-direction), but none of 
and 
is the solution for the full heat equation (1). Contrary, 
and 
are the solutions to Eq. (1), defined on the mesh (5) – (7). Thus, this second idea leads to an erroneous result and must be avoided.
3. Finally, one can take the values 
and 
(i.e., the values at time moment 
)for each entrance of these coefficients in scheme (2) – (4). This will correspond to a “truly” non-linear equation (1) and one has to treat somehow this nonlinearity.
In what follows, we are going to consider the application of the Newton – Raphson method to the system (2) – (4) keeping in mind the 3rd idea discussed above. Let us suppose that the system (2) – (4) can be regarded as an ordinary system of non-linear equations, to which the Newton – Raphson method is applicable. Under such assumption, one can write for the system (2) – (4) the following stiffness equation:
, (11)
where 
- is the jacobian matrix
, (12)
, (13)
, (14)
, (15)
, (16)
, (17)
, (18)
, (19)
, (20)
, (21)
, (22)
, (23)
, (24)
, (25)
, (26)
, (27)
, (28)
, (29)
where 
, 
and 
are the initial guesses to the values of 
, 
and 
; 
- is the number of iteration.
Then,
. (30)
Each element in the Eq. (30) is a 
matrix. Thus 
- is the 9-diagonal 
matrix.
Now, when the (12) and (14) are found, one can solve (11) with respect to (13). Note that in such formulation this is no more an ADI method, because one has to treat the matrices that are not 3-diagonal. Moreover, one can not reduce this system to an ADI Douglas – Rachford scheme. Indeed, if one sums up the equations of the system (11) – (14), one would not obtain a Newton – Raphson equation for an implicit finite difference scheme. And an attempt to eliminate the intermediate quantities 
and 
from the system (11) – (14) would lead to an equation that has nothing in common with the ADI method.
Thus, we conclude that the above mentioned ideas Num. 2 and Num. 3 are erroneous for the solution of a heat equation with the ADI method, and the only justified method for the account of variable coefficients is the “method of frozen coefficients”. In practical applications either the method of “frozen coefficients” [(J.E. Dendy, 1977)] or some approximation of nonlinear coefficients (e.g. – Douglas and Dupont approximation [(Vijinder K. Bangia, Curtis Pennett, Albert Reynolds, Rajagopal Raghavan, 1978)], [(J. Douglas, T. Dupont, 1970, v. 7)]) is used.
The presented analysis would be incomplete without the discussion of numerical stability of the used finite difference schemas.
It is well known that the scheme (2) – (4) with the coefficients 
and 
being taken at the time moment k can be rewritten in the following form [(Peaceman, 1977)]::
, (31)
where 
-is the unknown function (temperature field), which is defined in the following mesh:
(32)
where
, (33)
, (34)
are the indexes of the nodes of spatial mesh along the directions 
and 
respectively. 
-are the numbers of nodes in the spatial mesh along the directions 
and 
respectively. 
- are the values of spatial steps along the directions 
and 
respectively. 
- is the value of time step.
-is the finite difference time derivative of the unknown function,
(35)
- is the value of heat capacity in time moment 
. That is, in accordance with the principle of “frozen coefficients”, is a function of spatial coordinates and does not depend on temperature now.
-is the finite difference operator that corresponds to a spatial gradient.
- is the finite difference operator, (36)
- is the value of thermal conductivity in time moment 
. That is, in accordance with the principle of “frozen coefficients”, is a function of space variables and does not depend on temperature now.
Eq. (31) can factorized in the following way:
(37)
Or, the same,
(38)
The scheme (38) is numerically stable if the absolute values of the eigenvalues of operator
(39)
are less than unity at each time level. This generally imposes some constraints on the values and functional forms of heat capacity 
and thermal conductivity 
. To the best of our knowledge, no proof of numerical stability and convergence has been given for the scheme (31) – (39) yet. Moreover, we are not aware of published results on the application of a 3D Douglas – Rachford scheme to the solutions of a non-linear (quasi-linear) heat equation. We attribute this to high difficulty one encounters when analyzing the numerical stability of the algorithm. As it is has been noted by J. Douglas “numerical experience with previous alternating direction methods indicates that these applications (application of the ADI methods to non-linear problems - authors note) should be successful, although no proofs have been constructed” [(Douglas, 1962)].
References:
Douglas, J. (1962). Alternating direction methods for three space variables // Numerische Mathematik. 4 (1, pp. 41-63).
J. Douglas, T. Dupont. (1970, v. 7). Galerkin methods for parabolic problems. SIAM J. Numer. Anal. , 576 - 626.
J.E. Dendy, J. (1977). Alternating direction methods for nonlinear time-dependent problems. SIAM J. Numer. Anal. , 14, 313-326.
Lees, M. (1966). A linear three-level difference scheme for quasilinear parabolic equations. J. Math. Comput. , v. 20, pp. 516 - 522.
Lees, M. (1961). Alternating Direction and Semi-Explicit Difference Methods for Parabolic Partial Differential Equations // Numerische Mathematik. 3, pp. 398 - 412.
Peaceman, D. (1977). Fundamentals of Numerical Reservoir Simulation. Amsterdam: Elsevier.
Vijinder K. Bangia, Curtis Pennett, Albert Reynolds, Rajagopal Raghavan. (1978). Alternating direction collocation methods for simulating reservoir performance. Technical Conference and Exhibition of the Society of Petroleum Engineers of AIME. Houston.