Method of lines for parabolic functional differential equations on cylindrical domains

The numerical method of lines is a technique for solving partial differential equations by discretising in all but one dimension. In this paper the solution of the approximate problem is extended outside the domain using the boundary condition. This leads to functional differential-algebraic equations. Sufficient conditions for the well-posedness, stability and convergence of the resulting method of lines are given.


Introduction
The method of lines (MOL) is a computational approach for solving PDE problems of the form 2. The resulting system of semi-discrete ordinary differential equations (ODE's) is integrated in time.
From the numerous literature concerning the numerical method of lines for classical differential equations we mention the monographs [8,25,27]. MOL for nonlinear parabolic functional differential equations with initial boundary conditions of the Dirichlet type are investigated in [30], where the error estimates implying the convergence of MOL on a rectangular domain is given. Shakeri and Dehghan [24] present two forms of MOL for one-dimensional nonlocal hyperbolic partial differential equation. MOL is used in [4] to obtain numerical solutions to a quasilinear parabolic inverse problem. Parabolic inverse problems can be reduced to a system of ODE's by fourth order compact scheme (see [18]). The aim of the paper is to construct a method of lines for nonlinear parabolic functional differential equations with general initial boundary conditions and with a non-rectangular domain. To deal with a cylindrical domain, we can proceed in several ways. For instance, we can consider the polar coordinate system, like in [25,27] or the Cartesian coordinate system. In the later case we have at least two possibilities. We can add extra points either at the boundary of the domain (see [14,17]), or outside the domain (see [19]). In this paper we consider the second possibility. In our scheme we obtain additional points using the reflection with respect to the boundary. We extend the solution of the approximate problem outside the domain by a transformation generated by the boundary condition. This leads to differential-algebraic equations (DAE). The theoretical analysis of DAE's and some appropriate numerical methods for initial and boundary value problems can be found in [3,15]. In [6] the authors study numerical solutions of DAE's. Soltaniana et al. [26] present a homotopy perturbation method to solve DAE's. The numerical solution of DAE's by the Rosenbroch onestage scheme with complex coefficients are investigated in [2]. However, there are considered systems which depend only on present time. Equations with aftereffects of various kinds (such as delays) are called functional differential equations (FDE) The theory of FDE's can be found in [7]. Semiexplicit numerical methods of the Rosenbroch type for functional differential-algebraic equations in the whole space were studied in [10]. In [13] the authors, on the basis of Newton's method, propose a fast quasilinearization numerical scheme, coupled with Rothe's method, for fully nonlinear parabolic equations. Comparison theorems for ODE and DAE systems are investigated in [11].
We admit some advantages of such scheme: The new approach is applicable for n-dimensional spatial variables, even though the numerical experiments are implemented for a two-dimensional spatial case. The weak point of the proposed method is a limitation on the boundary of the area of the spatial variable. Namely, we assume that the domain is bounded, convex with the boundary of class C 2+α , 0 < α < 1. We require this assumption to construct a suitably regular extension of the functions beyond the domain Q. Fakhar-Izadi and Dehghan [5] deal with irregular domains for a weakly parabolic partial integro-differential equation. They propose a spectral method to find numerical solutions. It seems that spectral methods are more time consuming compared with methods on regular meshes which obey a kind of maximum principle. Our technique may be used in many physical problems such as chemical reaction [9] and heat conduction [22], where a flow across the boundary surface is proportional to the difference between the surrounding density and the density inside surface. In these papers authors consider complex approximate methods, exclusively for the differential equations without any delay. In the article [23] one examines non-linearities and delays with Robin boundary conditions. In our paper we deal with difference schemes for functional differential equations with Robin boundary conditions, which is a much wider class that delays. We can also use our technique to a quasilinear parabolic system of PDE's, more precisely in the chemotaxis model considered in [29]. Recent advances in computational mathematical biology confirm our interest in theoretical investigations, compare: [16]. If the domain is not convex, as it often happens in realworld applications (chemotaxis, reaction-diffusion etc.), it is still possible to use our techniques by means of some smooth mappings which transform an irregular domain to a disc. For instance, a bean-shape region can be just projected to a disc.
The paper is organized in the following way. In Sect. 2 we set up the notation and terminology. Section 3 contains auxiliary lemmas. In Sect. 4 our main results are stated and proved. In the last section numerical experiments are presented.

Formulation of the problem
Let ⊂ R n be a bounded, convex domain with boundary ∂ of class C 2+α . Write where T > 0 and T 0 ∈ R + = [0, +∞). For a function u : Q 0 ∪ Q → R and for t ∈ [0, T ] we define a function u t : Q 0 → R by For any metric spaces X and Y we denote by C(X, Y ) the space of all continuous functions defined on X and taking values in Y . In the case Y = R we write C(X ). Given with the initial boundary conditions where n = n(x) is the unit outward normal on ∂ . Equation (1) with a particular right hand side can be interpreted as a reactiondiffusion equation, which is widely used as a model describing various physical, chemical and biological problems, see [1,16,20]. Note that the right hand side of the equation contains a functional variable. Therefore, we can also consider differential equation with deviated variables or differential integral equations, as it is shown in our examples 1 and 2.
We transform the boundary condition (2) by considering an extension of a function u outsideQ : Given u :Q → R and denoting r ( For any smooth function u this extension implies boundary condition (2) on [0, T ] × ∂ . If β ≡ 0, then this extension is a mirror reflection with respect to the boundary ∂ , generated by Neumann's conditions.

Discretization
We construct a regular mesh on R n in the following way. Let h = (h 1 , . . . , h n ), h i > 0 be the steps of the mesh. For m ∈ Z n , m = (m 1 , . . . , m n ), we denote nodal points in the following way: The difference operators δ (2) ..,n , are defined in the following way: Let us consider the interpolating operator In [12] (page 85) we find another extrapolation method. It is easy to see that for (t, x (m) ) ∈ [0, T ] × h , the algebraic equations for (t, x (m) ) ∈ [−T 0 , T ] × * h , and with the initial condition for (t, The method of lines (5)- (7) can be written as the abstract differential-algebraic problem with the initial condition K is generated by [a i j ] and ∂ q j f , while ϕ h by f . The choice is not unique. Here one can choose K dependent on a i j for the components with the indices not greater than p and ϕ h (t, u t ) dependent on f .

Assumption A[C]
The initial function ψ h ∈ C([−T 0 , 0], R P ) satisfies the consistency condition

Lemma 1 Suppose that A[K], A[σ ] is satisfied andφ is bounded and continuous.
Then problem I p× p 0 0 0 Proof Problem (11) can be written in the following form In this theorem we only consider the case when the right hand side is dependent only on t. Hence from the boundedness, DA-irreducibility and continuity of K there exist a unique solution of problem (12) (see [28] Thm VII). The component z A is uniquely determined by z D . The consistency condition A[C] guarantees the continuity of the solution.

Lemma 2 Suppose that A[σ ], A[K], A[C] are satisfied and |g i (t, ξ)| ≤ σ (t, ξ ) for i ≤ p and g i (t, ξ) = 0 for i > p. Then
where z is a solution of with the initial condition (9) and γ = ψ h γ = ˜ .
Proof Ifγ , γ ≡ 0, then according to [11] z ≡ 0. If we assume thatγ , γ ≥ 0 then we consider the comparison system and B i is the solution of the ODE Proceeding analogously like in [11] we have x i = x 1 B i B 1 , i = 2, . . . , P and It follows from the comparison principle for ODE's that which completes the proof.
The following lemma is crucial in the proof of Theorem 1.  (8,9).

Lemma 3 Suppose that A[σ ], A[K], A[C] are satisfied and
Proof Consider an iterative method, which starts from a prescribed function u 0 ∈ C([−T 0 , T ], R P ) satisfying (8) with some error˜ (t) andγ = ˜ . Now we consider the linear system of equations for each k = 0, 1, . . .
with the initial condition (9). It follows from Lemma 1 that there is exactly one solution of the above problem. Applying Lemma 2 to the differences u k+l − u k for k, l ∈ N we conclude that {u k } k∈N is the Cauchy sequence.

Stability and convergence
We are now in position to state our main stability and convergence result for the MOL corresponding to (1,2,3). We need the following assumptions on the functions f , β, a i j , and the steps h of the mesh.
Assumption A w, q, the same property have the derivatives ∂ q j f and they are bounded, a i j : Q → R are bounded and continuous in t for i, j ∈ {1, . . . , n}, -∂ q j f , a i j and steps h satisfy the relations (CFL) Proof Consider the iterative method. We choose an arbitrary function with the initial and boundary conditions (6)- (7).
We apply the Hadamard mean value theorem to the difference We have a i j t, x (m) δ (2) i j u (m) We substitute the formulas for δ, δ (2) in the above equation. The matrix K consists of elements which are linear combinations of a i j and ∂ q j f . According to Lemma 2 we show that the matrix K satisfies all conditions of A[K]. Put It follows from assumption A that is a convex combination of u(t, x (m) ) and first two conditions of A[K] are satisfied. Note that at least one coefficient in J h is positive. Thus last two conditions of A[K] are satisfied. It follows from Lemma 1 that there is exactly one solution of problem (13).
Next we analyse the difference u k+l − u k for k, l ∈ N. Put l u = u k+l − u k . We have We apply the Hadamard mean value theorem Applying Lemma 2 we have where lim k→∞ ω k (t) = 0 uniformly on [0,T ], for eachT ∈ (0, T ) (see [21]). We conclude that {u k } k∈N is the Cauchy sequence.

Theorem 2 Suppose that A[σ ], A are satisfied and
-u is a solution of (1,2,3) andũ is a solution of (5,6,7) such that Remark 1 We assume that sign a i j is constant. One can omit this assumption by considering J − , J + as sets which depend on x (m) .
Proof (Theorem 2) Let h : Q h → R be defined by the relation a i j t, x (m) δ (2) i j u (m) (t) It follows that there isγ h such that From the definition of J h we have Applying the Hadamard mean value theorem we have Hence the above equation can be written in the following form According to Lemma 2 we show that the matrix K satisfies all conditions of A[K]. Put It follows from assumption A that

Numerical examples
We apply the results presented in Sect. 3 to a differential equation with deviated variables and to a differential integral problem. We consider our numerical examples on the cylinder [0, T ] × B 1 , where B 1 is the unit ball centered at (0, 0). with initial and boundary conditions u(0, x, y) = 1 on B 1 , ∂u ∂n + 2t tan t (x 2 + y 2 + 1) u = 0 on 0,

Example 1 Consider the differential integral problem
The solution of the above problem is known, u(t, x, y) = cos(t (x 2 + y 2 + 1)).
In Table 1  The solution of the above problem is known, u = e −t (x 2 +y 2 ) .
In Table 2 we give experimental values of the maximal error for h 0 = 0.01 and h 1 = h 2 = 0.125, where (h 0 , h 1 , h 2 ) are steps of the mesh with respect to (t, x, y).
The above examples are carried out for two-dimensional spatial variables. This is done only for our convenience of implementation. The theory presented in our paper is not limited with respect to the dimension of spatial variables. Both coefficients of the derivatives of the unknown function and the functions appearing on the right hand side of the equation and the initial and the boundary conditions satisfy the assumptions imposed in our paper. The computed results have been compared with the exact solution to show the required accuracy of the method. The computation time is 0.38sec for the first example, and 24.82 s for the second example. The presented experiments illustrate the convergence of the proposed method.