On the a posteriori error analysis for linear Fokker-Planck models in convection-dominated diffusion problems

This work is aimed at the derivation of reliable and efficient a posteriori error estimates for convection-dominated diffusion problems motivated by a linear Fokker-Planck problem appearing in computational neuroscience. We obtain computable error bounds of the functional type for the static and time-dependent case and for different boundary conditions (mixed and pure Neumann boundary conditions). Finally, we present a set of various numerical examples including discussions on mesh adaptivity and space-time discretisation. The numerical results confirm the reliability and efficiency of the error estimates derived.


Introduction
Convection-diffusion systems are widely used in natural sciences and financial mathematics for describing physical flows and stochastically changing systems, respectively. The Fokker-Planck (or Kolmogorov) equation is one of the examples of this class illustrating the velocity of a particle in the Brownian motion, along with the Black-Scholes equation governing the price evolution, and the Navier-Stokes equation that describes the motion of viscous fluid substances. There are many more problems appearing in physics or biology that are concerned with modelling of semiconductors or of decision-making processes in neuroscience that can be targeted by this type of systems. The considered equation is sometimes also called the drift-diffusion equation.
This paper derives and discusses a posteriori error estimates of the functional type for convection-diffusion problems, which are motivated by the decision-making Fokker-Planck model problem discussed by Carrillo et al. [15] appearing in computational neuroscience. The mathematical model discussed in [15] is associated with stochastic dynamical systems, modelling the evolution of the decision-making (average firing rates) of two interacting neurone populations, see also [13] and [17]. In [15], the existence of a unique solution for stationary as well as evolutionary linear Fokker-Planck equations is discussed and proved under the assumption that the (regular enough) flux or drift is incoming in the bounded domain. The proofs are mainly based on the work by Wloka [56] as well as Evans [20]. In [56], a more general result is proved, which can be then applied in order to show existence and uniqueness in the time-dependent case too. Carrillo et al. [15] discusses this issue together with the large time behaviour of the problem and proves for that the convergence of the solution to the stationary state.
The main contribution of this work is the derivation of two-sided functional a posteriori error estimates (the so-called majorant and minorant) for stationary as well as time-dependent linear convection-diffusion problems. The ideas of functional type error bounds derivation were introduced by S. Repin in the 90's, see, e.g., [39,38,44,40,35]. Naturally, there are alternative approaches in a posteriori error control such as the residual-based or the goal-oriented estimation techniques, see, e.g., the books [54,3]. Moreover, we refer to the works [30,55] dedicated to convectiondiffusion problems. Compared to other methods, the functional type a posteriori estimation techniques provide a general (only functional based) framework to derive independent and guaranteed upper bounds. The initial work addressing elliptic convection-diffusion problems was done in [41], where the convection-dominated aspect was considered in detail (see, e.g., [41,Section 4.3.]). The derivation of majorants for the time-dependent reactionconvection-diffusion problem has been discussed in [46]. However, the numerical aspects as well as properties of the majorant applied to this class of problems has not been studied so far. With the current work, we fill this gap and consider various aspects of two-sided error bounds in applications. In particular, we present a set of In addition, we define the following weighted L 2 norm: v ε,Ω := ε Ω |v(x)| 2 dx 1/2 . (1) We define the H 1 spaces andṼ with underlying standard L 2 spaces. Here, the first space includes the homogeneous Dirichlet boundary conditions and the second one contains functions with zero mean in Ω. Moreover, we introduce the space Y div (Ω, Γ N ) := y ∈ L 2 Ω, R d div y ∈ L 2 (Ω), y · n ∈ L 2 (Γ N ) of vector-valued functions satisfying the identity Ω div (y w) dx = Ω (div y w + y · ∇w) dx = Γ N (y · n) w ds, ∀w ∈ V 0 .
Finally, the classic trace inequality is given by

