Analyse numérique des inéquations variationelles (in French), R. Glowinski, J. L. Lions and R. Trémolieres, Dunod, Paris, 1976. No. of pages: 558. price: 390FF

A new algorithm is proposed for the approximate realization of an implicit difference scheme for the one-phase Stefan problem formulated in the form of a variational inequality. This algorithm is theoretically based on the estimates of movement rate for the phase transition mesh boundary and on the estimates of the boundary perturbation impact on the solution of a difference scheme. The approximate realization of a difference scheme consists in its replacement by two problems in overlapping mesh subdomains, one of them being a linear mesh problem and, the other, a nonlinear problem. To solve each of these problems, a method is then proposed to decompose the domain into small subdomains. The Stefan problem and numerical methods of its solution are investigated in a large number of publications (see [5-7, 11-13] and the references therein). Many authors investigated, among other things, the convergence and accuracy of difference schemes. This paper treats a very simple implicit difference scheme on a uniform mesh for the one-phase Stefan problem. To find its solution at each time step, it is necessary to solve a nonlinear problem, namely, a mesh obstacle problem. The paper suggests a new method of the approximate solution of this problem, which is essentially as follows. The original mesh domain is decomposed into a set of overlapping subdomains of small dimensions. In the subdomains, mesh problems are solved independently. The linear difference schemes constitute the larger part of these problems, while mesh obstacle problems, the smaller part. The linear problems can be solved by efficient direct or iterative methods. To solve nonlinear problems, one can use different iterative methods [7], the splitting method [9,10] in particular. In Section 1, we formulate the one-phase Stefan problem in the form of a variational inequality and discuss the known properties of its solution. In Section 2, we construct a difference scheme and investigate the movement rate for the phase transition mesh boundary. We establish here the estimate of the maximal distance at which this mesh boundary can be moved in one time step. In Section 3, we carry out an approximate splitting of the nonlinear mesh problem for a fixed time step into a linear and a nonlinear subproblems in two overlapping mesh subdomains. In this case, the mesh subdomain for the nonlinear subproblem is located in the neighbourhood of the mesh boundary of the phase transition from the preceding time step. In Section 4, we consider the nonlinear mesh subproblem. To solve it, following [8] we use an iteration-free method of decomposing a domain into small subdomains and give its justification. In Section 5, we outline the method and consider certain implementation aspects. 488 Yu. A. Kuznetsov and A. V. Lapin 1. A VARIATIONAL FORMULATION OF THE ONE-PHASE STEFAN PROBLEM Let Ω = (0,1) χ (0,1) be the unit square with the boundary Γ = (Jf=1Fi, where Γ1 = [xi = 0}, F2 = (x2 = 0} u (x2 = 1}, F3 = {xl = 1}. We seek functions 0(x, t) and /(x) satisfying the following relations (see [5,13] and others): -—Δθ = 0, χεΩ, t>l(x) dt 0 = 0, χθΩ, t^l(x) V0.V/=-/c, Χ6Ω, ί = /(χ) 0 = 0!(χ2), xeF1} ί > 0 (1.1) 0 = 0, xeF3, ί > 0 de Λ ^ — = 0, xeF2, ί > 0 0(x,0) = 0, Χ6Ω where η is the vector of the normal to Ω, θt is a prescribed positive function, and k = const > 0. The problem formulated is a mathematical model of the process of melting of ice with zero temperature, which occupies the domain Ω at the time moment t = 0. Here, θ — 0(x, t) is the temperature, t = /(x) is the equation for the boundary of the water-ice phase transition, k is the specific melting heat. At every time moment 0 < ί < T(for T > 0) the domain Ω consists of two subdomains: Ω+(ί) = {χεΩ| θ(χ, ί) > 0} occupied by water Ω0(ί) = {χθΩ| 0(x, t) = 0} occupied by ice. The boundary 3Ω+(ί)ηΩ separating these subdomains is the phase transition boundary. Following [4], define the function Ρ 0(*, Jl(x) "(*' = > ο, ,</(*). τ) dr, t > /(x) (1.2) Assume that 9leW\l(r]}r^C(Yl\ and define V = {veH(Si)\ v(x) = 0 on F3} K(t) = {veV\v(x)^QmQ,, v(x) = 61(x)t on FJ. (1.3) The function Μ is a solution to the following variational inequality: for almost all ie[0, T]: u(t)eK(t) and

Domain decomposition method to realize an implicit difference scheme for the one-phase Stefan problem

New algorithmes are proposed for the approximate realization of implicit mesh schemes (constructed via f.e.rn. or f.d.m.) for a class of singular parabolic equations including, among others, one-phase and two-phase Stefan problems. These algorithmes are theoretically based on the estimates of movement rate for the mesh free boundary and on the estimates of the boundary perturbation impact on the solution of a mesh scheme. The approximate realization of a mesh scheme consists in its replacement by two problems in overlapping mesh subdomains, one of them being linear mesh problem and, the other, a nonlinear problem. To solve each of these problems, a method of decomposition of the domain into small domains can be then used.

Domain decomposition method to realize implicit mesh schemes for some singular parabolic equations

Direct methods are constructed for solving quadratic programming problems with the tridiagonal Jacoby M-matrix and simple constraints. These methods are employed in block relaxation and splitting algorithms for solving two-dimensional mesh obstacle problems. The convergence of the algorithms is investigated and the results of numerical experiments are presented. We consider a difference approximation of the obstacle problem min<|Vw| dx— /udx>, X = {u^W2\\^)\\u(x)^Q9 xeQ}. (0.1) As a result, we arrive at a quadratic programming problem with a sparse symmetric positive definite matrix. To solve this problem in a one-dimensional case where the matrix is tridiagonal, we construct a direct method which requires O(N) arithmetic operations with N unknowns, i.e. an asymptotically (in N) optimal method. In a multidimensional case, we employ block relaxation and splitting methods based on successive solution of one-dimensional mesh problems. The quadratic programming problem minU(Aj^-(/,jok K = {yeR|y,>0, i= 1,...,N} (0.2) y*K t 3 with the symmetric positive matrix A can be solved by well-known methods. Among them we can single out conjugate gradient methods developed for the sparse matrix [11] as well. In fact, these methods consist in determining a set of indices CD' = {i \\ yt > 0 } at which the equalities (Ay —f\\ = 0 hold. The conjugate gradient method is used to solve at each step a subsystem of the simultaneous equations Ay =/. Conjugate gradient methods converge in a finite number of iterations only due to the finite number of subspaces in the space R but efficient estimates of the number of iterations, required for attaining a prescribed accuracy, are unknown. Methods are also available ([2,10] etc.) for solving a complementarity system [equivalent to (0.2)] tt^O, (Ay-f^O, 3^3;-A = 0, i=l,. . . ,N (0.2') which exploit the fact that the matrix A belongs to classes (P) and (Z)[4], i.e. A is an M-matrix [16]. Unlike the conjugate gradient methods, these methods provide for a special construction of approximations to the set ω', with the exact solution y to problem (0.2) found at a step q < N. t Originally published in Russian as Preprint No. 140 of the Department of Numerical Mathematics of the USSR Academy of Sciences, Moscow, 1986. 450 A. V. Lapin and D. O. Solovyev Despite roundoff errors, these methods converge in a finite number of iterations, and their efficiency is determined by the efficiency of solution of systems of linear algebraic equations (Ay —/)f = 0 (ie<w, k = 0,1,.. .)· In the case where A is a banded matrix, it is sufficient at each step to compute the values y only at some points of the set ω*, which are adjacent to the boundary. This fact is taken into account in the algorithm constructed in Section 1 of the present paper for mesh approximation (0.2) of one-dimensional problem (0.1). To find the exact solution to the problem, this algorithm requires O(N) arithmetic operations (N is the number of unknowns). For the sake of comparison, note that for this problem the complete solution of the system of equations by the method [10] requires O(N) arithmetic operations. In Section 1 we also consider a one-dimensional mesh problem with two constraints Φί < J>i ̂ </Ί· (* = 1> · · · >N)· Such problems can be solved by known iterative methods [14] (in addition to the conjugate gradient methods mentioned above) which generate monotonic sequences of approximations to the exact solution (without convergence rate estimates). We suggest a direct method for this problem; each step of the method consists in solving the problem with one of the constraints and requires O(N) operations. The number of such steps can generally be a magnitude of order JV, so that the rough estimate of the total number of arithmetic operations required for finding the exact solution is of the order N. However, this number was small and independent of N in all numerical experiments carried out by the authors. Efficient algorithms constructed in Section 1 for solving one-dimensional problems can be employed for solving multidimensional problems. Thus, in Section 2 we analyze, for mesh approximation of problem (0.1), the splitting and block relaxation methods in which blocks are made up by choosing values of the mesh function on the one-dimensional mesh line. The convergence of these methods is implied by the results of [1,9]. The monotonicity is established of a sequence of iterations of the block relaxation method with the parameter σ = 1 and the convergence rate estimate is obtained for a version of the splitting method. All the algorithms suggested have been numerically tested. 1. ONE-DIMENSIONAL MESH PROBLEMS WITH OBSTACLES INSIDE THE DOMAIN 1.1. Let A be a symmetric M-matrix [16], i.e. A' >0 and α0·<0 for i+j. Let us consider a quadratic programming problem

A. V. Lapin

Papers

Solution of mesh obstacle problems