Exact Solutions and Continuous Numerical Approximations of Coupled Systems of Diffusion Equations with Delay

In this work, we obtain exact solutions and continuous numerical approximations for mixed problems of coupled systems of diffusion equations with delay. Using the method of separation of variables, and based on an explicit expression for the solution of the separated vector initial-value delay problem, we obtain exact infinite series solutions that can be truncated to provide analytical–numerical solutions with prescribed accuracy in bounded domains. Although usually implicit in particular applications, the method of separation of variables is deeply correlated with symmetry ideas.


Introduction
In real-world problems, there are time lags between actions and responses. When these time delays are much shorter than the scale of observation, and they do not significantly affect the dynamics of the system, these problems might be satisfactorily modeled using ordinary or partial differential equations (PDEs). However, there are many situations where the presence of delays can not be safely ignored, requiring the use of modeling tools such as delay differential equations (DDEs) and partial delay differential equations (PDDEs) [1,2], including, among others, problems in life sciences [3], population dynamics [4], and control engineering [5].
For transport phenomena, diffusion, and heat conduction problems, it has been long pointed out that classical models, derived from Fourier or Fick laws and resulting in parabolic partial differential equations, implied infinite speed of propagation [6][7][8][9], and alternative models including the presence of delays have been proposed, finding increasing interest and wider areas of applications in recent years (see [10,11] and references therein).
Exact and analytical-numerical solutions for the generalized diffusion equation with delay, where a delay term is added to the classical model, have been previously obtained [12,13]. In the present work, we consider coupled systems of generalized diffusion equations with delay, written in matrix form as u t (t, x) = Au xx (t, x) + Bu xx (t − τ, x), t > 0, 0 ≤ x ≤ l, with solutions X n (x) = sin(λ n x), corresponding to the sequence of eigenvalues λ n = nπ l , n = 1, 2, . . . ; and the corresponding temporal initial-value vector delay differential problems T n (t) + λ 2 n (AT n (t) + BT n (t − τ)) = 0, t > 0, (6) T n (t) = B n (t), −τ ≤ t ≤ 0, where B n (t) are the Fourier coefficients in the expansion of the initial function ϕ(t, x) in terms of the eigenfunctions X n (x), i.e., so that Then, the formal series will be a candidate for the exact solution of problem (1)-(3). Problems (6) and (7) are of the general form with F(t), ψ(t) ∈ R M and A, B ∈ R M×M . We consider first the auxiliary matrix problem where G(t) ∈ R M×M and I is the identity matrix, whose solution is given in the next lemma.

Lemma 1.
Consider problem (13) and (14) with A invertible. Write C = A −1 B and let Then, the solution of (13) and (14) in the interval [(m − 1)τ, mτ] is given by Proof. It is clear that G(t) is a well-defined continuous function, as in each interval [(m − 1)τ, mτ] it is a sum of continuous functions, and values at the ends of the intervals agree. It is also immediate to check that the matrix functions Q k (t) defined in (15) satisfy Thus, for t ∈ [0, τ], one has Moreover, for t ∈ [(m − 1)τ, mτ], with m > 1, one has The next lemma gives an expression for the solution of problem (11) and (12) in terms of the function G(t) defined in Lemma 1.

Lemma 2.
Consider problem (11) and (12) with A and I + C invertible. For a differentiable initial function ψ(t), its solution is given by F(t) = ψ(t), for −τ ≤ t ≤ 0, and, for (m − 1)τ < t ≤ mτ and m ≥ 1, by Proof. It is clear that F(t) is continuous, as it is defined in terms of continuous functions. It is also immediate from (20) that F(t) satisfies (11) for t ∈ [(m − 1)τ, mτ] when m > 1. Thus, we only need to check that (11) is also satisfied for t ∈ [0, τ]. In this case, since G(t − τ − s) = I for s > t − τ, one has A more explicit expression for the solution is presented in the next theorem.
Theorem 1. The solution of problem (11) and (12), with conditions as in Lemma 2, is given by F(t) = ψ(t) for −τ ≤ t ≤ 0, and, for (m − 1)τ < t ≤ mτ and m ≥ 1, by where the second summation is assumed to be empty for m = 1.
Proof. It is immediate by substituting in (21) the expression of G(t) given in (17), taking into account that t − τ − s > (m − 1)τ ⇐⇒ s < t − mτ, and cancelling out some terms.

Remark 1.
Although the matrix functions Q k (t) have been defined in (15) recursively, they can be written explicitly as iterated integrals, When A and B commute, it is not difficult to check that they are given by the following compact expression without integrals, In particular, if A and B are diagonal, or with the appropriate change of variables when they are simultaneously diagonalizable, problems (11) and (12) consist of M independent scalar problems, and it can be checked that the expressions given by Theorem 1 for each component of F(t) agree with those given in [12,13] for the corresponding scalar problems.