Fokker-Planck model problems
The state Fokker-Planck boundary value problem (BVP) reads as follows where f ∈ L 2 (Ω). In the following, we consider the problem (12)-(13) for two cases, which differ with respect to (w.r.t.) the boundary conditions imposed on it. In the first problem, we impose mixed BCs of the form on (12)- (13), and in the second one, Robin BCs as follows Here, n denotes the vector of unit outward normal to ∂Ω, and the function F satisfies the following regularity requirements and the incoming flux conditions F · n < 0 on Γ N , and F · n < 0 on Γ, (18) in the case of mixed and Robin BCs, respectively. The time-dependent version of the problem (12)- (13) with Robin boundary conditions is discussed in section 4 as well. It reads as follows: find u ∈ H 1,1 (Q T ) satisfying where f ∈ L 2 (Q T ), and u 0 ∈Ṽ .
Function F satisfies (17) for a.a t ∈ (0, T ) and the incoming flux condition F · n < 0 on S T .

Two-sided estimates for the static case
The current section is dedicated to the derivation of the majorant and minorant of the distance between any approximation function v and the exact solution u of static or time-dependent Fokker-Planck problem. Presented upper bounds can be viewed as a generalisation of the majorant studied in [41,Section 4.3.] and [46]. The derived lower bound of the error is a completely original result.

Error majorant in case of mixed boundary conditions
By testing (12) with a function η ∈ V 0 , we arrive at the generalised formulation of (12)-(13) with boundary conditions (14)-(15): find u ∈ V 0 satisfying the integral identity where the non-symmetric bilinear form a(·, ·) and the right-hand side (RHS) are given by respectively. Assume that function v ∈ V 0 is the approximation of solution u and the distance e = u − v ∈ V 0 between them is measured in terms of the norm where ν := ε is a positive weight parameter, and δ 2 := 1 2 divF and χ 2 := − 1 2 F · n 1 /2 are positive weight functions.
By transforming (24) and assuming η = e, we arrive at the following equality Substituting the identity into (27), we obtain the relation known as 'error-identity' for e ∈ V 0 Following the ideas exposed in monographs [35,41,33], we introduce an auxiliary vector-valued function y ∈ Y div (Ω, Γ N ) satisfying the identity (5). By means of substituting (5) into (29), we arrive at From now on, let r eq , r d , and r ΓN denote the residuals in the RHS of (30) r eq := f − div(Fv) + divy, r d := y − ε∇v, and r ΓN : which also denote the residuals of the equations (12), (13), and (14), respectively. In order to estimate the RHS of (30), we apply the Hölder inequality, i.e., For treatment of the terms r eq and r ΓN , we introduce weight functions µ, θ from the space L ∞ [0,1] (Ω) in order to obtain robust majorants w.r.t. drastic changes in terms of divF and F · n. Then, the terms with residuals r eq and r ΓN can be estimated as follows: and The resulting majorant follows from the collection of (32)-(34) obtained by means of Young's inequality: where β and ζ are positive parameters. Therefore,

Error majorant in case of Robin boundary conditions
Next, we consider the state Fokker-Planck BVP (12)- (13) with Robin BCs (16). By testing (12) with a function η ∈Ṽ , we arrive at the generalized formulation of (12)-(13) with (16): find u ∈Ṽ satisfying the integral identity a(u, η) = f, η Ω for all η ∈Ṽ , where a(u, η) and f, η Ω are defined in (25). Assume that function v ∈Ṽ is the approximation of solution u and the error e = u − v ∈Ṽ is measured in terms of the norm (26). By transforming a(u, η) = f, η Ω and assuming η = e, we arrive at the error-identity which is similar to (29). Analogously to the previous subsection, we introduce an auxiliary vector-valued function y ∈ Y div (Ω, Γ), which yields the relation In order to estimate the RHS of (38), we apply the Hölder, Young, Poincaré, and trace inequalities leading to where

Error minorant in case of mixed boundary conditions
In the following, we derive a lower bound of the approximation error in terms of a different quantity.

