# Convergence of a lattice numerical method for a boundary-value problem with free boundary and nonlinear Neumann boundary conditions.

1. Introduction. We consider the boundary-value problem,

(1.1) [[partial derivative].sub.t]c(t, x) = [partial derivative].sub.xx]c{t, x), x [member of] ([rho](t),L), t [member of] (0,T),

(1.2) [[partial derivative].sub.x] c (t, L) = G(c(t, L)), t [member of] [0,T],

(1.3) [[partial derivative].sub.x] c(t, [rho](t)) = g(c(t, p(t))), t [member of] [0,T],

(1.4) c(t, [rho](t))[??] = -g(c(t, [rho](t))), [rho](0) = [[rho].sub.0], [rho](T) = 0,

(1.5) c(0, x) = [phi](x), x [member of] [[[rho].sub.0], L], 0 < [[rho].sub.0] < L.

The following assumptions are made: the functions G(c) and g(c) are defined for c > 0; the continuous derivatives G'(c) < 0, g'(c) > 0 for c > 0 exist; the conditions G(1) > 0, G([infinity]) < 0, and g(1) = 0 hold; the function [phi](x) is defined for x G [[[rho].sub.0], L], p(x) > 1, G([phi](x)) > 0 for x [member of] [[[rho].sub.0], L], and [rho]" (x) exists and is continuous in [[[rho].sub.0], L]. Note that the initial-boundary compatibility condition is included in the problem: conditions (1.2), (1.3) hold for t = 0. The following consequence can be derived.

Proposition 1.1. The initial distribution [phi](x) is not constant.

Proof. If [phi](x) = const, then g([phi]) = G([phi]) = 0 due to (1.2) and (1.3) and thus [phi] = 1; but G(1) [not equal to]0.

Let us define the set [Y.sub.[rho]](T) c [R.sup.2]; it consists of all (t, x) such that t [member of] [0,T], x [member of] [[rho](t), L]. Its closure is [[bar.Y].sub.[rho]](T). To solve the problem, we must find the number T, [rho](t) [member of] [C.sup.1] ([0,T]), and c(t, x) [member of] C([Y.sub.[rho]](T)). These functions must satisfy (1.1)-(1.5). If (1.1) is satisfied in a weak sense, then the solution obtained is called the weak solution.

This problem is a generalization of the model of isothermal hydriding of a metal particle under constant pressure [1]. Equation (1.1) is the diffusion equation; without loss of generality we assume that the diffusivity is unity. The nonlinear Neumann boundary condition (1.2) connects the diffusion flux near the surface with sorption and desorption on the surface. The stoichiometric concentration in hydride is unity. The condition G(1) > 0 means that the pressure is above equilibrium with respect to stoichiometric hydride. The condition (1.3) connects the diffusion flux on the phase boundary with the concentration near it and describes the hydride formation process. For stoichiometric concentration hydriding is impossible. The Stefan-type condition (1.4) is the conservation on the free boundary. The initial concentration (1.5) is between the stoichiometric and the equilibrium values.

The main difficulties are the free boundary and rather general nonlinear boundary conditions. Analytical solution methods [2] can hardly be applied to such problems. A lot of research has been devoted to the numerics of free boundary problems (e.g., [3, 4, 5, 6]), but little attention has been paid to diffusion boundary-value problems with free boundary in finite domains and with nonlinear boundary conditions. The aim of this paper is to apply the well-known idea of an implicit difference scheme with time steps dependent on the velocity of the free boundary (see [2] and references therein) to this class of problems and to prove the convergence to a classical solution. Therefore, we prove the existence of the solution in a constructive way. We prove the maximum principle for the discrete lattice problems and use it to construct a convergent sequence of approximations of a solution. The maximum principle for parabolic PDEs and the corresponding discrete systems is well-known [7, 8], but in our case it also holds on the boundary, i.e., boundary conditions also do not allow values of the solution to be too high. The restrictions are rather weak and have a physical meaning. Nevertheless, they are sufficient for the results obtained.

The structure of the paper is the following. First, we construct the difference scheme for this problem and prove a few statements for the lattice solution. Next, we obtain the sequences of continuous approximations to the free boundary and the concentration by linear interpolation. Then, we show that the sequence for the free boundary converges uniformly to some function [rho](t), [rho](T) = 0 for some T; also the sequence for the concentration converges in C([[bar.Y].sub.[rho]](T)). Next, we show that these functions are actually the classical solution of the problem. Therefore, we prove the existence of the classical solution of the nonlinear boundary-value problem with free boundary and additionally justify the difference scheme for the problem. The idea is from [2, 7].

The solution c(t, x) must be positive (because it represents a concentration); moreover, c(t, x) > 1 because the concentration in hydride cannot be below the stoichiometric. We show that the solution indeed has these properties; see Theorem 3.4 below.

2. Difference approximations. Let us divide [[[rho].sub.0], L] in M pieces of equal length h = (L - [[rho].sub.0])/M. In the sequel h [right arrow] 0 means M [right arrow] [infinity]. Let I = [??]/h[??] (integer part), K = I - M, and [[delta].sub.h] = L - Ih. Choose any sequence [K.sub.n] such that [K.sub.0] = K, for n [greater than or equal to] 0 either [K.sub.n+1] = [K.sub.n] or [K.sub.n+1] = [K.sub.n] + 1. Let us denote [k.sub.n] = [K.sub.n] - n.

Now let us consider any spatially uniform lattice with the steps h and [T.sub.n], 0 [less than or equal to] n [less than or equal to] N, 0 [less than or equal to] i [less than or equal to] I. Here N is the minimal N such that [[summation].sup.N.sub.0] [T.sub.n] [greater than or equal to] T, T > 0 is some given time. Let the nodes of this lattice be (n, i); these are the points ([t.sub.n], [x.sub.i]), [x.sub.i] = [[delta].sub.h] + ih, [t.sub.0] = 0, [t.sub.n+1] = [t.sub.n] + [[tau].sub.n], n > 0.

Let [[bar.D].sub.N] be the lattice subset, such that, if (n, i) [member of] [[bar.D].sub.N], then 0 [less than or equal to] n [less than or equal to] N and 1 = [k.sub.n], [k.sub.n] + 1, ..., I. In other words, each succeeding layer contains either the same number of nodes (if [K.sub.n+1] = [K.sub.n] + 1) or one left node more (if [K.sub.n+1] = [K.sub.n]). Let [D.sub.N] be the open subset; if (n, i) [member of] [D.sub.N] then 0 < n < N and i = [k.sub.n] + 1, [k.sub.n] + 2, ..., I - 1. We need one more lattice set [[??].sub.N]; it is 5N without the corner nodes (0, K), (0, I).

We denote the value of the lattice function f at the node (n, i) by [f.sup.i.sub.n]. Let us approximate the derivatives,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now replace the derivatives in (1.1)-(1.3) and (1.5) by the approximations,

(2.1) [[partial derivative].sub.[tau]] [c.sup.i.sub.n] = [[partial derivative].sub.hh][c.sup.i.sub.n], (n, i) [member of] [D.sub.N],

