Reduction of a boundary value problem for a system of diffusion-reaction equations to problem for a finite interval

The boundary value problem on a semi-infinite interval for a system of ordinary differential equations of the diffusion-reaction type is considered. To apply a difference scheme for the numerical solution of such a problem, it is necessary to reduce it to a problem for a finite interval. For this perpose we offer to use well-known method of construction of the stable variety of solutions that satisfy the limit condition at infinity. The estimates of the error arising when the problem is reduced to a finite interval are obtained.


Introduction
Various convective-diffusion processes are modeled on the basis of differential equations and of systems with a small parameter before higher derivatives. Solutions of such problems have large gradients in the boundary layer. This must be considered when building computational algorithms [1]. In mathematical modeling of a number of physical processes, for example, the transfer of impurities or chemical reactions, the limiting boundary conditions are set at infinity. To solve numerically such problem, it is necessary to reduce this problem to a finite domain. In this paper, we consider the problem of reducing the problem to a finite interval in the case of a linear system of ordinary differential equations with a limit boundary condition at infinity.
We consider a boundary value problem for a linear system of differential equations of diffusionreaction type on a semi-infinite interval. We use the method of reducing a problem to a finite interval on the basis of extracting a variety of solutions satisfying a given condition at infinity. The method of transferring conditions from a singular point was proposed by A.A. Abramov [2] and developed in a number of works by A.A. Abramov, N.B. Konyukhova and their co-authors in [3], [4] and in other works. It is shown how the developed approach can be applied in the case of a parabolic equation on the half-line. Now we introduce the notation. The Euclidean norm for vectors and the spectral norm for matrices are used. It is assumed that a matrix is nonnegative if its eigenvalues do not lie in the left half-plane. By C and C i , we mean positive constants that independent of the parameter ε and grid steps.

A boundary value problem for a system of equations on a semi-infinite interval
Consider the boundary problem for the system of equations reaction-diffusion type on a semiinfinite interval We assume that the vector function F(x) of order n and the matrix P (x) of order n × n are sufficiently smooth in the variable x. Let the matrix P (x) be positive definite, In part, the problem (1)-(2) was investigated in [5], where the following stability estimate is proved We investigate the method of reduction of the problem (1)-(2) to finite interval based on the extraction of the variety of solutions satisfying the limit condition at infinity [4].
So, we define the desired variety as the set of solutions of a system of equations of the first order: where the matrix G(x) is the solution of the singular Cauchy problem for the matrix Riccati equation: and the vector function β(x) is a solution of singular Cauchy problem The equations (5), (6) are defined so that the solution of equation (4) is a solution to the original equation (1). Let us dwell on the properties of the solution of the problem (5). In accordance with [5] G(x) ≥ √ αI if x ≥ 0. In accordance with [6], the following lemma is true.
Then from the conditions: Here, under the inequality U (x) ≥ 0 it is assumed that all the eigenvalues of this matrix are not negative. The proof of this lemma is based on the use of the maximum principle, which takes a place in accordance with lemma 1. Now we investigate the stability of the solution of the problem (5) to perturbation of the matrix P (x).
The matrix U (x), as the difference of two symmetric matrices, is symmetric. Lemma 2 implies the assertion of the lemma. Thus, in the case of a symmetric matrix P (x), the problem (5) has the only solution that is stable to perturbation of the matrix P (x), which leaves the matrix symmetric.
Then the conditions (7) are satisfied and therefore Ψ(x) ≥ 0, x ≥ L, which proves the lemma. Now we define the Cauchy problem with respect to the extracted variety (4): By analogy with lemma 4, we can prove the following lemma.