Theorem 1
For v ∈ V 0 , the following inequality holds: Proof: The supremum can be estimated as follows Substituting identity (28) then yields the following estimate: which finally leads to equation (41).
Remark 1 Note that the new norm differs from the norm ||| e ||| 2 (ν,δ,χ) only by the weighting parameters that are given byν = ν 2 andδ 2 = δ 2 + |F| 2 2ε . The relation can also be written as follows ||| e ||| 2 (ν,δ,χ) = ||| e ||| 2 (ν,δ,χ) + When maximising the minorant functional, the auxiliary function w ∈ V 0 must be chosen from the richer approximation space in comparison to the approximation v ∈ V 0 . Otherwise, the minorant vanishes. Analogously to the majorant functional, the minorant includes only known data and, therefore, is fully computable. By an appropriate selection of w, we can find a lower bound arbitrarily close to the exact error. In particular, if w = u − v, the minorant coincides with the error norm. We consider now the parabolic initial-boundary value problem. After multiplying (19) by a test function η ∈Ṽ 0 , we arrive at the generalized formulation of (19)- (22): find u ∈Ṽ 0 satisfying the integral identity We present the functional error estimate, which provides a guaranteed upper bound of the deviation e = u − v between the generalised solution of (46) and function v ∈Ṽ 0 (generated by any numerical method) measured in terms of the norm where ν and ζ are positive weight parameters, and θ, χ are weight functions. The initial step in the derivation of both upper estimates is the transformation of (46) into the integral identity By substituting η = e and using the relation we obtain Taking into account the equality we rewrite the LHS of (50) as follows 1 2 e(·, T ) 2 Ω − e(·, 0) 2 We follow the ideas from [45] and introduce an arbitrary vector-valued function y ∈ Y div (Q T ) leading to the equation The residuals in the RHS of (53) correspond in a natural manner to the equations (19), (20), and (22), respectively. We proceed by estimating the terms in (53) using the Hölder, Young, Friedrichs as well as trace inequality, i.e., where functions α i > 0, i = 1, 2, 3. Finally, we arrive at the estimate where parameter δ ∈ (0, 2].

Lower estimates for the time-dependent case
Minorants of the error between the approximate and exact solutions for the evolutionary reaction-convectiondiffusion problem provide useful information while testing upper error-bounds (when u is not available). The quality of the majorant is evaluated by considering its ratio w.r.t. the minorant. The first results related to lower error bounds for evolutionary problems were derived and numerically tested in [34]. The result presented in the current section can be viewed as generalisation of Section 3 in [34].
Theorem 2 For any v, η ∈Ṽ 0 the following estimate holds: where Proof: It is not difficult to see that we find that which is exactly defined as [e] 2 M,Q T in (59). On the other hand, by using (46), we see that for any η the functional generates a lower bound for the norm [e] 2 M,Q T . Thus, we arrive at (59).

Remark 3
Since the majorants and minorants are equal to norms, see (35), (41), (47) together with (58), and (59), they are nonnegative. In [15], positivity for the solution in the static as well as time-dependent case is proved under the condition that F ∈ C 1 (Ω, R 2 ) and F · n < 0 on Γ N or F · n < 0 on Γ in the case of mixed or Robin BCs, respectively. For the time-dependent case, the nonnegativity is provided for given nonnegative initial values u 0 . Existence of a minimizing sequence for the majorant as well as a maximizing sequence for the minorant, see for instance [33], and corresponding infimum of the majorant and the supremum of the minorant, leads to the estimates of the accuracy of the approximation w.r.t. exact solution. However, it does not automatically imply positivity (or non-negativity). Together with the results on the existence of a minimizing sequence for the majorant as well as a maximizing sequence for the minorant, see for instance [41,Section 3.6], the infimum of the majorant and the supremum of the minorant both lead to the exact (nonnegative) solution.
Remark 4 Note that the model problems (12)- (13) and (19)- (22) are written in form of a first-order system of conservation laws supplied by a proper constitutive relation for the flux. This could yield to the extension of this work's functional-based error analysis to the use of a dual mixed weak formulation of the PDE system, which would be of great interest in the implementation of computations for more realistic applications, where the approximation of the total flux (Fu− p) is often more important than the approximation of u itself. Studies concerning the functionalbased error analysis devoted to a posteriori error estimates of mixed approximations for the Poisson equation with Dirichlet and mixed Dirichlet/Neumann boundary conditions can be found in [42] and [43], respectively. Moreover, in the latter work it is shown that after an appropriate scaling of the coordinates and the equation, the ratio of the upper and lower bounds for the error in the product norm never exceeds the constant 3. Such an a posteriori error approach is computationally very cheap and can also be used for the indication of the local error distribution. As an application, [43] considers the linear elasticity problem, where the reconstruction of the stresses (dual variable) is equally important to the reconstruction of the displacement (primal variable).