Remark 2.
In scalar problems, diffusion coefficients are always positive, so it is common to assume in diffusion vector problems that the corresponding matrix coefficients are positive defined (see, e.g., [23]). In the next sections, where exact and numerical solutions for the coupled diffusion problems (1) and (3) are derived, we will assume the weaker condition of A being positive stable, i.e., having eigenvalues with positive real part, similar to the condition assumed in [16] for diffusion problems without delay. In this section, in Lemma 2 and Theorem 1, we have only required A to be invertible, and similarly for I + C, or equivalently A + B, which would be the coefficient when τ = 0. As will be indicated in the last section, the solution of (11) and (12) given in Theorem 1 might find application in different problems, not necessarily of diffusion type, so that only conditions guaranteeing the obtention of compact, closed form solutions of (11) and (12) have been assumed.
When A is singular, matrices Q k (t) can still be defined, by replacing the definition of Q 1 (t) in (15) with the nonintegrated form However, even in the much simpler case of commuting coefficients, infinite sums seem unavoidable, as the expression corresponding to (25)

Exact Infinite Series Solution
The solution of the general initial-value vector DDE problem given in Theorem 1 provides expressions for the functions T n (t) in (10), by taking in (11) and (12) The corresponding functions Q k (t) for each of these problems will be denoted Q n k (t). Writing p = π/l, the candidate series solution of problem (1) where In the next theorem, we show that under a condition on the eigenvalues of A and for sufficiently regular initial functions, the three series in (29)-(31) converge uniformly and can be differentiated termwise with respect to t and twice with respect to x, and the candidate series (29) is a classical solution of problem (1)-(3). Specifically, we will assume the following conditions on the initial function ϕ(t, x): Before proving this theorem, we present in the next lemmas bounds for the norms of matrix exponentials and of the functions Q n k (t). In what follows, denotes a vector norm or a compatible norm for matrices, and we will assume that the conditions of Theorem 2 hold. Lemma 3. For a matrix A, let σ(A) be the set of its eigenvalues, and α(A) = max{ (z)|z ∈ σ(A)}. Then, for each ν > α(A), we can find a constant K ν such that e At ≤ K ν e νt , ∀t ≥ 0. (33) Proof. Consider the Jordan decomposition A = PJP −1 , so that J = D + E, where D is diagonal, with eigenvalues of A, and E is an upper triangular nilpotent matrix. Letting P(t) = ∑ m−1 k=0 Et k /k!, one has ( [24], p. 396) that e At ≤ e σ(A)t P(t). Thus, for any ν > α(A), there is T > 0 such that e tα(A) P(t) ≤ e νt for t > T, and we can take K ν = max{max{P(t)|t ∈ [0, T]}, 1}.

Lemma 4.
Let α = A and γ = C . Fix ν > 0 such that −ν > α(−A), and write K for a constant K −ν > 1 satisfying the condition of Lemma 3. Then, Q n k (t) can be written in the form where the matrices S n k (t) admit the bounds Proof. We proceed by induction on k. For k = 1, one has Q n 1 (t) = (e −n 2 p 2 At − I)(I + C), and we can write Q n 1 (t) = −(I + C) + S n 1 (t), where S n 1 (t) = e −n 2 p 2 At (I + C), so that which is of the form (35). For k > 1, from (15) and assuming the induction hypothesis, one has Q n k (t) = − t 0 e −n 2 p 2 A(t−s) n 2 p 2 AC (−1) k−1 C k−2 (I + C) + S n k−1 (s) ds where Thus, We can now proceed to the proof of Theorem 2. We note that the conditions (3) assumed for the initial function ϕ(t, x) ensure convergence of the Fourier expansions of ϕ t (t, x) and ϕ xx (t, x).
Proof of Theorem 2. From Lemma 2, since ν > 0, it follows that for each k and any finite positive t, the functions Q n k (t) are bounded. In addition, for each t there is N such that Q n k (t) is decreasing for n > N.
Consider now the series obtained by termwise differentiating u(t, x), as given in (28) and (29)-(31), with respect to t. We denote this series u t (t, x), as it will be proved that it converges uniformly in (0, ∞) × (0, l). For each of the subseries in (29)-(31) and for each k, and letting aside some constant terms, one gets infinite series of the form for u 1 (t, x), and of the form for u 2 (t, x) and u 3 (t, x). From the conditions assumed in (3) for the initial function ϕ(t, x), it follows the uniform converge and continuity of the series and for x ∈ (0, l). Hence, since Q n k (t) are bounded and decreasing for n large, it follows that u t (t, x) converges uniformly in [(m − 1)τ, mτ] × (0, l). It is immediate to check the continuity at the connecting intervals for m > 1. For m = 1, at t = 0 it is easy to check the continuity of u(t, x), but in general this is not the case for u t (t, x), unless a special condition is required on the initial function ϕ(t, x).
Similar arguments apply for the series resulting from termwise twice differentiating u(t, x) with respect to x, since they have the same form as those previously discussed for u t (t, x).

