Computable estimates of the distance to the exact solution of the evolutionary reaction-diffusion equation

We derive guaranteed bounds of distance to the exact solution of the evolutionary reaction-diffusion problem with mixed Dirichlet-Neumann boundary condition. It is shown that two-sided error estimates are directly computable and equivalent to the error. Numerical experiments confirm that estimates provide accurate two-sided bounds of the overall error and generate efficient indicators of local error distribution.


Problem statement
Let Ω ∈ R d (d = 1, 2, or 3) be a bounded connected domain with Lipchitz continuous boundary ∂Ω, which consists of two measurable non-intersecting parts Γ D and Γ N associated with Dirichlet and Neumann boundary conditions, respectively. Let We consider the classical reaction-diffusion initial boundary value problem where n denotes the vector of unit outward normal to ∂Ω, and f (x, t) ∈ L2(Q T ), ϕ(x) ∈ L2(Ω), g(x, t) ∈ L2 (0, T ; L2(Γ N )) .
Assume that v ∈H 1 (Q T ) presents a certain numerical approximation of u. Our goal is to deduce accurate and explicitly computable estimates of the distance between u and v. For this purpose, we use the norm where ν, θ, and ζ are certain positive weights (weight-functions), which balance three components of the error. They can be selected in different ways so that (9) presents a common form of a wide spectrum of different measures. A fully computable and guaranteed upper bound of [u − v] 2 (ν, θ, ζ) is derived in Theorem 1 with the help of the method originally introduced in [15]. In [17], the method was applied to problems with convection, and in [12] the guaranteed error majorants were derived for the Stokes problem. In Section 2, we combine this method with the technique suggested in [14] for the stationary reaction-diffusion problem, which makes it possible to obtain the efficient error majorants for problems with drastically different values of the reaction function. Theorem 1 presents such an estimate for the problem (1)- (4). In Section 3, we derive a guaranteed and fully computable lower bound of the error (9) (Theorem 2). Sections 4 and 5 are devoted to practical applications of the estimates. In them, we discuss numerical results obtained for several typical examples, which confirm the efficiency of two-sided bounds.

Error majorant
Let e(x, t) := (u − v)(x, t) denote the deviation of v ∈H 1 (Q T ) from the exact solution u. From (8), it follows that Since e ∈H 1 (Q T ), we can set η = e and, using the relation This relation is a form of the 'energy-balance' identity in terms of deviation. It plays an important role in subsequent analysis. Now, we introduce an additional variable y ∈ Y div (Q T ), where Y div (Q T ) := y(x, t) ∈ L2(Ω, R d ) div y(x, t) ∈ L2(Ω), y · n ∈ L2(Γ N ) for a.e. t ∈ (0, T ) .
Theorem 1 (i) For any v ∈H 1 (Q T ) and y ∈ Y div (Q T ) the following inequality holds where C FΩ is the constant from the Friedrichs' inequality C tr is the constant in the trace inequality Here is a real-valued function taking values in [0, 1]; and α 1 (t), α 2 (t), α 3 (t) are positive scalar-valued functions satisfying the relation (ii) For any δ ∈ (0, 2], γ ∈ 1 2 , +∞ , and a real-valued function µ(x, t) taking values in [0, 1], the lower bound of the variation problem is zero, and it is attained if and only if v = u and y = A∇u.
Proof. (i) We transform the right-hand side of (12) by means of the relation which yields where By means of the Hölder's inequality, we find that and where ν 1 appears due to (6). Let µ(x, t) be a real-valued function taking values in [0, 1]. Then, we estimate the term I f as follows: In [14], this decomposition was used in order to overcome difficulties arising in the stationary problem if λ is small (or close to zero) in some parts of the domain (a more detailed study of this form of the majorant can be found in [11,9]). Combining (23), (23), and (25), we obtain The second term on the right-hand side of (26) can be estimated by the Young-Fenchel inequality where γ is an arbitrary real parameter from 1 2 , +∞ . Analogously, Here, α 1 (t), α 2 (t), and α 3 (t) are functions satisfying (18). The estimate (14) follows from (27)-(30).
On the other hand, if M 2 (δ, γ, µ) (v, y) = 0, then v satisfies the initial and boundary conditions, and for a.e. (x, t) ∈ Q T the following relations hold: From (31), it follows that The identity (32) is equivalent to (8). Hence, v = u, and we see that M 2 (δ, γ, µ) (v, y) vanishes if and only if This set of requirements is fulfilled if v coincides with the exact solution of the problem (1)-(4), i.e., e = u − v = 0 and y coincides with A∇u.
Remark 1 If g is a relatively simple function, e.g., piecewise affine function, then the function y may be selected such that g − y · n = 0 for a.e. (s, t) ∈ S N , and the constant C tr does not appear in the estimate.
Remark 2 An important question, which should be discussed in the context of a posteriori error estimation, concerns indication of local errors. We note that the majorant is presented by the sum of integrals, i.e., it automatically generates a sum of local quantities, which can be used as markers of local errors. In the numerical tests below, we show the efficiency of these error indicators.