Numerical examples
We have presented a generalized error control method, namely, the functional error estimates, for convectiondominated diffusion problems in Sections 3 and 4. They are derived independently from the numerical method used to reconstruct the approximation. As a result, one does not need to adapt the obtained error estimates to the stabilization technique chosen in the solver but can just use the reconstructed approximation. The current section presents numerical results for the problems (12)- (13) and (19)- (22) in the static and timedependent case, respectively. The examples were computed for various parameters and boundary conditions. A particularly interesting setting (discussed throughout this section) is the convection-dominated problem, namely, when ε F L ∞ . In that case, the solution contains a boundary layer, which can be either of the width O(ε) (the so-called regular boundary layer) or of the width O( √ ε) (parabolic boundary layer). The detailed discussion addressing properties for both of these layers can be found in [50].
The stabilised form of (24) can be written as follows: where s K (u, η) depends on the approach used and δ K is a non-negative stabilisation parameter, which can vary from one method to another. In this work, we use the Streamline Upwind Petrov-Galerkin (SUPG) introduced by Brooks and Hughes [12], where The stabilisation parameter δ K is defined depending on the cell's Peclét number, i.e., where · ∞,K is the norm in [L ∞ (K)] d , h K denotes the diameter of a finite element K, and F prescribes the direction of the convective flow. Then, where δ 0 , δ 1 > 0 are appropriately chosen constants. Set δ K := h K 2 |F| , where |F| stands for the magnitude of the vector. All the numerical results obtained in the section are carried out in Python and C++ in the framework of The FEniCS Project [52,32].  Table 1: Example 1. Error order of convergence (e.o.c.) for approximations v ∈ P 1 and y ∈ RT 1 as well as v ∈ P 2 and y ∈ RT 1 w.r.t. refinement steps. Example 1 Consider the one-dimensional classical example on the unit interval Ω = (0, 1): with the solution u = 1−e a x/ε −e a /ε known to have boundary layer in the neighborhood of the right part of the boundary x = 1. We select a = 1.0 and ε = 1e−2. Its approximation v is reconstructed for two cases with Lagrangian finite elements (FEs) of order 1 (v ∈ P 1 ) and of order 2 (v ∈ P 2 ). The correctness of the numerical solver is confirmed by the error order of convergence (e.o.c.) for the uniform refinement procedure executed for both v ∈ P 1 and v ∈ P 2 (see Table 1). The same results are illustrated in Figure 1. In both cases, the initial mesh contains of 8 elements and 9 nodes, and we proceed for 12 global refinement steps. As expected, the error behaves as O(h) and O(h 2 ) in the first and the second case, respectively. Here, the flux y is reconstructed with Raviart-Thomas elements of the first order denoted by RT 1 as follows The auxiliary function w for the minorant reconstruction is approximated by P 3 -Lagrangian FEs. The values of the majorants and minorants as well as their efficiency indices w.r.t. several refinement steps are presented in Table  1. It is easy to see that both upper and lower bounds stay rather sharp w.r.t. the error measures [e] M and [e] M , respectively.
Example 2 Next, we consider the two-dimensional example defined on the unit square domain Ω = (0, 1) 2 : where ε = 1e−3 and F = (−x, 0) T . In this case, the exact solution is not known. However, by following the analysis presented in [47], one can predict the location of the boundary layer. The approximation v is reconstructed by     P 1 -elements, whereas the majorant is reconstructed with an auxiliary function y ∈ RT 1 and the minorant with an auxiliary function w ∈ P 3 . Let us first consider the numerical experiments obtained by globally refining the mesh. The order of convergence of M and M is confirmed in Figure 2. Due to their high quantitative efficiency, the lines of the error estimates and corresponding errors almost coincide on the plots. Table 3 illustrates values of both minorant and majorant as well as their efficiency indices w.r.t. the error measures that they bound.
Next, we test the performance of an adaptive strategy of the refinement. It is based on the local distribution of the indicator, following from the majorant functional. We choose a widely-used bulk marking criterion for selecting elements to be refined denoted by M BULK (θ), see [18]. Here, the bulk parameter is set to θ = 0.3. Table 5 provides the values for the efficiency indices of the majorant and minorant, which are getting sharper as the refinement steps proceed. The evolution of the meshes is illustrated in Figure 3. Here, the left column presents the meshes constructed on refinement steps 3-4, whereas the right column illustrates the distribution of the local error values on the elements of the mesh. The elements with the highest and lowest errors are indicated by the range of colours between red and blue, respectively. From the meshes, one can also see that refinement based on the majorant resolves the boundary layers at x = 0.
Example 3 Next, we consider yet another convection-dominated example on the unit square Ω = (0, 1) 2 :     where A = 10 0 0 1 , f = x + y, and F = − cos( π 3 ) x, sin( π 3 ) y T . Again, let the approximation v be reconstructed with P 1 -elements, y by RT 1 -elements, and w ∈ P 3 (the optimal error convergence is presented in Figure 4b). The approximate solution obtained with an adaptive refinement is illustrated in Figure 4a, where one can see that the solution convects towards the boundary y = 1. By the current example, we aim to show that even for the anisotropic A and heterogeneous F the performance of the error estimates remains rather sharp. In particular, Table 5 confirms that the efficiency indices of the minorant and the majorant lie close to 1 as the refinements evolve.   Table 6: Example 4. Error order of convergence for approximations v ∈ P 1 and y ∈ RT 1 for bulk parameters θ = 0.3 and θ = 0.4.