Lemma 5
The following estimate of stability is true Let us prove that the equation (4) extracts the variety of solutions of the equation (1) tending to zero at infinity.
Proof. Consider the problem Considering that β i (s) → 0, s → ∞, from the last relation it is easy to get that Let w = (Z, Z). We multiply each term of the last equation by Z(x) and get We account that W(x) is bounded and Transforming inequality (10) and integrating from 0 up x, we obtain It is easy to conclude that w(x) → 0, x → ∞. Therefore, Z(x) → 0 if x → ∞. This proves the lemma. Using (4), we transform the problem (1)-(2) to the problem for a finite interval [0, L] :  Proof. For each of the considered problems the estimates of stability (3), (9), (13) are true, respectively. The stability estimates correspond to the uniqueness of the solution of the problems under consideration.
Let W(x) be the solution of the problem (8). Now we prove that W(x) is a solution of the problem (1)- (2). In accordance with lemma 6, the limit boundary condition (2) is satisfied for the function W(x). The equation (1) is satisfied by W(x) according to the choice of the matrix G(x) and the vector function β(x). So, W(x) is a solution of the problem (1)- (2). Similarly, we can verify that W(x) is solution of the problem (11)-(12) with x ≤ L. It proves the lemma.
In accordance with lemma 8, the problem (1)-(2) can be exactly transformed to the problem on a finite interval: Now we investigate the stability of the problem (14)-(15) to perturbation of the coefficients G(L) and β(L), which can be calculated only approximately. Then From the fulfillment of these relations it follows that Ψ(x) ≥ 0 for all x. This proves (16) with C = (||Ũ(L)|| + 1)/ √ 2α. The theorem is proven. So, by (14)-(15) the problem (1)-(2) is precisely transformed to a finite interval. It remains to find the coefficients G(L) and β(L) from corresponding singular Cauchy problems. We'll find these coefficients by the expansion in the series [3], [4].
Substituting G m (x), β m (x) into the equations (5), (6), we get recurrent formulas It is easy to verify that the asymptotic approximation of G(x) in (17) is the solution of the problem (5) with the perturbed matrixP (x), similarly β m (x) is the solution of problem (6) with the perturbed vector functionF(x), with the error of order O(ε m+1 ).
Expansions in powers of x −1 . Applying such approach, we don't require the smallness of ε. May be ε = 1. Suppose that in (1) for x ≥ L the following decomposition is valid Functions G(x) and β(x) will be searched as Substituting these expansion into the matrix equation (5) and equating the coefficients with the same powers x, we get the recurrence formula: Let us construct the solution of the problem (6). Using the decomposition for G(x), we get the recurrence formula Then valid error estimates are To find G(x) based on asymptotic expansions it is required to solve the matrix equation AX + XB = F. Such methods were investigated in [7] and in other works.
So, the coefficients G(L) and β(L) for the problem (14)-(15) can be found with some given accuracy. In accordance with theorem 1, the solution of reduced problem is stable to arising errors.

Error estimation of classical artificial boundary condition
In a number of papers, the boundary condition is transferred from infinity to the point L in the form U(L) = 0 . Let us estimate the error of such a transfer.
Let us move from the problem (1)-(2) to the problem: Proof. Let us Z = U −Ũ. Then Z(x) is the solution to the problem Multiplying the equation and the boundary conditions scalarly by Z(x), we obtain relations for the scalar function w(x) = (Z(x), Z(x)) : Then It follows from the fulfillment of these conditions that Ψ(x) ≥ 0 for all x ∈ [0, L]. This proves the lemma. According to the lemma, the error decreases exponentially with decreasing x.

Parabolic equation on the half-axis
Consider the problem for a parabolic equation Functions b, f, φ 1 , φ 2 , are smooth enough and consistent so that the solution u(x, t) is a twice continuously differentiable function in all interior points of the area G. We suppose that Then the estimate max |u(x, t)| ≤ max |f (x, t)|/γ + max |φ 1 (t)| + max |φ 2 (x)| is true. We approximate the equation (19) with the step τ in time The problem (22)-(23) can be written in the form (1)-(2), while Further, on the basis of the developed approach the problem (22)-(23) can be reduced to a problem of the form (14)-(15) for a finite interval.
The reduced problem can be written as According to (24) by (26) the problem (19)-(20) reduced to a finite interval with the error of order O(τ ). The matrix P (x) from (25) is a diagonally dominant M-matrix, and therefore is positivedefinite. Consequently, the obtained estimates of the errors that arise when the problem is set on a finite interval are also valid in the case of the problem (26).
In the problem (26) it remains to replace the derivative (u i ) (x) using a finite-difference relation and on each i-th level in time to solve the system of linear equations.
From the analysis of the problem (5) with P (x) from (25) it follows that the matrix G(x) is lower triangular. In this case the matrix equation AX + XB = F, arising in the construction of an asymptotic solution, can be solved explicitly: where p = 0, 1, 2, . . . , i − 1 and i = 1, 2, . . . , n.
To use recurrent formulas for decomposition (17) we need to calculate the matrix root G 0 (x) = P (x). Such methods were investigated, for example, in [8], [9]. Taking into account that matrix P (x) is lower triangular, we can write an exact formula to find X = √ P . For an arbitrary element P i,i−p , with p = 0, 1, 2, . . . , i − 1 can be written i k=i−p X i,k X k,i−p = P i,i−p .
Then we get