Continuous Numerical Solutions
In this section, we will obtain bounds on the errors of continuous numerical solutions of problem (1)-(3), computed by truncating to N terms the exact series solution defined in Theorem 2.
Using the decomposition of Q n k (t) given in (34), it is immediate that the series in (29)-(31) can be written in the form where S i (t, x), i = 1, 2, 3, are the respective infinite sums in (29)-(31) corresponding to the terms S n k (t). Hence, the errors resulting from approximating u(t, x) by using the expressions for u i (t, x), i=1,2,3, given in (44)-(46) but computing only a finite number of terms, N, in the infinite sums S i (t, x) i = 1, 2, 3, are To bound these errors, we will use bounds for S n k (t) derived from (35), in terms of incomplete gamma functions ( [25], p. 174). Letting µ k = max{(α/ν) j , j = 0, . . . , k − 1}, from (35) it follows that We note that, from the conditions (3) assumed for the initial function ϕ(t, x), the Fourier coefficients of ϕ(t, x) and ϕ t (t, x) decay as O(n −2 ), so we can find constants H and H 1 such that, for s ∈ [−τ, 0], Hence, writing γ 1 = (I + C) −1 , one has since Γ(k, v) is decreasing with respect to v and ∑ ∞ n=N+1 1/n 2 ≤ 1/N. Similarly, for R N 2 , one has and with the change of variable v = n 2 p 2 ν(t − mτ − s), one gets since ∞ 0 Γ(m, v)dv = Γ(m + 1) and ∑ ∞ n=N+1 1/n 4 ≤ 1/(3N 3 ). Therefore, we have proved the following theorem. Theorem 3. Consider problem (1)-(3), let u(t, x) be the exact series solution given in Theorem 2, and let u N (t, x) be the approximation obtained when the infinite series' in (44)-(46) are replaced by the corresponding partial sums with N terms. Then, Consequently, for any δ > 0 and given a prescribed a-priori Next, we present an example showing the practical feasibility of computing the numerical solutions u N (t, x) and the error bounds given in Theorem 3. In this example, we used the infinity norm; computations were performed using Maple R , and graphics were prepared using Matlab R .
Example 1. Figure 1 shows the numerical solutions computed with N = 20 for problem (1)-(3) with parameters τ = 1, l = 1, and In the next figure (Figure 2), total error bounds and individual contributions to the error bounds by the three terms in (56) are presented.
In this example, the Fourier coefficients of the initial function satisfy more stringent conditions than those assumed in (51), since one has B n (s) = 12(−1) n (s+1) Thus, taking H = H 1 = 12/π 3 , one has B n (s) ≤ H/n 3 and B n (s) ≤ H 1 /n 3 . Hence, the bounds given in (56) can be refined for this example by substituting N in the denominator of the first two terms by 2N 2 , and 3N 3 in the denominator of the third term by 4N 4 .
As seen in Figure 2, the order of the global error bound is essentially determined by the third term in (56).

Conclusions
In this work, we have presented exact and continuous numerical solutions for coupled systems of diffusion equations with delay, extending previous results for scalar problems [12,13]. These solutions are expressed in terms of the parameters of the problem, and may be used to analyze how the system behavior depends on those parameters.
As shown in Theorem 3, the analytic-numerical solutions proposed in this work may provide numerical approximations with a-priori prescribed errors in bounded domains. The error bounds given in Theorem 3 for the first two terms, R N 1 and R N 2 , depending on gamma incomplete functions, are essentially exponentially decaying with N; and the third term, R N 3 , is O(N −3 ). As was shown in [13], the rate of convergence can also be made exponential for this last term, by using suitable polynomial approximations to the initial function. In any case, the error bounds guaranteed by Theorem 3 are usually far from being sharp, and convergence can be much faster, especially for problems with highly regular initial functions. Additionally, when the initial function consists of a finite combination of the eigenfunctions sin(npx), the exact solution defined in Theorem 2 reduces to a finite sum.
The expression given in Theorem 1 for the solution of a general initial value vector DDE problem may find applications in a wide variety of problems, beyond its use in this work. For instance, it would provide the basis to deal with related problems for coupled system with delay, as the type of reaction-diffusion equations with delay considered in [13]. More generally, since the matrix coefficients are not required to commute, this expression can be applied to obtain constructive solutions for problems involving scalar higher-order linear DDE, as they can be converted into a first-order system of the type considered in Theorem 1. The exact solution given in Theorem 1 could also provide the basis for extending previous works dealing with random scalar delay problems [26,27] to general vector delay problems, or with the construction of numerical schemes based on exact solutions [28][29][30].
Author Contributions: Conceptualization, E.R. and F.R.; methodology and formal analysis, all authors; writing-original draft preparation, F.R.; writing-review and editing, all authors. All authors have read and agreed to the published version of the manuscript.