Static convection-dominated problems
Next, we discuss the following set of example problems: consider the reaction-convection-diffusion equation (which basically yields from (12)-(13)) with homogeneous Dirichlet BCs. However, in contrast to the previous examples, we assume now that reaction and convection are independent of each other.
Example 4 Assume that d = 2, a = 1, 0 T , ε = 1e−3, λ = 1, and the exact solution is defined by the function y) with a boundary layer in the neighbourhood of x = 1. Its approximation v is reconstructed by P 1 elements (see Figure 5a). Since the uniform refinement does not provide expected orders of convergence, we apply an adaptive refinement with bulk marking criterion M BULK (θ). We choose two different parameters, i.e., θ = 0.3 and θ = 0.4, to confirm optimal convergence order for M (see Table 6 and Figure 5b). The values of the total majorants and minorants as well as corresponding efficiency indices are presented in Table 7 w.r.t. several refinement steps. It is easy to see that even though the efficiency indices are high on the initial steps, both majorant and minorant become rather sharp towards the 9-th refinement step. The evolution of the meshes corresponding to the different θ parameters is presented in Figure 6.     According to the theory of asymptotic expansion [23], the solution of the original problem depends continuously on ε. Using this fact, the authors of [16,51] use the so-called multilevel-homotopic-adaptive finite element method (MHA FEM) with respect to the diffusion constant ε, such that grids are iteratively adapted to a better approximate solution as the diffusion parameter gets decreased. The method gets initialised withε = ε 0 (e.g., ε 0 = 1), and the mesh of the ad hoc diffusion problem −ε ∆ u + a · ∇ u + λ u = 0 gets refined adaptively. On the next iteration, the value ofε is reduced, and the mesh adaptation is performed for the problem with the updatedε. Iterations continue until the desired value of ε is reached, i.e., until conditionε > ε is satisfied. For the reader's convenience, the steps described above are summarised in Algorithm 1 (see also [51]). We use the MHA FEM in combination with the streamline diffusion FEM, where mesh adaptation is driven by the local distribution of the error indicator following from the majorant. Such an approach can also serve as a preprocessing for a problem so that one can generate problem-adapted meshes once realistic problems are concerned.
Example 5 Assume now that d = 2, a = 1, 0 T , ε = 1e−4, λ = 0, and the RHS f as well as and the boundary data g = 0 are chosen such that the solution of the problem is defined by Figure 7 illustrating the approximate solution on the uniformly refined mesh). It exhibits an exponentially regular boundary (of the width O(ε)) near {x = 1} perpendicular to the convection direction 1, 0 T (see also [51]). We test two different scenarios of updating the parameter ε. The first one contains 4 adaptive refinement steps for each of the levels ε, and the diffusion parameter decreases by 10 on each level. Table 8 contains the numerical results obtained by such an approach. Another way of 'switching' from level to level is refining until h min becomes small enough (mesh becomes fine enough) to represent the current boundary layer O(ε). The numerical testing of this approach is summarised in Table 9. We see that with the first approach, we are able to reach the error 16e−02, on the homotopic level ε = 1e−3. Therefore, choosing the number of adaptive refinements based on the approximation or mesh properties seems to be the more advantageous approach. Such a strategy can be used as an efficient mesh-generation technique to provide an initial setting for strongly convection-dominated problems (see Figure 8). Here, each presented mesh corresponds to the final mesh on each of the ε levels.
We also present convergence rates for the L 2 -norm of the errors that has been addressed in [51]. The obtained Refine mesh K h using refinement procedure R, i.e.,  Table 10).  Table 11.     Table  9.  Mesh obtained on the different homotopy levels using error indicators following from M µopt for v ∈ P 1 , y ∈ P 2 with bulk marking θ = 0.1.
However, in our experiment this accuracy is obtained for 2.8 times less d.o.f. then in [51]. We also present the mesh obtained by the adaptive refinement (driven by the indicator, which yields from the majorant), see Figure 10b. From the comparison of Figure 10b to Figure 9 in [51,Example 5], it is easy to see that in our case the refinement is rather concentrated on the region with a boundary layer rather than the region, where solution is smooth.