Error minorant
Minorants of the deviation from the exact solution are well studied for elliptic problems, which have a variational form (see [13,10] and the literature cited therein). They provide useful information and allow us to judge on the quality of error majorants). Below, we derive computable error minorants for the evolutionary problem (1)-(4).
Theorem 2 For any v, η ∈H 1 (Q T ) the following estimate holds: Proof: It is not difficult to see that M(e) := sup where [e] 2 (ν, θ, ζ) is defined in (34). On the other hand, by using (8), we see that for any η the functional generates a lower bound of the norm [e] 2 (ν, θ, ζ) . Thus, we arrive at (34).

Incremental form of the estimates
and T N1×..×N d be a mesh selected on Ω. Therefore, Θ K×N1×...×N d = T K × T N1×...×N d denotes the mesh on Q T . Computational methods developed for parabolic type of problems often use incremental (semi-discrete) schemes. Below, we discuss forms of the majorants and minorants adapted to this class of methods. They lead to estimates, which evaluate errors on each interval [t k , t k+1 ] and accumulate the overall error. In this section, we assume that A = I, λ = λ(x), and S T = S D (these assumptions are made in order to simplify the notation).
We set where Assume that v is computed by a simple semi-discrete approximation method on interval [t k , t k+1 ] of length τ = t k+1 − t k (see, e.g., [19,1,4,8] Next, we define which yields Analogously, we deduce similar relation for the minorant. We set Then, on Q t k+1 the minorant is presented as where M 2 (k−1) is constructed on Q t k and Remark 3 We note that the presented incremental forms of the estimates (42)-(43) and (45)-(50) are valid only for the first order time discretization scheme. The corresponding higher order schemes can be derived by similar arguments, but the estimates will have a more complicated form. However, this subject is beyond the framework of this paper and will be considered in a subsequent publications. Also, we note that in the case of an oscillating right-hand side, we can combine the current estimate with the majorant or minorant of the modeling error (see, e.g., [16]).

Numerical tests
Example 1 We begin with a relatively simple problem where Ω = (0, 1), T = 10, ∂Ω = Γ D , u = 0 on S D , The quality of the error estimates is measured by efficiency indexes: . Table 1 shows these quantities for different values of the parameter δ (approximate solution was computed on the uniform mesh Θ K×N1 = Θ 40×40 ). We see that within the interval [0.5, 1.5] the efficiency of the majorant is not very sensitive with respect to the parameter δ. Results of other examples suggest similar conclusions. Therefore, in subsequent tests we set δ = 1 in order to obtain the optimal value of the majorant.
Growth of the error log [e] 2 (1, 0, 1) and the majorant log M 2 in logarithmic scale is depicted on Fig. 1. We see that the majorant reproduces the error quite accurately.  (1, 0, 1) Figure 1: Example 1. The error and the majorant with respect to time.  [u] 2 and the corresponding minorant M 2 [u] 2 . We see that the computed two-sided bounds of the error are efficient and guaranteed for any t k ∈ [0, T ].
If the weights of the error terms are selected in the special way, then the majorant and the minorant generate two-sided bounds for the same error norm (9). For any κ ∈ (0, 2 − δ), we have Then, the left-hand side is estimated from above by the majorant M 2 (see Theorem 14) and the right-hand side is estimated from below by the minorant M 2 (see Theorem 34) with the corresponding weights. Henceforth, we set δ = 1, then weighted error norms are denoted by [e] 2 (ν,θ,ζ) , whereν = √ 1 − κ,θ = √ κ CFΩ , andζ = 1. In minorant, , and κ 4 = 2. By changing κ, we obtain different weighted norms of the error, which have computable two-sided error bounds. In Table 3, we present the efficiency index of these two-sided bounds for different κ. In Fig. 3, we depict the behavior of    and κ = 10 −3 . The efficiency of the majorant (minorant) depends on the selection of y ∈ Y div (Q T ) η ∈H 1 (Q T ) . For elliptic problems, this question is well studied, and there exist methods which are able to reconstruct y using the computed flux A∇v (see, e.g., [18,20,5]). We follow the same technique and obtain a suitable y (η) by local minimization     Table 5 presents analogous results for the minorant, which were obtained by computing integrals in (34) on each time-step cylinder.
For the next set of numerical tests, we assume that the reaction term is positive and behaves as the Gaussian function, i.e., Then, the right-hand side of the problem is changed to Fig. 4, we see that for the certain σ λ the reaction λ changes rapidly from very small values (in one part of Ω) to relatively big values (in another part). The estimate (14) was derived specially for such type cases, and we use this example in order to verify it. Consider a simplified form of (14) with δ = 1 and β = const, which implies the majorant Minimization of the right-hand side of (54) with respect to µ is reduced to the auxiliary variational problem: findμ ∈ L ∞ (Ω) such that It is easy to find that for a.e. (x, t) ∈ Q T the minimizer satisfies the condition (57) Table 6 shows the efficiency of M 2 (μ) and M 2 (0) for different σ λ . In the left part of it, the results correspond to the case where y is reconstructed by piecewise affine approximations, and in the right part we expose the results obtained if y is taken from a reacher space (which includes piecewise quadratic functions). We see that M Example 2 Consider the same problem as in Example 1 but with ϕ = x sin(3 π x), A = I, λ(x, t) = ρ (t 2 + 1)(x + 10 −3 ), where ρ is a positive constant, and f = e t x (1 + 9π 2 ) sin(3 π x) − 6 π cos(3 π x) + λ(x, t) sin(3 π x) e t . The exact solution is u = x sin(3 π x) e t . Table 7 presents the efficiency of M The exact solution is u(x, t) = sin(π x) (t cos(t) + 1) (cos(π x) + 1). [u] 2 , the majorant M 2 [u] 2 , and the efficiency index for different meshes, where N 1 denotes the amount of intervals with respect to the space coordinate x and K with respect to the time t. We can see that the efficiency index stays on the approximately same level for all considered meshes, therefore the majorant does not deteriorate in the process of mesh refining. It is worth remarking that results exposed in Table  8 are quite typical, and similar behavior of the error majorant was observed in many other numerical tests.  Table 8: Example 3. The relative error, majorant, and its efficiency index with respect to different meshes Θ K × N 1 for t = T . In Fig. 6, we depict growth of the error and majorant in the logarithmic scale (for mesh Θ K×N1 = Θ 40×40 ). The gap between two curves reflects the impact of the term m 2 f (see Table 9) and corresponds to the efficiency index I M eff ≈ 1.87, . . . , 1.96. The efficiency of the majorant can be improved by choosing higher order approximations for the flux reconstruction (e.g., we can add element-wise based bubble functions; corresponding results are presented in columns 7 and 8 of Table 10). Table 10 illustrates minimization of M 2 with respect to y and corresponding efficiency indexes on every local time-cylinder Q k . For example, consider the row of Table 10 related to t k = 6.15.

The corresponding error [e]
2 (1, 0, 1) = 1.33e−03 and the estimate provided by the majorant with a simple (patch averaging) reconstruction of the flux is M 2 = 6.42e−03 I M eff = 2.20 . If we apply a more sophisticated procedure and reconstruct flux by minimizing the majorant with respect to values of y associated with the given spatial mesh on the time-layer t k = 6.15, then we obtain the efficiency index I M eff = 2.18. If we use a twice finer spacial mesh, then the index decreases up to 1.90. Obtained results show that for this particular example simple flux reconstructions generate sufficiently accurate estimates.
It is worth noting that, in general, the efficiency of the majorant depends on the several factors. First, it can be improved by using advanced reconstructions of the flux (e.g., by adding extra degrees of freedom, bubble functions, etc). However, this approach may lead to a limited effect, if approximations are coarse with respect to the time variable and/or the time steps are large (the term v t in the balancing term may be a piecewise constant function, which cannot reflect the behavior of f ).
where   The 'reliability term' m 2 f is necessary to provide a guaranteed upper bound, but the major part of the error is usually encompassed in m 2 d . The conclusion, which follows from our experience is that the term m 2 d (which is an integral over Ω) can be considered as an efficient indicator of element-wise error. Fig. 7 presents distribution of the error indicator m 2 d and error over time-step cylinder Q 30 (for zero reaction function). Assume now that λ(x) =   ϕ(x) = sin(π x) sin(3 π y) + sin(3 π x) sin(π y), A = I, and f = sin(π x) sin(3 π y) + sin(3 π x) sin(π y) cos(t) + 10 π 2 sin(t) + 10 π 2 t 3 + 3 t 2 + 10 π 2 .
In Table 11, we compare (1, 0, 1) [u] 2 with M 2 [u] 2 and its efficiency index (for the approximate solution computed on the mesh Θ K×N1×N2 = Θ 100×50×50 ). We can see that the majorant based on 'cheap' (local patch averaging) reconstruction of the flux y = G(∇v) provides a quite realistic upper bound of the error. However, in more complicated problems an optimization of the majorant with respect to y may be useful. This procedure yields sharper upper bounds but requires more computational efforts (concerning the corresponding methods based on multigrid, isogeometric elements, and other methods, see [20,5,9]).
Next goal is to investigate the accuracy of the error indicator defined in (59)-(60). We analyze two different measures, which can be called 'weak' and 'strong' and are discussed in details in [9]. The first measure is studied in the context of a certain marking procedure M, which maps element-wise error into a boolean array, i.e., it deals with the values 0 and 1 only. The corresponding 'weak' measure M weak ∈ [0, 1] is defined by the percentage of correctly marked elements.
To analyze the quality of the weak measure, we consider bulk marking procedure M θ , where θ ∈ (0, 1) (see, e.g., [2]). Fig. 9 illustrates M θ for θ = 0.2 and 0.4 (it has been performed for the actual error and for the error indicator m 2 d ). The results obtained for the error (Fig. 9 left) and for the indicator (right) are almost identical. We also have obtained quite small values of the weak error measure for bulk parameters θ = 0.
To understand whether or not the error indicator is quantitatively sharp and reproduces the error distribution accurately, we consider the histograms depicted in Fig. 10, which are constructed by the procedure suggested in [9]. We assume that all element-wise errors are ranked in the decreasing order with respect to values of the true error distribution [e] 2 (1, 0, 1) , and renumber all the elements accordingly, so that the element with the largest error is numbered 1. Then, we depict errors in this new order. The distribution of the element-wise errors generated by m 2 d is depicted in the same way. In Fig. 10, we consider true and indicated error distributions (the approximation computed on a regular mesh with 2 500 elements). If m 2 d is accurate in the strong sense, then the corresponding histogram (on the right) must resemble the histogram generated by the true error. Therefore, Fig. 10 shows that in this example the indicator is indeed sharp in the strong sense.  and f = sin(2 π x) sin(π y) cos(t) − sin(t) + 5 π 2 cos(t) + sin(t) + sin(π x) sin(2 π y) sin(t) + t cos(t) + 5 π 2 t sin(t) + 1 . (64) The exact solution is u = sin(π x) sin(2 π y) (t sin(t) + 1) + sin(2 π x) sin(π y) (cos(t) + sin(t)). We investigate the efficiency of bulk marking for θ = 0.5 and different time-layers for approximate solution constructed on the mesh Θ K×N1×N2 = Θ 200×50×50 (see Fig. 11). Here, M weak (M 0.5 ) = 5.08e−02 for Q 40 , M weak (M 0.5 ) = 5.4e−02 for Q 80 , and M weak (M 0.5 ) = 3.70e−02 for Q 160 .
The histograms from Fig. 12 are constructed by the same method as in the previous example (for 2 500 elements). They confirm quantitative efficiency of m 2 d . Again, we see that the majorant provides an efficient estimation of the global error as well as indication of element-wise errors. In this and also other experiments, we have observed that the quality of error estimation is good if the solution is smooth in space and rather monotonic in time, and it becomes less accurate if the solution admits large gradients with respect to the spacial variables and large time derivatives. In our experiments, this was mainly due to rather coarse (piecewise affine) approximations in time used for v and y. This fact constrains accuracy of the error estimation (especially in the context of the reliability term m 2 f ). However, choosing richer spaces for the reconstruction of v and y will lead to sharper estimates.  x (1) x (2) v