(2.2) [[partial derivative].sub.h][c.sup.I-1.sub.n] = G([c.sup.I.sub.n]), 0 [less than or equal to] n [less than or equal to] N,

(2.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

(2.4) [c.sup.i.sub.0] = [phi](ih), i = K, ..., I.

This is a system of algebraic equations for the unknowns [c.sup.i.sub.n], (n, i) [member of] [[bar.D].sub.N]. The sequence [k.sub.n] and the steps h, [[tau].sub.n] are given.

3. The maximum principle for the lattice problem. Here we prove a few properties of the solution to the lattice problem (2.1)-(2.4). The technique is similar to the maximum principle for parabolic PDEs.

Theorem 3.1. (Maximum principle). Let the lattice function [c.sup.i.sub.n] be defined in [[bar.D].sub.N], satisfy the system (2.1) in 5N, and its maximum or minimum be achieved at the node ([n.sup.*], [i.sup.*]). Then either [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Proof. Assume that the maximum (strict or not) is achieved at an inner node [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Let us use (2.1), its left-hand side is obviously nonnegative, while the right-hand side is nonpositive. Thus [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Due to the fact that ([n.sup.*], [i.sup.*]) is the maximum node, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The same argument applied to the nodes [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Continuing the argument, we get [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the function is constant on the whole layer. As [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the same maximum is also obtained on the layer [n.sup.*] - 1 and so the function is also constant on it. Continuing, we prove that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The proof for the minimum is similar.

The theorem says that either [c.sup.i.sub.n] = const in [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] or the node ([n.sup.*], [i.sup.*]) [member of] [[bar.D].sub.N] \ [D.sub.N].

Corollary 3.2. Let the lattice function [c.sup.i.sub.n] be defined in [[bar.D].sub.N] and satisfy the system (2.1), (2.4), its maximum or minimum (strict or not) be achieved at a node ([n.sup.*], [i.sup.*]), and h be small enough. Then either [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Proof. Due to (2.4) and Proposition 1.1, the function [phi](x) is not constant; then [c.sup.i.sub.0] = [phi](ih) is also not constant provided that h is small enough. Thus [c.sup.i.sub.n] is not constant in [[??].sub.n] for any n. Theorem 3.1 provides the rest.

Theorem 3.3. Let the lattice function [c.sup.i.sub.n] be defined in [[bar.D].sub.N] and satisfy (2.1)-(2.4). Then [c.sup.i.sub.n] < A in [[bar.D].sub.N], the number A > 0 is such that G(A) = 0.

Proof. It is sufficient to show the inequality for the maximum. The constant A is finite because G([infinity]) < 0. As G([phi]) > 0 and G is decreasing, [phi](x) < A. The maximal value of [c.sup.i.sub.n] cannot be achieved at a node [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The initial distribution [phi](x) > 1 and thus [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] cannot be the maximal value. Assume that the maximal value of [c.supl.i.sub.n] is achieved at a node ([n.sup.*], I). The right-hand side of (2.2) is negative if [c.sup.I.sub.n] > A; therefore [[partial derivative].sub.h][c.sup.I-1.sub.n] < 0. This means that [c.sup.I.sub.n] < [c.sup.I-1.sub.n] and thus [c.sup.I.sub.n] cannot be maximal. If [c.sup.I.sub.n] = A then the maximal value is achieved also at the inner node (n, I - 1). This is impossible due to Corollary 3.2.

Note that the bound A depends only on G(x) and is independent of h, [[tau].sub.n].

Theorem 3.4. Let the lattice function [c.sup.i.sub.n] be defined in [[bar.D].sub.N] and satisfy (2.1)-(2.4). Then [c.sup.i.sub.n] > 1 in [[bar.D].sub.N].

Proof. It is sufficient to show the inequality for the minimum. For n = 0 it is given that [phi](x) > 1. The minimal value C [less than or equal to] 1 cannot be achieved at a node ([n.sup.*],I). Assume the contrary; (2.2) yields [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (because G(C) > 0) and thus [right arrow]. So [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] cannot be the minimum. Assume that the minimal value C < 1 is achieved at a node [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. From (2.3) we get [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and thus a contradiction. Let us consider the case of the minimal value [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. From (2.3) it follows that the minimum is also achieved at the inner node: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. This is impossible due to Corollary 3.2.

Note that the lattice solution [c.sup.i.sub.n] > 1 in [[bar.D].sub.N] and thus it is positive.

Theorem 3.5. (Uniqueness) Let [c.sulp.i.sub.n] and [w.sup.i.sub.n] be two solutions to system (2.1)-(2.4) in [[bar.D].sub.N]. Then [c.sup.i.sub.n] = [w.sup.i.sub.n] in [[bar.D].sub.N].

Proof. Consider the lattice function [w.sup.i.sub.n] = [c.sup.i.sub.n]. We need to prove that [w.supl.i.sub.n] = 0 in [[bar.D].sub.N]. Suppose the contrary: let [w.sup.i.sub.n] > 0 at some node (n, i).

The function [w.sup.i.sub.n] satisfies (2.1) due to linearity. Therefore, Corollary 3.2 guarantees that the positive maximum is achieved on the boundary of [[bar.D].sub.N]. But [w.sup.i.sub.0] = 0, so the maximum cannot be achieved if n = 0.

Note that [w.sup.i.sub.n] satisfies the condition

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

If the maximum is achieved at a node (n, [k.sub.n]), then the left-hand side is nonpositive while the right-hand side is greater than zero (remember that g' > 0). The contradiction means that a maximum is impossible.

Finally, [w.sup.i.sub.n] satisfies the condition

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

If the maximum is achieved at the node (n, I), then the left-hand side is nonnegative, while the right-hand side is less than zero (G' < 0). A maximum is impossible.

The contradiction implies that no positive maximum of [w.sup.i.sub.n] can be achieved in [[bar.D].sub.N]. Thus [c.sup.i.sub.n] = [u.sup.i.sub.n] in [[bar.D].sub.N].

Let us consider the lattice set [[bar.D]'.sub.N]: if (n, i) [member of] [[bar.D]'.sub.N] then 0 [less than or equal to] n [less than or equal to] N, [k.sub.n] [less than or equal to] i [less than or equal to] I - 1. In other words, the set contains one node less at each layer compared to [[bar.D].sub.N]. Also we will consider the subset [D'.sub.N] with the nodes (n, i), 0 < n < N, [k.sub.n] + 1 [less than or equal to] i [less than or equal to] I - 2. The set [[??]'.sub.N] is equal to [[bar.D]'.sub.N] without the corner nodes (0, K) and (0, I - 1).

Let [c.sup.i.sub.n] satisfy the system (2.1)-(2.4) in 5N. The lattice derivative [[partial derivative].sub.h][c.sup.i.sub.n] is defined as a lattice function in [[bar.D]'.sub.N] and satisfies (2.1) in [D'.sub.N]. Thus it obeys Theorem 3.1: if its minimum or maximum is achieved at the node ([n.sup.*], [i.sup.*]), then either [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Theorem 3.6. Let [c.sup.i.sub.n] satisfy (2.1)-(2.4) in [[bar.D].sub.N]. The lattice derivative [[partial derivative].sub.h][c.sup.i.sub.n] is bounded: [absolute value of [[partial derivative].sub.h][c.sup.i.sub.n] < B in [bar.D]'.sub.N]. The constant B does not depend on h and [[tau].sub.n]. If additionally [phi]' (x) > 0 on [[[rho].sub.0], L] then [[partial derivative].sub.h][c.sup.i.sub.n] > 0 in [[bar.D]'.sub.N].

Proof. Extremal values of [[partial derivative].sub.h][c.sup.i.sub.n] cannot be achieved in [D'.sub.N] due to the maximum principle. Therefore, we only need to check if the values are bounded at the nodes (0, i), (n, I - 1), and (n, [k.sub.n]). Due to (2.4) [[partial derivative].sub.h][c.sup.i.sub.0] = [phi]'([x.sub.i] + [theta]), 0 [less than or equal to] [theta] [less than or equal to] h, is bounded for all i = K, ..., I - 1. From (2.2) and Theorems 3.3 and 3.4, we have [[partial derivative].sub.h][c.sup.I-1.sub.n] = G([c.sup.I.sub.n]) [member of] (0, G(1)). Finally, (2.3) and Theorems 3.3 and 3.4 imply [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Let us construct the lattice subset [[bar.D]".sub.N] [subset] [[bar.D]'.sub.N] in the following way: a node (n, i) [member of] [[bar.D].sub.N] belongs to [[bar.D]".sub.N] if the node (n - 1, i) belongs to [[bar.D].sub.N]. One can see that in [[bar.D]".sub.N] the index n > 0 (i.e., the layer n = 0 is not contained); also the nodes (n, [k.sub.n]) are not contained if [K.sub.n] = [K.sub.n-1]. Therefore this subset is generated by some sequence [K".sub.n], n [greater than or equal to] 1. Denote [k".sub.n] = [K".sub.n] - n. Let us consider also the subsets [[bar.D].sub.N] and [[??]".sub.N] similarly to [D.sub.N] and [[??].sub.N].

Let [c.sup.i.sub.n] satisfy the system (2.1)-(2.4) in [[bar.D].sub.N]. The lattice derivative [[partial derivative].sub.[tau]][c.sup.i.sub.n] is defined as a lattice function on [[bar.D]".sub.N] and satisfies (2.1) in [D".sub.N]. Thus it obeys Theorem 3.1; if its minimum or maximum is achieved at the node ([n.sup.*], [i.sup.*]) then either [n.sup.*] = 1 or [i.sup.*] = I or [i.sup.*] = [".sub.n], or it is constant in [[??].sub.N].

Theorem 3.7. Let [c.sup.i.sub.n] satisfy (2.1)-(2.4) in [[bar.D].sub.N]. The lattice derivative [[partial derivative].sub.[tau]][c.sup.i.sub.n] is bounded: [absolute value of [[partial derivative].sub.[tau]][c.sup.i.sub.n]]< Z in [[bar.D]".sub.N]. The constant Z does not depend on h and [[tau].sub.n].

Proof. Extremal values of [[partial derivative].sub.h] [c.sup.i.sub.n] cannot be achieved in [D".sub.N] due to the maximum principle. Therefore, we only need to check if the values of [[partial derivative].sub.[tau]] [c.sup.i.sub.n] are bounded at the nodes (n, I), (n, [k".sub.n]), and (0,i).

Consider two nodes (n + 1, I) and (n, I), n [greater than or equal to] 0. Due to (2.2),

[c.sup.I.sub.n+1] - [c.sup.I-1.sub.n+1] = hG'([c.sup.I.sub.n+1]), [c.sup.I.sub.n] - [c.sup.I-1.sub.n] = hG([c.sup.I.sub.n].

Substract the second expression from the first one and divide by [[tau].sub.n],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

As G' < 0, the positive maximum and negative minimum of the function [[partial derivative].sub.[tau]][c.sup.i.sub.n] cannot be achieved for i = I, otherwise the contradiction appears.

Now consider two nodes (n + 1, [k".sub.n+1]) and (n, [k".sub.n]), n [greater than or equal to] 0. There can be two cases. First, let us study the case [K".sub.n+1] = [K".sub.n] + 1. Due to (2.3),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Substract the second expression from the first one and divide by [T.sub.[tau]],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

As g' > 0, the positive maximum and negative minimum of the function [[partial derivative [tau]] [c.sup.i.sub.n] cannot be achieved at the node (n, [k".sub.n]) if [K".sub.n] = [K".sub.n-1] + 1 .

Second, we study the case [".sub.n+1] = [K".sub.n]. Due to (2.3)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Substract the first expression from the second one,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now let us transform the first bracket on the left-hand side,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and the second bracket on the left-hand side,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The bracket on the right-hand side equals the first one on the left-hand side. Therefore,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The big bracket equals [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (we have used (2.1)). Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Then due to boundedness, [[partial derivative].sub.h][c.sup.i.sub.n] [less than or equal to] B,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now we see that if [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], is the maximal positive value, then it is less than [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Otherwise, the bracket in the right-hand side is negative and thus

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A similar argument shows that the minimal possible value is also bounded.

Now let us study the case n =1. Let [[epsilon].sup.i] = [c.sup.i.sub.1] - [c.sup.i.sub.0] and consider (2.1) for n = 1, K [less than or equal to] i [less than or equal to] I - 1:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

But [c.sup.i-1.sub.0] - 2[c.sup.i.sub.0] + [c.sup.i+1.sub.0] = [h.sup.2][phi]"([x.sub.i]) + o([h.sup.2]) and thus the first term in the right-hand side is bounded. Besides [[partial derivative].sub.[tau]][c.sup.i.sub.1] = [[epsilon].sup.i]/[[tau].sub.0]. Therefore,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Here R(i) is bounded independently on h and [[tau].sub.0]. Let the function [[partial derivative].sub.[tau]][c.sup.i.sub.n] reach a positive maximum at the node (1, [i.sup.*]); then [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] also reaches the positive maximum. But the first term in the right-hand side is negative; so [[tau].sup.-1.sub.0] [[epsilon].supl.i] [less than or equal to] R (i) and thus is bounded. A similar argument applies to the negative minimum.

Let us sum this all up. Minimal and maximal values of the lattice function [[partial derivative].sub.[tau]] [c.sup.i.sub.n] in [[bar.D]".sub.N] are achieved at the nodes ([n.sup.*], [i.sup.*]) for either [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. But these extremal values are bounded independently on the steps; therefore there exists Z such that [absolute value of [[partial derivative].sub.[tau]][c.sup.i.sub.n]] [less than or equal to] Z in [[bar.D]".sub.N].

Corollary 3.8. The second lattice derivative [[partial derivative].sub.[hh]] [c.sup.i.sub.n] is bounded independently of the steps at nodes where it is defined.

Corollary 3.9. If the initial distribution has a positive second derivative [phi]"(x) > 0 and the step h is sufficiently small, then [[partial derivative].sub.[tau]] [c.sup.i.sub.n] > 0 in [[bar.D]".sub.N].

Corollary 3.10. Let [phi]'(x) > 0 and [phi]"(x) > 0 and the step h be sufficiently small. Then 0 < [[partial derivative].sub.[h]] [c.sup.i.sub.n] < B in [[bar.D]'.sub.N] and B = [[partial derivative].sub.h][c.sup.I-1.sub.0].

Proof. As [[partial derivative].sub.[tau]][c.sup.i.sub.n] > 0, then also [[partial derivative].sub.[hh]] [c.sup.i.sub.n] > 0 at nodes where it is defined. This implies, in particular, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and therefore the maximal value of [[partial derivative].sub.h][c.sup.i.sub.n] cannot be achieved at a node (n, [k.sub.n]). As [[partial derivative].sub.[tau]] [c.sup.i.sub.n] > 0, [c.sup.i.sub.n] increases; due to monotonicity of the function G the right-hand side of (1.2) decreases. Thus if the maximal value of [[partial derivative].sub.h][c.sup.i.sub.n] is achieved at a node (n, I), then n = 0. The only possibility for the maximum is the node (0, I - 1).

4. Choosing the time step. Up to this point we considered the given lattice; the spatial step h, the time steps [[tau].sub.n], the sequence [K.sub.n] (it generates the subset [[bar.D].sub.N]), and the time T were fixed. But [[tau].sub.n] and [K.sub.n] have to be determined. Some equations must be added to the system (2.1)-(2.4) to find not only unknown [c.sup.i.sub.n] but also unknown [T.sub.n] and [K.sub.n].

Let us consider the lattice analogue of the Stefan condition (1.4),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

We have already proved that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and thus the shift of the free boundary is negative: [[rho].sub.n+1] - [[rho].sub.n] < 0. In order to make the absolute value of this shift equal h (up to the error of approximation), we need to choose the special [T.sub.n],

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

However, it can turn out that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is close to unity, so [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is too small and then [T.sub.n] is too large. If [[tau].sub.n] does not tend to zero as h [right arrow] 0, then there is no approximation.

From (1.4) it follows that the boundary moves slowly when the concentration near the boundary is close to unity: [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Therefore let us assume that for small [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], the boundary does not move. Let us consider [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] small if

(4.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The left-hand side is positive. Let the step h obey the inequality h [less than or equal to] 1 . The reason is to guarantee that the initial distribution is not small.

Summing this up, we choose the steps [[tau].sub.n] as follows,

(4.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

Here [g.sup.-1] is for the inverse function, [LAMBDA] = const is independent of h. Note that bounds for the time step [[tau].sub.n] > h/g(A), [[tau].sub.n] < [bar.[tau]] = 0([square root of h]) hold.

In the first, case the boundary shifts left one step h during a time step; in the second case, it remains motionless. In other words, in the first case we add a node to the new layer, while in the second case we do not. Therefore,

(4.3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The system (2.1)-(2.4), (4.2), and (4.3) is called the system (*).

Proposition 4.1. Let [c.sup.i.sub.n] satisfy the system (*) in 551. Then [[partial derivative].sub.h][c.sup.i.sub.n] = const in [[bar.D].sub.N].

Proof. Due to h < 1 the initial value [phi]([[rho].sub.0]) is not small and thus [K.sub.1] = [K.sub.0] = K. Thus the layer n = 1 has one node more compared to the layer n = 0. Assume that the lattice derivative [[partial derivative].sub.h][c.sup.i.sub.n] is constant in [[bar.D]'.sub.1]. Then due to (2.4) [[partial derivative].sub.h][c.sup.i.sub.n] = [phi]'([[rho].sub.0]) in [[bar.D].sub.1]. In particular, [[partial derivative].sub.h][c.sup.K.sub.0] = [[partial derivative].sup.K-1.sub.1] and from (2.3) we have g([c.sup.K.sub.0]) = [phi]'([c.sup.K-1.sub.1]). As g(c) is monotonic, this means [c.sup.K.sub.0] = [c.sup.K-1.sub.1]. Thus [c.sup.K.sub.1] = [c.sup.K.sub.0] + [phi]'([[rho].sub.0])h. In the similar way, [c.sup.K+i.sub.1] = [c.sup.K.sub.0] + [phi]'([[rho].sub.0])(i + 1)h. Thus [c.sup.K.sub.1] = [c.sup.K.sub.0] + [phi]'([[rho].sub.0])(N - K + 1)h. As the derivative is constant, [c.sup.N.sub.0] = [c.sup.K.sub.0] + [phi]'([[rho].sub.0])(N - K)h [not equal to] [c.sup.N.sub.1]. However, due to (2.2) G([c.sup.N.sub.0]) = G([c.sup.N.sub.1]) and, thus, [c.sup.N.sub.0] = [c.sup.N.sub.1] because G(c) is monotonic. This contradiction shows that [[partial derivative].sub.h][c.sup.i.sub.n] cannot be constant in [[??]'.sub.1] and thus in [[bar.D]'.sub.N] for any N [greater than or equal to] 1.

5. Solving the system. We are going to present an algorithm for the solution of the system (*). Thus we will prove that it has a solution. A solution is a set [T.sub.n] for 0 [less than or equal to] n [less than or equal to] N - 1, [K.sub.n] for 0 [less than or equal to] n [less than or equal to] N, [c.sup.i.sub.n] for 0 [less than or equal to] n [less than or equal to] N, [K.sub.n] [less than or equal to] i [less than or equal to] I for given natural M; h < 1, I, K, N, and T are uniquely determined.

Let us develop a sweep method to obtain [c.sup.i.sub.n] when the [c.sup.i.sub.n-1] are already known. Denote [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Express all [c.sup.i.sub.n], i = [k.sub.n] + 1, ..., I - 1 via [c.supl.i+1.sub.n] and X linearly with unknown coefficients,

(5.1) [c.sup.i.sub.n] = [a.sub.i][c.sup.i+1.sub.n] + [b.sub.i] + [d.sub.i] X.

Substitute (5.1) into (2.1) instead of [c.sup.i-1.sub.n] 1 for i = [K.sub.n] + 2, ..., I - 1 . The result is

[c.sup.i.sub.n] = [c.sup.i+1.sub.n] + [S.sub.n][c.sup.i.sub.n-1] + [b.sub.i-1] + [d.sub.i-1]X/2 + [S.sub.n] - [a.sub.i - 1].

Here [S.sub.n] = [h.sup.2]/[[tau].sub.n-1]. Then

(5.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This are the recurrent sequences. To get the initial values consider (2.1) for i = [K.sub.n] + 1 ,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now express [c.sup.i.sub.n], i = [k.sub.n] + 1 , ..., I - 1 via [c.sup.i-1.sub.n] similarly,

(5.3) [c.sup.i.sub.n] = [A.sub.i][c.sup.i-1.sub.n] + [B.sub.i] + [D.sub.i]Y.

Analogously,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

with initial values [A.sub.I-1] = [D.sub.I-1] = [(2 + [S.sub.n]).sup.-1], [B.sub.I-1] = [A.sub.I-1][S.sub.n] [c.sup.I-1.sub.n-1].

Proposition 5.1. The following inequalities holdfor [K.sub.n] + 1 [less than or equal to] i [less than or equal to] I - 1 : [a.sub.i] [member of] (0,1), [a.sub.i+1] > [a.sub.i], [d.sub.i] [member of] (0, 0.5), [d.sub.i+1] < [d.sub.i], [b.sub.i] > 0, [a.sub.i] < 1 - [[LAMBDA].sub.1] h 3/4 for some constant 0 < [[LAMBDA].sub.1] < 1. The same inequalities hold for [A.sub.i], [B.sub.i], [D.sub.i].

Proof. The first five inequalities are proven by induction; obviously they are true for i = [k.sub.n] + 1; if they hold for some i then (5.2) shows that they also do for i + 1.

The sequence [a.sub.i] increases and is bounded; therefore it has an upper bound. To calculate it find the fixed point a of the function [a.sub.i]([a.sub.i-1]), i.e., solve the equation [a.sup.2]-(2+[S.sub.n])a+1 = 0. It has a real root a from (0,1). If [a.sub.i] = a then [a.sub.i+1] = a. As [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Let us estimate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

This finishes the proof. The proof for [A..sub.i], [B.sub.i], [D.sub.i] is similar.

Substitute (5.1) for i = I - 1 and (5.3) for i = [k.sub.n] + 1 into (2.2) and (2.3), respectively,

(5.4) Y(1 - [a.bu.I-1]) - hG(Y) = [b.sub.I-1] + [d.sub.I-1] X,

(5.5) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This is a system of two non-linear equations with two unknowns, X and Y.

We have proved that a solution to the lattice problem (2.1)-(2.4) with any steps h and [[tau].sub.n] is greater than unity and bounded independently of the step sizes. Therefore we can suppose without loss of generality that the function g(c) grows faster than c: let g(c)/c [right arrow] [infinity] as c - [infinity]. For technical purposes we may need G(c) for negative c; suppose that it remains continuous and monotonic and thus G(c) > 0 for c < 0.

Theorem 5.2. The system of equations (5.4) and (5.5) has a solution; if h is sufficiently small, then the solution is unique.

Proof. Let us express Y as a function of X defined by (5.5). We see that Y(0) < 0 because g(1) = 0 and increases (and thus g(0) [less than or equal to] 0), [B.sub.i] > 0, [D.sub.i] > 0. On the other hand, Y(+[infinity]) = +[infinity] because 1 - [A.sub.i] > 0. Moreover, Y(X) grows faster than X.

Now substitute Y(X) into (5.4) and consider the continuous function

F(X) = Y(X)(1 - [a.sub.I-1]) - hG(Y(X)) - [b.sub.I-1] - [d.sub.I-1]X.

If F([X.sup.*]) = 0, then the pair ([X.sup.*], Y([X.sup.*])) is a solution to the system of equations (5.4) and (5.5). To prove that F(X) has a zero, note that F(0) < 0 (because [a.sub.i] < 1 , G is decreasing, [b.sub.i] > 0, [d.sub.i] > 0) and F([infinity]) = [infinity] (because [a.sub.i] < 1 and Y(X) grows faster than X). Hence, there exist [X.sup.*] [greater than or equal to] 0 and [Y.sup.*] = Y([X.sup.*]), such that equations (5.4) and (5.5) are satisfied.

To prove the uniqueness, let us show that the derivative F' (X) > 0,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

To prove that F' (X) > 0, it is sufficient to show that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

These inequalities are proved similarly, so let us prove the second one,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Suppose that h is small enough so that 2(1 - [[LAMBDA].sub.1]h 3/4)2 > 1. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The function f(z) = (1 - z) 1/z is decreasing. Therefore, f(z) < /(0) = exp(-l). Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for any [alpha] > 0. Hence, for sufficiently small h, we have [d.sub.I-1] < [[LAMBDA].sub.1]h 3/4.

Note that a small h is not necessary for the uniqueness of the solution.

Here we present the algorithm:

1. Let n = 0. Define [K.sub.0] = K. For [K.sub.0] [less than or equal to] i [less than or equal to] I calculate [c.sup.i.sub.0] using (2.4).

2. Find [[tau].sub.n] from 4.2). obviously h < [[tau].sub.n] < [[bar.[tau]].sub.n].

3. Detemine [K.sub.n+1] from (4.3) and calculate [k.sub.n+1] = [K.sub.n+1] - (n + 1).

4. If [k.sub.n+1] = 0 then stop: the problem is solved because the boundary has reached zero. Assign T = [summation] [[tau].sub.n], N.

5. Increase n by one.

6. Solve equations (5.4) and (5.5) simultaneously to obtain X and Y.

7. Using (5.1) and the computed [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] calculate [c.sup.i.sub.n] one by one for i = I - 1 , I - 2, ..., [k.sub.n] + 1 . It is also possible to use (5.3) .

8. Go to step 2.

We have not proved that the determined X > 1 and Y > 1 ; but this is true because the constructed solution obeys Theorem 3.4.

Theorem 5.3. Provided that h is sufficiently small, the algorithm terminates after a finite number of steps.

Proof. Suppose that the algorithm never stops; the boundary moves only in a finite number of steps. Thus, after the finite number of steps the boundary does not move ([k.sub.n] = const for sufficiently large n). Let us see how the amount of matter changes,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The boundary does not move, so [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is small, and thus

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

It is necessary that [[partial derivative].sub.h][c.sup.I-1.sub.n] - [phi]'([[rho].sub.0])[square root of h] [right arrow] 0, otherwise the amount [I-1.summation over (i=[k.sub.n+1])] [c.sup.i.sub.n] h grows to infinity and some [c.sup.i.sub.n] become greater than A for large n. Thus if n is large,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The last inequality holds for sufficiently small h and some [epsilon] > 0. Hence, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

On the other hand, if [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] remain small, then [[partial derivative].sub.[h]] [c.sup.i.sub.n] for all i = [k.sub.n], ... ,1 - 1 become small for sufficiently large n. This follows from the maximum principle. So we can assume that [absolute value of [[partial derivative].sub.h][c.sup.i.sub.n]] < [phi]'([[rho].sub.0])[square root of h] for i = [k.sub.n], ..., I - 1. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

If h is sufficiently small, this is below any [epsilon] > 0. Thus we have the contradiction.

It is easy to show that the time T(h) = [summation] [[tau].sub.n] is also finite and, moreover, bounded for all h: T(h) [less than or equal to] [bar.T] for some [bar.T].

6. Order of approximation. Let us suppose that the classical smooth solution p(t), c(t, x) > 1 to the boundary-value problem (1.1)-(1.5) exists. We need to see how large the error is if we substitute it into the system (*).

First, we see that, for all n, [[tau].sub.n] [right arrow] 0 as h [right arrow], and thus, the rate of convergence is at least [bar.[tau]] = O([square root of h]). Consider the equation [[partial derivative].sub.[tau]][c.supl.i.sub.n] - [[partial derivative].sub.hh][c.sup.i.sub.n] = [[PSI].sup.i.sub.n]. If [c.sup.i.sub.n] is a solution to (2.1) then *[[PSI].sup.i.sub.n] = 0 in [D.sub.N]. Substitute the values of the exact solution c([t.sub.n], [x.sub.i]) at nodes of the lattice subset [[bar.D].sub.N] into the left-hand side. Then in general [[PSI].sup.i.sub.n] [not equal to] 0 is some lattice function called the discrepancy. Let [PSI] denote the maximum of the absolute value of the [[PSI].sup.i.sub.n] on [D.sub.N]. Using the Taylor expansion for the exact solution and the diffusion equation (1.1), we obtain [PSI] = 0{[h.sup.2], [bar.[tau]]) = 0([square root of h]). These are standard arguments of the theory of numerical methods, so the details are omitted. In a similar way, we substitute the exact solution into (2.2) and (2.3). The discrepancy turns out to be O(h). The initial condition (2.4) is satisfied precisely.

Let us recall how we approximated the Stefan condition (1.4),

(6.1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

and zero otherwise. Substitute the exact solution [rho](t), c(t, [rho](t)) with a discrepancy [[PSI].sub.n] and apply the Stefan condition (1.4) ,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

and zero otherwise. In both cases [[PSI].sub.n] = O ([square root of h], + [bar.[tau]]) = O ([square root of h]). Thus the total discrepancy and the order of approximation of the scheme is O ([[square root of h]).

The order of approximation cannot be made better than i/^byjust choosing the condition (4.1) in the form g^") < p'(p0)h1-e with e > 0.5. If we do so then [T.sub.n] < 5 = O(he). The discrepancy of the Stefan condition and thus the total discrepancy will become worse, O([h.sup.1-[epsilon]]).

7. Convergence of approximations. In the remainder of the paper, we assume that h is sufficiently small.

If an h is chosen, one can determine a solution using the presented algorithm. The solution is the set of T, [[tau].sub.n] for 0 [less than or equal to] n [less than or equal to] N - 1, [[bar.D].sub.N], and [c.sup.i.sub.n] in [[bar.D].sub.N]. Theorem 3. 5 guarantees that if the domains 55 n are the same for two solutions, they coincide.

Construct the piecewise linear continuous time function [[rho].sub.h](t) by connecting the points [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (linear interpolation). Being piecewise linear [[rho].sub.h] (t) has a piecewise constant derivative between the nodes; let us estimate it. The difference [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is either zero (if [K.sub.n+1] = [K.sub.n] + 1) or -h (if [K.sub.n+1] = [K.sub.n]). The difference [t.sub.n+1] - [T.sub.n] = [T.sub.n]. We know that [T.sub.n] [greater than or equal to] h/g(A). Thus the derivative between the nodes belongs to [-g(A), 0] for any h. The function [[rho].sub.h](t) [member of] [0, [[rho].sub.0]] for all h and is defined for t [member of] [0, T(h)]. Let us expand it to [0, T] being continuous and constant for t > T.

Now let us consider a sequence [h.sub.j] [right arrow] 0. The functions [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are uniformly bounded and equicontinuous; due to the Arzela-Ascoli theorem a uniformly convergent subsequence can be selected (but may be not unique). Denote this subsequence again by [h.sub.j] and its limit by [rho](t). The functions [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] are all nonnegative, uniformly bounded by [[rho].sub.0] < L, and non-increasing; thus [rho](t) is also nonnegative, bounded by [[rho].sub.0], and does not increase. Also [rho](T) = 0 for some T [member of] [0, [bar.T]. This is the continuous free boundary.

Now define the lattice function [c.sup.i.sub.n] at the nodes (n, i) [member of] [[bar.D].sub.N]: let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] at other nodes (n, i) [member of] [[bar.D].sub.N] let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The continuous approximations [[psi].sub.h] (t, x) in the rectangle N = [0, T] x [0, L] are constructed in the following way. Let [[psi].sub.h] ([t.sub.n], [x.sub.i]) = [c.sup.i.sub.n] at the nodes (n, i). For each n we define [[psi].sub.h] ([t.sub.n], x) to be a polynomial of the fifth order [P.sub.n,i] (x) on each segment [[x.sub.i], [x.sub.i+1]] for i = [k.sub.n] ..., I - 2. The six coefficients of the polynomial are obtained from six equations: [P.sub.n,i] ([x.sub.i]) [c.sup.i.sub.n] = [P.sub.n,i] ([x.sub.i+1]) = [c.sup.i+1.sub.n] [P.sub.n,i] ([x.sub.i])' = [[partial derivative].sub.h] [c.sup.i.sub.n], [P.sub.n,i] ([x.sub.i+1])' = [[partial derivative].sub.h] [c.sup.i+1.sub.n], [P.sub.n,i] ([x.sub.i])" = [[partial derivative].sub.[hh]] [c.sup.i.sub.n], [P.sub.n,i] ([x.sub.i+1])" = [[partial derivative].sub.hh] [c.sup.I.sub.n]. On the segments [[x.sub.I-1], [x.sub.I]] we have only four equations because [[partial derivative].sub.h][c.sup.I-1.sub.n] and [[partial derivative].sub.[hh]] [c.sup.I.sub.n] are not defined. Let us add the equations [P.sub.n,I-1] ([x.sub.I])' = [[partial derivative].sub.h][c.sup.I-1.sub.n] and [P.sub.n,I-1] ([x.sub.I])" = 0. The equations for the coefficients are linear with a regular matrix (its determinant can be easily calculated). Finally we define [[psi].sub.h] (t,x) being linear with respect to t for t [member of] [[t.sub.n], [t.sub.n+1]], x [member of] [[x.sub.i], [x.sub.i+1]],

[[psi].sub.h] (t,x) = [P.sub.n,i](x) + t - [t.sub.n] / [[tau].sub.n] x ([P.sub.n+1,i] (x) - [P.sub.n,i] (x)).

Obviously [partial derivative].sub.x] [[psi].sub.h] (t,x) and [[partial derivative].sub.xx] [[psi].sub.h] (t,x) are continuous and uniformly bounded, and [[partial derivative].sub.t] [[psi].sub.h] (t, x) and [[partial derivative].sub.xt] [[psi].sub.h] (t, x) exist between the nodes and are also uniformly bounded (see Theorems 3.6 and 3.7 and Corollary 3.8); thus

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Due to the fact that [[psi].sub.h] (t, x), [[partial derivative].sub.x] [[psi].sub.h] (t, x), [[partial derivative].sub.x] [[psi].sub.h] (t, x), and [[partial derivative].sub.xx] [[psi].sub.h] (t, x) are uniformly bounded at the nodes, the functions [[psi].sub.h] (t, x) and [[partial derivative].sub.x] [[psi].sub.h] (t, x) are uniformly bounded and equicontinuous (even equi-Lipschitz) in N.

Now take the sequence [h.sub.j][right arrow] 0 chosen above and consider the corresponding sequence [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. These functions are uniformly bounded and equicontinuous; by the Arzela-Ascoli theorem a uniformly convergent subsequence can be chosen. Let its limit be [psi] (t, x). It is continuous. Denote the corresponding subsequence of [h.sub.j] again by [h.sub.j]. It is clear that [psi] (t, x) [greater than or equal to] 1, is bounded, and satisfies (1.5). Moreover, the derivatives [[partial derivative].sub.x] [[psi].sub.h] (t, x) are themselves uniformly bounded and equicontinuous. Applying the Arzela-Ascoli theorem, we learn that continuous [[partial derivative].sub.x] [psi] (t, x) exists and term-by-term differentiation is allowed.

Proposition 7.1. The function [psi] (t, x) is Lipschitz continuous in N.

Proof. Let [??] = (t, x) and bound

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Here [epsilon] > 0 is any number. We have used Lipschitz continuity of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] and their uniform convergence (assuming that j is large enough). Passing to the limit as [epsilon] [right arrow] 0 finishes the proof. []

Theorem 7.2. The free boundary p(t) has negative Lipschitz continuous derivative for t [member of] [0, T] and the Stefan condition (1.4) holds.

Proof. Rewrite the Stefan condition (1.4) in integral form,

(7.1) [rho](t) = [[rho].sub.0] [[integral].sup.t.sub.0] g(c([xi]), [rho]([xi]))) / c ([xi], [rho]([xi])) d [xi].

Substitute the approximations [[psi].sub.h] (t, x) and [[rho].sub.h] (t) into (7.1),

(7.2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Here [[PSI].sub.h] is the discrepancy. We need to prove that [[PSI].sub.h] [right arrow] 0 as h [right arrow] 0.

First, let us study [[rho].sub.h] (t). Let [T.sub.n] [less than or equal to] t [less than or equal to] [t.sub.n+1]. We defined [[rho].sub.h] ([t.sub.n]) by linear interpolation, so that

[[rho].sub.h] ([t.sub.n]) = [[rho].sub.0] + [n.summation over (m=1)] [[partial derivative].sub.[tau] [[rho].sub.m] [[tau].sub.m-1].

Let [[mu].sub.n] be the set of natural m [less than or equal to] n, such that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] is not small, i.e., [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] [greater than or equal to] [phi]' ([[rho].sub.0]) [square root of h]. Let [[bar.[mu]].sub.n] be its complement, i.e., m = 1, ..., n that are not in [[mu].sub.n]. Now apply (6.1),

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The sum over m [member of] [[bar.[mu]].sub.n] has positive terms and each is at most [phi] ([[rho].sub.0]) [square root of h]; thus the sum is 0([square root of h]). Note that [[rho].sub.h] (t) = [[rho].sub.h] ([t.sub.n]) + O ([[tau].sub.n]) = [[rho].sub.h] ([t.sub.n]) + 0 ([square root of h)].

Let us consider

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

We have the continuous functions under the integrals; so that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Let us rewrite o ([[tau].sub.m-1]) = [[tau].sub.m-1] [[omega].sub.m-1], [[omega].sub.m] [right arrow] 0 as h [right arrow] 0. Then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The last sum tends to zero as h [right arrow] 0; denote it by W. Note that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Let us substitute the obtained expressions into (7.2),

O ([square root of h]) = O ([square root of h]) + W + [[PSI].sub.h].

Thus [[PSI].sub.h] [right arrow] 0 as h [right arrow] 0 and the uniform limit of (7.2) exists. Therefore, the obtained continuous solution [rho] (t), [psi] (t, x) satisfies the integral Stefan equation (7.1). Its left-hand side p(t) is not only continuous, but also has the continuous derivative for t [member of] (0, T) . Differentiating by t we see that the solution [rho] (t), [psi] (t, x) also satisfies the original Stefan condition (1.4). Moreover, the derivative [??] is Lipschitz because [psi] is and [??] = g ([psi]) / [psi], where the right-hand side is smooth. []

Corollary 7.3. The free boundary [rho] (t) has the inverse function [[rho].sup.-1] (x).

8. Weak solution to the Dirichlet boundary-value problem. Let the set [Y.sub.[rho]] (T) [subset] [R.sup.2] contain all (t, x) such that t [member of] (0, T), x [member of] (p(t), L), let [[bar.Y].sub.[rho]] (T) be its closure. We are going to prove that [psi] (t, x) is the weak solution to problem (1.1)-(1.5) in [Y.sub.[rho]] (T). Consider the Dirichlet boundary-value problem,

(8.1) [[partial derivative].sub.t] c(t, x) = [[partial derivative].sub.xx] c(t, x), (t, x) [member of] [Y.sub.[[rho]] ( T),

(8.2) c(t, L) = [psi] (t,L), c(t, [rho] (t)) = [psi] (t, [rho] (t)), t [member of] [0,T],

(8.3) c(0, x) = [phi] (x), x [member of] [[[rho].sub.0], L], 0 < [[rho].sub.0] < L.

We are going to define a weak solution to this problem (similar to [7, 9, 10]).

Definition 8.1. A weak solution of the problem is a continuous function c(t, x) that has the weak derivative [[partial derivative].sub.x] c in [[bar.Y].sub.[rho]] (T) and satisfies conditions (8.2) and (8.3) and the integral identity

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

for each continuous in [[bar.sub.[rho]] (T) function v(t, x) with the weak derivative [[partial derivative].sub.x] v; the function v is such that v(0, x) = v(T, x) = 0, v(t, p(t)) = v(t, L) = 0.

Let us prove that the constructed function [psi] is the only weak solution.

Proposition 8.2. The constructed solution [psi] (t, x) to the problem (1.1)-(1.5) is a weak solution to the problem (8.1)-(8.3).

Proof. Boundary and initial conditions (8.2) and (8.3) are obviously satisfied. The identity is proved by considering it for the [[psi].sub.h] up to the small errors and passing to the limit. []

Theorem 8.3. The weak solution to (8.1)-(8.3) is unique.

Proof. Suppose [10] that there are two solutions [c.sub.1] (t, x) and [c.sub.2] (t, x); their difference u(t, x) = [c.sub.1] - [c.sub.2] also satisfies the identity; besides it satisfies the homogenous boundary conditions: u(0, x) = 0, u(t, [rho] (t)) = u(t, L) = 0. Consider the function v (t,x) such that v(T, x) = 0, [[partial derivative].sub.t] v (t, x) = -u(t, x). Substitute v(t, x) into the identity,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Let us transform the right-hand side,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Adding the identities with these two forms on the right-hand side, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The two integrals with ^ as the integration variable are independent of t. Thus,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This can hold only if u = 0 in [[bar.Y].sub.[rho]] (T) and so the uniqueness is proved. []

As the solution to (8.1)-(8.3) is unique, there are no other solutions beside the function [psi] (t, x) constructed above .

9. Regularity. Let us show that the constructed solution [psi] (t, x) is actually classical, i.e., smooth, by applying results from [11]. First, we transform the domain to a rectangle by changing the spatial variable x in (8.1)-(8.3) in the following way: x = [rho] (t)+ y x (L - [rho] (t)), u(t,y) = c(t,x), y [member of] [0,1], [rho] the obtained free boundary. The problem in the new variables looks like

(9.1) [[partial derivative].sub.t] (t,y) = a(t,y) [[partial derivative].sub.yy] u(t,y) + b(t,y) [[partial derivative].sub.y] u(t,y), (t,y) [member of] D,

(9.2) u(t,y) = [psi] (t,y), (t, y) [member of] [bar.B] x S,

where a = [(L - [rho]).sup.-2], b = (1 - y) [(L - [rho]).sup.-1] [??], D = [bar.B] x [0,T] is a cylinder, B = (0,1), [bar.B] = [0,1], the boundary [partial derivative] B of B consists of two points, x = 0 and x = 1, and the boundary manifold S consists of two parts, [S.sub.0] = (0, T) x {x = 0} and [S.sub.1] = (0, T) x {x = 1}. The function [psi] is continuous in [[bar.Y].sub.[rho]] (T) and thus on [bar.B] x S also. Note that a does not depend on y. The notation is from [11].

Let us check the requirements of [11, Corollary 2, [section] 4, Chapter III]. Condition [bar.A] from [11] holds: the coefficients a, b are uniformly Holder continuous in 5; indeed, even Lipschitz continuous with the same Lipschitz constant. This follows from the fact that [??] is Lipschitz continuous. Condition B from [11] holds trivially: a(t, y) [greater than or equal to] [L.sup.-2] > 0. The other requirements of [11] also hold trivially.

From [11, Corollary 2, [section] 4, Chapter III], the existence of the unique classical solution u to the boundary-value problem (9.1)-(9.2) follows. The inverse change of variables c(t, x) = u(t, y(x)) provides us the classical solution of the problem (8.1)-(8.3) in [[bar.Y].sub.[rho]] (T). Being a classical solution, c(t, x) is then a weak solution and thus coincides with the weak solution [psi] (t, x) (because the weak solution is unique). We have proved that [psi] (t, x) satisfies the boundary Stefan condition (1.4). The boundary conditions (1.2) and (1.3) also hold; the reason is that they obviously hold at the nodes of the lattice and the approximations have continuous derivatives with respect to x. Thus between the nodes the boundary conditions hold up to the error O ([square root of d [x.sup.2]] + d [t.sup.2]]), where [square root of d [x.sup.2]] + d [t.sup.2]] is the distance to the nearest boundary node. Passing to the uniform limit as h - 0 finishes the proof. Therefore, we have proven the following theorem.

Theorem 9.1. The constructed pair of the functions [rho](t), [psi](t, x) is the classical solution to the problem (1.1)-(1.5).

The presented numerical method (algorithm) converges to a solution of the boundary-value problem (1.1)-(1.5) and can be used for solving the problem.

10. Numerical example. We illustrate the suggested method by a numerical example. The model is similar to that of [1]. Assume that the right-hand sides of the boundary conditions (1.2), (1.3), and the initial condition (1.5) are

G(c) = P - [beta] [c.sup.2], g(c) = [kappa]c(c - 1), [psi](x) = [[gamma].sub.1] x + [[gamma].sub.2].

Here P is the sorption flux density, which is constant provided that the temperature and the pressure are constant; the desorption flux density [beta] [c.sup.2] depends on the concentration c, the square law follows from the fact that two hydrogen atoms form a molecule [H.sub.2], [beta] < P is the desorption constant. The difference between sorption and desorption is the diffusion flux density, i.e., the right-hand side of (1.2). Hydride formation is described by the second formula, [kappa] is constant (this boundary condition differs from the appropriate one from [1], where we assumed c(t, [rho]) = 1). The initial distribution [psi] (x) is stationary for the diffusion equation, i.e., is its time-independent solution. The constants [[gamma].sub.i] are uniquely determined from the boundary conditions.

All assumptions hold, only g'(c) is negative for c < 0. 5. But this in not important because we have proved that c [greater than or equal to] 1 , so we can change the condition for c < 1 without any influence on the results. After the time T, i.e., when hydriding is over, we have the problem of saturation with fixed boundaries and nonlinear boundary condition (1.2). Another condition follows from the symmetry and looks like [[partial derivative].sub.x] c(t, 0) = 0, thus g(c) = 0. All our results hold for this problem also. The dimensionless parameters were L =1, [[rho].sub.0] = 0.63, P = 1.21 x [10.sup.5] , [beta] = 1.24 x [10.sup.4,] [kappa] = 50. To return to usual units it is enough to know L = 7 x [10.sup.-5] cm, the stoichiometric concentration in uranium 8.29 x [10.sup.22] atoms per [cm.sup.3], and the hydrogen diffusivity in U [H.sub.3] d = 7 x [10.sup.-13] [cm.sup.2]/s. The figures are in the usual units.

Figure 10.1 shows the flux densities at the boundaries; the upper curve is the flux density at x = L, the other is the flux density at the free boundary. It vanishes instantly when hydriding is over. Figure 10.2 displays the concentration at the free boundary; one can see that the assumption from [ 1] (that c(t, p) [approximately equal to] const while hydriding is not over) is reasonable.

[FIGURE 10.1 OMITTED]

[FIGURE 10.2 OMITTED]

In Figure 10.3 we present the concentration profiles at regularly distributed time instants. One can estimate the free boundary's velocity: it is almost constant. Also the concentration at x = L is almost constant. This is the reason why we did not present it in Figure 10.2.

[FIGURE 10.3 OMITTED]

REFERENCES

[1] I. A. CHERNOV, J. BLOCH, AND I. E. GABIS, Mathematical modelling of U [H.sub.3] formation, Internal J. Hydrogen Energy, 33 (2008), pp. 5589-5595.

[2] J. CRANK, The Mathematics of Diffusion, Clarendon Press, Oxford, 1975.

[3] C. W. CRYER, On the approximate solution of free boundary problems using finite differences, J. ACM, 17 (1970), pp. 397-411.

[4] B. S. JOVANOVIC AND B. POPOVICH, Convergence of a finite difference scheme for the third boundary-value problem, Comput. Methods Appl. Math., 1 (2001), pp. 356-366.

[5] E. BLOCH, Finite difference method for the solution of free boundary problems, Phys. Fluids, 12 (1969), pp. 11-129.

[6] G. BECKETT, A moving mesh finite element method for the solution of two-dimensional Stefan problems, J. Comput. Phys., 168 (2001), pp. 500-518.

[7] O. A. LADYZHENSKAYA, Boundary Value Problems of Mathematical Physics, Applied Mathematical Sciences, Vol. 49, Springer, New York, 1985.

[8] R. D. RICHTMYER AND K. W. MORTON, Difference Methods for Initial-Value Problems, Second ed., John Wiley & Sons, New York, 1967.

[9] L. C. EVANS, Partial Differential Equations, Graduate Studies in Mathematics, Vol. 19, American Mathematical Society, Providence, RI, 1991.

[10] V. P. MIKHAILOV, Partial Differential Equations, Mir, Moscow, 1983.

[11] A. FRIEDMAN, Partial Differential Equations of Parabolic Type, Prentice-Hall, Englewood Clifts, NJ, 1964.

I. A. CHERNOV ([dagger])

* Received February 20, 2008. Accepted for publication December 12, 2008. Published online on April 14, 2009. Recommended by Y. Achdou. The work has been financially supported by the programme "Modern Numerical and Informational Technologies of Solving Large Problems" of the Department of Mathematical Sciences of the Russian Academy of Sciences.

([dagger]) 1185910, Inst Appl Math Research, Karelian Research Centre RAS, Pushkinskaya 11, Petrozavodsk, Russia (IAChernov@yandex.ru).
COPYRIGHT 2009 Institute of Computational Mathematics
No portion of this article can be reproduced without the express written permission from the copyright holder.