Time-dependent problems
In the next subsection, we consider the time-dependent reaction-convection-diffusion problem that follows from (19)- (22). We assume, that the balance equation (19) is written in the form such that the reaction and convection are assumed to be independent of each other. Since the incremental timestepping method possesses several disadvantages, for instance, time-consuming sequentiality which leads to complications in the parallelisation, we use a variational space-time approach to solve (19)- (22). For our numerical examples, the generalised formulation (24) is discussed in the Petrov-Galerkin setting with different test and ansatz spaces (following the spirit of [49]): find such that ε (∇u, ∇η) Q T + (λ u, η) Q T + (a · ∇ u, η) Q T + (u t , η) Q T =: a(u, η) = (η) : for all η ∈ W := L 2 (0, T ; H 1 0 (Ω)). The stability condition is provided by the extension of [   holds with a positive constant c. After 'homogenisation' of the problem by splitting u =ū(x, t) +ū 0 , (x, t) ∈ Q T , whereū 0 is an extension of u 0 ∈ H 1 0 (Ω) to the rest of the space-time cylinder. Hence, we aim to findū ∈ V 0 such that Using the fact that V 0 ⊂ W , we can state existence and uniqueness of the solution of the weak formulation (70) based on [48,Theorem 3.7]. Let the discretisation spaces be V h0 ⊂ V 0 and W h ⊂ W such that the inclusion V h0 ⊂ W h is satisfied. Then, the discrete version of (70) can be presented as follows: findū h ∈ V h0 satisfying The unique solvability of (71) follows from the discrete stability condition which is proved in [49,Theorem 3.1]. Finally, we introduce the FE spaces for u h ∈ V h0 and w h ∈ W h , i.e., where P 1 (Q h ) is the FE space with a basis provided by piecewise linear and continuous functions (e.g., Lagrange polynomials). This brings us to the final approximation result forū ∈ V 0 andū h ∈ V h0 . According to [49,Theorem 3.3], by assuming that u ∈ H 2 (Q), we obtain the following energy error estimate with ac > 0 independent on h. More general estimate follows as a corollary [49,Corollary 3.4].
follows. Here,C is a generic positive constant, independent from the mesh parameter h. This theoretical result is confirmed by results of numerical tests (see Table 13 and 14).
Since the majorant is defined as an integral over the whole cylinder Q T , it comes rather natural to use it for the error control. It also does not dependent on the discretisation used, therefore can be applied to the solutions reconstructed with fully unstructured grids. Next, we demonstrate numerical results obtained by applying functional type error estimates in case of a space-time discretisation.
The corresponding solution is illustrated in Figure 11. In order to provide robustness of the majorant w.r.t. changing values of λ, analogously to the elliptic case, we consider the additional function µ ∈ L ∞ [0,1] (Q T ), which is used to balance the contribution of λ into the r eq term, i.e., After setting parameters α 1 = 1 + β, α 2 = 1 + 1 β , with β = const, and α 3 = 2 in (58), assuming that the initial conditions are satisfied exactly, i.e., e(0, x) = 0, the final form of the estimate can be written as follows     Figure 13: Example 7. The mesh obtained on the sixth refinement step using M µopt (first row) and M µ=1 (second row). Here, v ∈ P 1 , y ∈ P 2 with bulk marking θ = 0.6.
The minimisation of (77) w.r.t. µ is equivalent to the following variational problem: findμ ∈ L ∞ [0,1] (Q T ) such that Functional Υ(µ) admits its minimum in µ opt (x, t) = C 2 FΩ (1+β)λ βγ+C 2 FΩ (1+β)λ . Table 13 provides the efficiency indices of M µopt and M µ=1 for different cases (a)-(b). The results correspond to the approximation v ∈ P 1 and the auxiliary flux y ∈ P 2 . The majorant M µ=1 grows dramatically in the case (a) and diverges according to the efficiency indices for the cases (b) and (c). At the same time, M µopt performs efficiently for all considered λ. It stays efficient and robust even if the reaction function changes its values drastically in different parts of the domain. The advantage of using an error estimate with optimal weighting function µ is illustrated also by the meshes in Figure 13. Using M µ=1 in all cases (a)-(c) results in a refinement procedure generating meshes that do not make changes in the solution visible (see the lower row of Figure 13). The meshes presented in the upper row, however, reflect the approximate solution accurately.
If we assume that λ = 0 and a = 0 (F = 0 in the problem statement (19)- (22)), the test case reduces to the diffusion equation of parabolic type with right-hand side f = u 0 (cos 2t − 2t sin 2t) + π 2 cos πx (t cos 2 t + 1) . Assuming that f ∈ L 2 (Q) and u 0 ∈ H 1 0 (Σ 0 ) such a problem is uniquely solvable in V ∆x 0 , where and u depends continuously on t in the H 1 0 norm. If u, v ∈ V ∆x 0 are provided, the time-dependent diffusion equation (19) imposes the error identity (see [1]): The results corresponding to the performance of the majorant M µopt as well as E Id in comparison to the errors they bound are illustrated in Table 14. We see that for the setting (I) we get very close to the value one for the efficiency indices of M µopt , and E Id is always identical to the values of |||e||| 2 L . However, an approximation of v by linear affine elements does not provide the condition v ∈ V ∆x 0 required for application of E Id. Therefore, stagnation of |||e||| 2 L , as well as E Id, is expected. For the setting (II), we see that the strong norm of the error and the error identity begin to converge as O(h). We also illustrate the results for the setting (III) in order to show that the performance of M µopt improves for smoother approximations of y (see the 5-th column). As we have observed in the numerical experiments, considering even smoother approximations of the flux does not pay off in this case, i.e., the minimisation of the majorant becomes even more time consuming, whereas the efficiency indices improve only slightly. The error identity stays quantitatively identical to the error norm |||e||| 2 L and, since it is only dependent on the approximation v, it does not require any computational time overhead. At the same time, the restriction on the regularity of u, v ∈ V ∆x 0 still remains, which makes the application of the error identity rather restricted to a certain class of problems.
One of the ways to study the efficiency of the majorant when it comes to generating the mesh refinement procedure is to compare the numerical results for the refinement steps of adaptive meshes, where the refinement is based on the element-wise distribution of the majorant, with those meshes refined w.r.t. the the element-wise errors. Figure 14 provides such a comparison of meshes obtained on the refinement steps 4, 6, and 8. We see that the meshes obtained in the upper row are topologically close to the meshes in the second row.

Conclusions
We have studied functional type a posteriori error estimation for stationary and time-dependent linear convectiondiffusion problems. This work was motivated by a decision-making Fokker-Planck model problem appearing in computational neuroscience, which has been discussed in [15]. We derived guaranteed two-sided estimates for the elliptic as well as parabolic problem and extended the a posteriori error analysis earlier derived in [41,46] by studying the numerical aspects as well as properties of the two-sided error bounds applied to this class of problems.
We presented a set of various numerical examples for different parameters and boundary conditions showing the sharpness of the upper and lower bounds in practice. In particular, we addressed convection-dominated problems including adaptive mesh refinement based on a bulk marking criterion from [18]. Furthermore, we presented numerical experiments applying the multilevel-homotopic-adaptive finite element method in the stationary case and space-time discretisation for the time-dependent reaction-convection-diffusion problem. The numerical results in both cases, static as well as time-dependent, successfully demonstrated the capability of the two-sided error bounds to provide a basis for the implementation of reliable and efficient computational methods.   (second row). Here, v ∈ P 1 , y ∈ P 2 with bulk marking θ = 0.6.