1 Introduction↩︎

Mean field games (MFG) model the Nash equilibria of differential games for a large population of players. MFG were introduced by Lasry and Lions [1]–[3] and independently by Huang, Caines, and Malhamé [4]. In this paper, we consider second-order time-dependent MFG systems of the form \[\label{eq:MFG-system} \begin{align} -\partial_t u -\nu \Delta u + H(t,x,\nabla u) &= F[m] \;&& \mathrm{in }~ (0,T) \times \Omega, \\ \partial_t m -\nu \Delta m - \mathrm{div}\left( H_p(t,x,\nabla u) ~m \right) &= G \;&& \mathrm{in}~ (0,T) \times \Omega, \\ m(0) = m_0,\quad u(T) &= S[m(T)] \;&& \mathrm{in}~ \Omega, \\ m = 0, \quad u&=0 \;&&\mathrm{on}~ (0,T) \times \partial \Omega , \end{align}\tag{1}\] where \(\Omega\subset \mathbb{R}^d\) is a bounded domain that represents the state space of the game, where \(u\) denotes the value function, and where \(m\) denotes the density of players across the state space \(\Omega\). The main assumptions in this work regarding the problem data are that \(\nu>0\) is a constant, the Hamiltonian \((t,x,p)\mapsto H(t,x,p)\) is \(C^{1,1}\) with respect to \(p\), the coupling terms \(F\) and \(S\) are Lipschitz continuous and monotone operators on suitably chosen function spaces, the source term \(G\) is nonnegative in the sense of distributions, and the initial density \(m_0\) is nonnegative. More detailed assumptions on the problem data are stated in Section 2 below.

The numerical approximation of MFG presents several significant challenges, such as the nonlinear coupling of the forward and backward parabolic equations, the lack of coercivity of the spatial differential operators, and the need to preserve nonnegativity of the density at the discrete level. Finite element methods (FEM) for time-dependent MFG systems, such as 1 , and their steady-state counterparts, were analysed recently by Osborne and the first author in [5], [6], where the convergence of the methods was proved for MFG systems with nondifferentiable Hamiltonians. The FEM in these works are based on a continuous piecewise affine discretization with stabilization of the advective terms in space, and an implicit Euler discretization of the temporal derivatives with lumped masses, which are chosen to ensure a discrete maximum principle and the nonnegativity of the approximations of the density \(m\). The asymptotic near quasi-optimality and optimal convergence rates of the FEM were then proved in [7] for steady-state problems where the Hamiltonian \(H\) is \(C^{1,1}\) with respect to the gradient variable of the value function, and the couplings are monotone. Convergence rates of FEM for steady-state MFG where the Hamiltonian \(H\) is merely Lipschitz have also recently been shown in [8], based on the regularization analysis from [9]. Berry, Ley & Silva [10] have analysed convergence rates of FEM when the solutions are stable with respect to linearizations. This approach has recently been extended to spatial semidiscretizations for time-dependent problems by Berry in [11].

Whereas the works above concern the a priori analysis of FEM for MFG, we are interested here in a posteriori analysis, where the goal is to bound the error in terms of a computable error estimator that depends on the numerical solution and the computational meshes, without requiring a priori knowledge of the true solution or assumptions about its regularity. A posteriori error estimators are a central ingredient for adaptive mesh refinement strategies, which can be significantly more computationally efficient than uniform mesh refinements. We refer the reader to the textbook [12] for an introduction to a posteriori error analysis. We emphasize also that the a posteriori error analysis of parabolic problems features a number of challenges not present in the elliptic setting, especially regarding the efficiency of the estimators and the lack of temporal conformity of standard temporal discretizations, see [13]–[21] for a variety of approaches in the case of linear parabolic problems, and also to [22] for an introductory overview. The only work so far on a posteriori error analysis for numerical discretizations of MFG is [23], where reliable and efficient estimators were established for a broad class of stabilized FEM in the case of steady-state MFG systems with \(C^{1,1}\) Hamiltonians. The estimators consist of residual-based estimators plus an additional stabilization estimator that results from the stabilization terms in the method. It was also shown in [23] that, if the stabilization terms have some additional structure, including being locally affine-preserving, then the stabilization estimator can be bounded in terms of the jump estimators, so it does not need to be computed in practice.

The starting point for the a posteriori error analysis of 1 is the well-posedness of the continuous problem. It is known already that, under suitable assumptions such as the monotonicity of the coupling operators \(F\) and \(S\) and nonnegativity of \(m_0\) and \(G\), the MFG system is well-posed with a unique weak solution pair \((u,m)\in X\times Y\) where the Bochner–Sobolev spaces are defined by \(X\mathrel{\vcenter{:}}= L^2(0,T;H^1_0(\Omega))\) and \(Y\mathrel{\vcenter{:}}= X\cap H^1(0,T;H^{-1}(\Omega))\); we refer to [6] and [24], as well as Section 2 below, for details. In the first main result of this work, see Theorem 1 below, we prove the stronger quantitative property, namely the equivalence between error and dual norm of residual, which takes the form \[\left\lVert u-\overline{u}\right\rVert_X + \left\lVert m-\overline{m}\right\rVert_Y \eqsim \mathcal{R}(\overline{u},\overline{m}), \label{eq:intro-equiv}\tag{2}\] for all \((\overline{u},\overline{m})\in X\times Y\) such that \(\overline{m}\geq 0\) a.e.in \((0,T)\times \Omega\), where the residual functional \(\mathcal{R}(\cdot,\cdot)\), defined in 25 below, comprises the dual norms of the residuals of the weak forms of the HJB and KFP equations, see 22 , alongside the initial condition error of the KFP equation. The notation \(\eqsim\) in 2 means that the left- and right-hand sides are bounded from above and below by each other, up to hidden multiplicative constants that depend solely on the problem data. We note that one improvement in the analysis here for the time-dependent MFG system over that of the steady-state problem in [23] is that the equivalence result 2 is not restricted to functions \((\overline{u},\overline{m})\) in some neighbourhood of the solution, so the resulting a posteriori error bounds for time-dependent problems are applicable also on coarse meshes.

We then derive a posteriori error bounds for a broad family of discretizations of 1 , consisting of conforming piecewise affine FEM in space and implicit Euler discretizations in time, along with general abstract stabilizations to ensure nonnegativity of the density. This class includes the method of [6] as a particular example. Methods in this family lead to numerical approximations \((u _{h,\tau},m _{h,\tau})\in \mathbb{V}_1\times \mathbb{V}_2\), where the discrete spaces \(\mathbb{V}_i\) are discrete approximations spaces of functions that are piecewise constant in time with respect to the time-step partition of \([0,T]\), and continuous piecewise affine in space with respect to a mesh \(\mathcal{T}\) over \(\Omega\), see Section 4 below for further details. Since the approximation spaces \(\mathbb{V}_i\) are contained in the space \(X\) but not in \(Y\), the equivalence 2 bound cannot be immediately applied to the numerical solution \((u _{h,\tau},m _{h,\tau})\). Following [14], we address this challenge by defining a suitable extension of the norm on \(Y\) to the sum space \(\mathbb{V}_i+Y\), see 40 below, and we apply the equivalence result 1 to the pair \((u _{h,\tau},M _{h,\tau})\) where \(M _{h,\tau}\in Y\) is a suitably defined continuous piecewise affine in time reconstruction based on \(m _{h,\tau}\). These extended norms and reconstructions bridge the temporal conformity gap, enabling us to use the continuous equivalence 2 to bound the errors of the stabilized FEM solution \((u _{h,\tau},m _{h,\tau}) \in \mathbb{V}_1\times \mathbb{V}_2\). This leads to our main results on the a posteriori error bounds in Theorems 2, 3, and 4 below, where we show the reliability and efficiency of the estimators via bounds of the form \[\label{eq:intro:a-post-bounds} \begin{align} \left\lVert u-u _{h,\tau}\right\rVert_{\mathbb{V}_1+ Y} + \left\lVert m-m _{h,\tau}\right\rVert_{\mathbb{V}_2+ Y} \lesssim \eta(u _{h,\tau},m _{h,\tau}) + \text{ oscillation}, \\ \eta(u _{h,\tau},m _{h,\tau}) \lesssim \left\lVert u-u _{h,\tau}\right\rVert_{\mathbb{V}_1+ Y} + \left\lVert m-m _{h,\tau}\right\rVert_{\mathbb{V}_2+ Y} + \text{ oscillation}, \end{align}\tag{3}\] where \(\eta(u _{h,\tau},m _{h,\tau})\) is a computable estimator that includes the standard residual estimators, the temporal jump estimators that measure the lack of \(Y\)-conformity of the numerical approximations, stabilization estimators, and estimators for the error in the initial and final time conditions. As above, the notation \(\lesssim\) means that the inequality holds with a hidden multiplicative constant that is independent of the discretization parameters. The oscillation terms in 3 represent the data approximation error terms that typically arise in residual-type a posteriori error analysis. The estimator and data oscillation terms are defined in Section 5.1. Therefore, the bounds in 3 show that the resulting estimators are reliable, as well as globally efficient.

For the general class of abstract stabilizations for which we show 3 , the stabilization estimator is not locally computable. Therefore, in order to obtain a locally computable and locally efficient estimator, we refine the analysis for a narrower class of stabilizations, namely those that consist of mass lumping of the time derivative terms with the patchwise affine-preserving stabilization of the spatial terms. First, we show that the part of the stabilization estimator that results from mass lumping is bounded by the standard temporal jump estimator, see Proposition 5 below, under the hypothesis that \(h^2\lesssim \tau\) where \(h\) is the spatial mesh-size and \(\tau\) is the time-step size. Note that the condition \(h^2\lesssim \tau\) corresponds to the practically relevant case in computations, as the time-steps are usually not excessively small compared to the mesh-size. Furthermore, we show that, for patchwise affine-preserving stabilizations, the spatial stabilization terms are bounded by the spatial jump estimator (which is one of the components of the standard residual estimator), see Proposition 6 below. Consequently, we show in Theorem 7 below that the stabilization estimators can be removed, and the whole error is bounded (up to data oscillation terms) by the residual, temporal jump, initial, and final time estimators. Therefore, we obtain reliable and locally efficient a posteriori error estimators for this class of methods.

The remainder of this paper is organized as follows. Section 2 introduces the notation and fundamental properties of the MFG system 1 required for our analysis. Section 3 establishes the continuous-level equivalence between the approximation errors and the dual norms of the residuals. The general class of stabilized FEM schemes considered in this work is defined in Section 4. In Section 5, we construct the a posteriori error estimators and prove their reliability and efficiency. Finally, Section 6 derives stabilization-free error bounds for a class of stabilization schemes with additional structure.

2 Preliminaries↩︎

2.0.0.1 Basic notation

We denote \(\mathbb{N}\mathrel{\vcenter{:}}= \{1,2,3,\cdots\}\). For a Lebesgue measurable set \(\omega\subset \mathbb{R}^d\), \(d\in\mathbb{N}\), we let \(|\omega|_d\) denote its Lebesgue measure and \(\mathop{\mathrm{diam}}~ \omega\) its diameter. Let \(L^2(\omega)\) and \(L^2(\omega;\mathbb{R}^d)\) denote respectively the usual Lebesgue spaces of scalar, respectively \(d\)-dimensional vector fields, on \(\omega\). The inner products on \(L^2(\omega)\) and \(L^2(\omega;\mathbb{R}^d)\) are denoted commonly by \((\cdot,\cdot)_{\omega}\), with induced norm denoted by \(\left\lVert\cdot\right\rVert_\omega\); there will be no risk of confusion as the two cases will be distinguished by the arguments. Let \(L^{\infty}(\omega;\mathbb{R}^{d\times d})\) denote the space of essentially bounded \(d\times d\) matrix-valued functions on \(\omega\), equipped with the essential supremum norm \(\|\cdot\|_{L^{\infty}(\omega;\mathbb{R}^{d\times d})}\) that is induced by the Frobenius norm on \(\mathbb{R}^{d\times d}\). Let \(\Omega\) be a bounded, connected, open subset of \(\mathbb{R}^d\), \(d\in\mathbb{N}\), with Lipschitz boundary \(\partial \Omega\). For a given final time \(T>0\), let \(I \mathrel{\vcenter{:}}= (0,T)\) denote the bounded open time interval, and let \(Q_T = I \times \Omega\) denote the time-space domain. Let \(H^1(\Omega)\) and \(H^1_0(\Omega)\) denote the usual Sobolev spaces, cf. [25]. Note that the Poincaré inequality for the bounded domain \(\Omega\) implies that \(v\mapsto \left\lVert\nabla v\right\rVert_\Omega\) defines a norm on \(H^1_0(\Omega)\). Let \(H^{-1}(\Omega)\) denote the dual space of \(H_0^1(\Omega)\). The canonical norm on \(H^{-1}(\Omega)\) is defined by \[\left\lVert\Phi\right\rVert_{H^{-1}(\Omega)}\mathrel{\vcenter{:}}= \sup_{\substack{v\in H^1_0(\Omega)\\\left\lVert\nabla v\right\rVert_\Omega =1}} \langle \Phi, v\rangle_{} \quad \forall\Phi \in H^{-1}(\Omega),\] where \(\langle\cdot,\cdot\rangle_{}\) denotes the duality pairing between \(H^{-1}(\Omega)\) and \(H^1_0(\Omega)\).

2.0.0.2 Notation for inequalities

In the following, the notation \(x \lesssim y\) will indicate the existence of a constant \(C >0\) such that \(x \leq C y\). If there exist two constants such that \(x \lesssim y\) and \(y \lesssim x\) we will denote this relation by \(x \eqsim y\). We shall, when necessary, list the dependency of the hidden constant but not explicitly define it. Hidden constants will typically depend on the PDE data, the domain, and the quality of the domain discretization. Hidden constants will not depend on any discretization parameters.

2.0.0.3 Bochner spaces

In the following, we will consider several Bochner–Sobolev spaces [26]. Let \(X\mathrel{\vcenter{:}}= L^2(I;H^1_0(\Omega))\), which is a Hilbert space when equipped with the inner-product \((w,v)\mapsto\int_0^T(\nabla w, \nabla v)_\Omega\,\mathrm{d} t\) and norm \[\left\lVert w\right\rVert_X^2 \mathrel{\vcenter{:}}= \int_0^T \left\lVert\nabla w\right\rVert^2_\Omega\, \mathrm{d}t\quad\forall w\in X. \label{eq:X-norm-def}\tag{4}\] Let \(X^*\) denote the dual space of \(X\) and let \(\left\langle\cdot,\cdot\right\rangle_{X^* \times X}\) denote the natural duality pairing. We define the dual norm on \(X^*\) by \[\left\lVert f\right\rVert_{X^*} \mathrel{\vcenter{:}}= \sup_{w \in X \backslash \{0\}} \frac{\left\langle f,w\right\rangle_{X^* \times X}}{\left\lVert w\right\rVert_X} \quad \forall f \in X^*.\] It is known that the space \(X^*\) can be identified with the space \(L^2(I; H^{-1}(\Omega))\) [26].

We also define the space \(Y\mathrel{\vcenter{:}}= L^2(I;H^1_0(\Omega))~\cap~H^1(I;H^{-1}(\Omega))\), which consists of all functions \(v\in X\) that are weakly differentiable in time with time derivative \(\partial_t v \in L^2(I;H^{-1}(\Omega))\). It is known that the space \(Y\) can be equipped with the norm \(\left\lVert\cdot\right\rVert_Y\) defined by \[\label{eq:Y-norm-def} \left\lVert v\right\rVert_Y^2 \mathrel{\vcenter{:}}= \int_0^T \left\lVert\partial_t v\right\rVert_{H^{-1}(\Omega)}^2 + \left\lVert\nabla v\right\rVert^2_\Omega\;\mathrm{d}t+ \left\lVert v(T)\right\rVert_\Omega^2 \quad \forall v \in Y,\tag{5}\] where \(v(T)\) denotes the trace value of \(v\) at time \(T\) in \(L^2(\Omega)\). In particular, the norm \(\left\lVert\cdot\right\rVert_Y\) is well-defined since the space \(Y\) is continuously embedded into \(C([0,T]; L^2(\Omega))\). Furthermore, we have the embedding bound \[\max_{t \in [0,T]} \left\lVert v(t)\right\rVert_\Omega \leq \left\lVert v\right\rVert_Y \quad \forall v \in Y. \label{eq:Y-embedding}\tag{6}\] Additionally, we define \(Y_0\) as the space of functions in \(Y\) that vanish at initial time, i.e.\(Y_0 \mathrel{\vcenter{:}}= \{ v \in Y\;:\;v(0)=0 \}\).

2.0.0.4 Problem data

Let the diffusion coefficient \(\nu >0\) be constant. The initial distribution \(m_0\) is assumed to be in \(L^\infty(\Omega)\) with \(m_0 \geq 0\) a.e.in \(\Omega\).

The Hamiltonian of the underlying optimal control problem is defined by \[H(t,x,p) \mathrel{\vcenter{:}}= \sup_{\alpha \in \mathcal{A}} \left[ b(t,x,\alpha) \cdot p - f(t,x,\alpha) \right] \quad \forall (t,x,p) \in \overline{Q_T} \times \mathbb{R}^d, \label{eq:H-def}\tag{7}\] where it is assumed that the control set \(\mathcal{A}\) is a compact metric space, and the control-dependent drift \(b:[0,T] \times \Omega \times \mathcal{A} \rightarrow \mathbb{R}^d\) and control-dependent running cost \(f:Q_T \times \mathcal{A} \rightarrow \mathbb{R}\) are uniformly continuous on \(\overline{Q_T} \times \mathbb{R}^d\). It follows from these assumptions that the Hamiltonian \(H\), defined in 7 above, is convex and Lipschitz continuous in its third argument, i.e. \[\label{eq:H-cts-condition} \left\lvert H(t,x,p) - H(t,x,q) \right\rvert \leq L_H \left\lvert p-q \right\rvert \quad \forall (t,x,p,q) \in \overline{Q_T} \times \mathbb{R}^d \times \mathbb{R}^d,\tag{8}\] where \(L_H\mathrel{\vcenter{:}}= \left\lVert b\right\rVert_{C\left(\overline{Q_T} \times \mathcal{A}; \mathbb{R}^d \right)}\). We further assume that \(H\) is \(C^{1,1}\) on \(\mathbb{R}^d\) with respect to its third argument, i.e.\(H_p\) exists for all arguments, and there exists a constant \(L_{H_p}\) such that \[\label{eq:dpH-cts-condition} \left\lvert H_p(t,x,p) - H_p(t,x,q) \right\rvert \leq L_{H_p} \left\lvert p-q \right\rvert \quad \forall (t,x,p,q) \in \overline{Q_T} \times \mathbb{R}^d \times \mathbb{R}^d.\tag{9}\] The Lipschitz continuity of \(H\) in 8 implies that \[\label{eq:dpH-bounded-condition} \left\lvert H_p(t,x,p) \right\rvert \leq L_H \quad \forall (t,x,p) \in \overline{Q_T} \times \mathbb{R}^d.\tag{10}\] Let \(D_H(\cdot,\cdot): \overline{Q_T} \times \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}\) denote the Bregman divergence of the Hamiltonian with respect to its third argument, defined as \[\label{eq:Bregman-div-def} D_H(t,x,p,q) \mathrel{\vcenter{:}}= H(t,x,p) - H(t,x,q) - H_p(t,x,q) \cdot (p-q) \quad \forall (t,x,p,q)\in \overline{Q_T} \times \mathbb{R}^d \times \mathbb{R}^d.\tag{11}\] As a consequence of the convexity of the Hamiltonian and the Lipschitz continuity condition 9 , we have the following well-known inequality, which can be found for instance in [27], \[\label{eq:Bregman-inequality} \left\lvert H_p(t,x,p)-H_p(t,x,q) \right\rvert^2 \leq 2 L_{H_p} D_H(t,x,p,q) \quad \forall (t,x,p,q) \in \overline{Q_T}\times \mathbb{R}^d\times \mathbb{R}^d.\tag{12}\] To abbreviate the notation, for a function \(v\in X\), let \(H[\nabla v]\in L^2(0,T;L^2(\Omega))\) denote the composition of \(H\) with \(\nabla v\), i.e.\(H[\nabla v](t,x)\mathrel{\vcenter{:}}= H(t,x,\nabla v(t,x))\) a.e.\(t\in (0,T)\), a.e.\(x\in \Omega\). In a similar manner, for a pair \((w,v)\in X\times X\), let \(H_p[\nabla v]\) and \(D_H[\nabla w,\nabla v]\in L^2(0,T;L^2(\Omega))\) be defined by \(H_p[\nabla v](t,x)\mathrel{\vcenter{:}}= H_p(t,x,\nabla v(t,x))\) and \(D_H[\nabla w,\nabla v](t,x)\mathrel{\vcenter{:}}= D_H(t,x,\nabla w(t,x),\nabla v(t,x))\) for a.e.\(t\in (0,T)\), a.e.\(x\in \Omega\).

We now specify the assumptions on the coupling operators \(F\) and \(S\). We suppose that there exists a Hilbert space \(Z\) such that the nonlinear coupling operator \(F\colon Z \to X^*\) is Lipschitz continuous, i.e.there exists a constant \(L_F >0\) such that \[\begin{align} \left\lVert F[v_1] - F[v_2]\right\rVert_{X^*} &\leq L_F \left\lVert v_1 -v_2\right\rVert_{Z} \quad \forall v_1,v_2 \in Z. \label{eq:F-cts-condition} \end{align}\tag{13}\] Similarly, we assume that the final time coupling operator \(S:L^2(\Omega) \to L^2(\Omega)\) is Lipschitz continuous, i.e.there exists a constant \(L_S>0\) such that \[\left\lVert S[v_1] - S[v_2]\right\rVert_\Omega \leq L_S \left\lVert v_1 - v_2\right\rVert_\Omega \label{eq:S-cts-condition} \quad \forall v_1,v_2 \in L^2(\Omega).\tag{14}\] Furthermore, we assume that the spaces \(X\) and \(Y\) are continuously embedded in \(Z\), with the embedding of \(Y \hookrightarrow Z\) moreover compact. We assume that the coupling operators \(F\) and \(S\) satisfy a monotonicity condition: \[\begin{align} \label{eq:F38S-monotonicity-condition} 0 \leq \int^T_0 \left\langle F[v_1]-F[v_2],v_1-v_2\right\rangle_{} \mathrm{d}t+ \left(S[v_1(T)]-S[v_2(T)],v_1(T)-v_2(T)\right)_{\Omega}, \end{align}\tag{15}\] for all \(v_1, v_2 \in Y\). Note the monotonicity condition 15 is only required on the space \(Y\), even for coupling terms defined over a larger Hilbert space \(Z\).

Let \(G \in L^2(I; H^{-1}(\Omega))\) be of the form \(G = g_0 - \nabla \cdot g_1\), where \(g_0 \in L^{r}(I;L^{s}(\Omega))\) and \(g_1 \in L^{2r}(I;L^{2s}(\Omega;\mathbb{R}^d))\) for some indices \(r,s \in (1,\infty]\) satisfying \[\label{eq:G-indices-conditon} \frac{1}{r} + \frac{d}{2s} < 1.\tag{16}\] Furthermore, we assume that the source term \(G\) is nonnegative in the sense of distributions, i.e.\(\int_0^T \left\langle G,w\right\rangle_{} \mathrm{d}t\geq 0\) for all \(w \in X\) which satisfy \(w \geq 0\) a.e. in \(Q_T\).

2.0.0.5 Weak formulation

A weak formulation of the MFG system 1 is as follows: find \((u,m) \in X \times Y\) such that \(m(0) = m_0\) and \[\tag{17} \begin{align} & \int_0^T \left[ \left\langle\partial_t v,u\right\rangle_{} + \left(\nu \nabla u,\nabla v\right)_{\Omega} + \left(H[\nabla u],v\right)_{\Omega} \right]\mathrm{d}t\tag{18} \\ & =\int_0^T \left\langle F[m],v\right\rangle_{}\mathrm{d}t+ \left(S[m(T)],v(T)\right)_{\Omega}, \notag \\ &\int_0^T \left[ \left\langle\partial_t m,\phi\right\rangle_{} + \left(\nu \nabla m,\nabla \phi\right)_{\Omega} + \left(m H_p[\nabla u],\nabla \phi\right)_{\Omega} \right]\mathrm{d}t = \int_0^T \left\langle G,\phi\right\rangle_{}\mathrm{d}t, \tag{19} \end{align}\] for all \((v,\phi) \in Y_0 \times X\). Notice that the weak formulation 17 is obtained by integration-by-parts in time of the temporal derivative term in the HJB equation. Under the hypotheses on the problem data given above, in [6], the existence and uniqueness of a weak solution of 17 was shown for the case when \(Z = L^2(I; L^2(\Omega))\), without requiring differentiability of \(H\). Following the proof of [24], it is clear that existence of a solution also holds when \(F\) is defined on the more general space \(Z\) that satisfies the assumptions stated above. Note that the uniqueness of the solution follows from the nonnegativity of \(m_0\) and \(G\), and the monotonicity condition 15 on the coupling terms \(F\) and \(S\). Furthermore, the density \(m\) is nonnegative a.e. in \(Q_T\), and is additionally essentially bounded, i.e. \[\begin{align} \left\lVert m\right\rVert_{L^\infty(Q_T)} &\leq M_\infty, \label{eq:M95infty} \end{align}\tag{20}\] where \(M_\infty\) is a constant that depends only \(d\), the indices \(r\) and \(s\) that satisfy 16 , on \(\nu\), \(L_H\), the measure of \(\Omega\), the \(L^r(I;L^q(\Omega))\)-norm of \(g_0\), and the \(L^{2r}(I;L^{2q}(\Omega;\mathbb{R}^d))\)-norm of \(g_1\). A proof of 20 can be found in [28]. Finally, we note that the hypotheses that \(F\) and \(S\) respectively take values in \(X^*\) and in \(L^2(\Omega)\) imply that the value function \(u\) that solves 18 is also in the space \(Y\), i.e.the temporal derivative \(\partial_t u\) exists in \(L^2(I;H^{-1}(\Omega))\), and \(u\) satisfies \(u(T)=S[m(T)]\) in \(L^2(\Omega)\) and \[\label{eq:HJB95time95strong95form} \int_0^T \left[-\left\langle\partial_t u,v\right\rangle_{} + (\nu \nabla u,\nabla v)_\Omega + (H[\nabla u],v)_\Omega \right]\mathrm{d}t = \int_0^T \left\langle F[m],v\right\rangle_{} \mathrm{d}t \quad \forall v \in X.\tag{21}\] Both formulations 18 and 21 of the HJB equation will be used in the following analysis.

3 Equivalence between errors and residuals↩︎

Our first contribution toward deriving computable a posteriori error bounds is to prove, in Theorem 1 below, the equivalence between the norms of the difference between the true solution \((u,m)\) and some general \((\overline{u},\overline{m})\) and the norms of the residual operators of the MFG system.

We now define the residual operators associated to the HJB and KFP equations. Let \(R_1^X\colon X\times Y_0 \to Y_0^*\) and \(R_2^Y\colon X\times Y_0\to X^*\) be defined by \[\tag{22} \begin{align} \left\langle R_1^X(\overline{u},\overline{m}),v\right\rangle_{Y^*_0 \times Y_0} &\mathrel{\vcenter{:}}= \int_0^T \left\langle F[\overline{m}],v\right\rangle_{} \mathrm{d}t+ \left(S[\overline{m}(T)],v(T)\right)_{\Omega} \tag{23}\\ & -\int_0^T \left[ \left\langle\partial_t v,\overline{u}\right\rangle_{} + \left(\nu \nabla \overline{u},\nabla v\right)_{\Omega} + \left(H[\nabla \overline{u}],v\right)_{\Omega} \right]\mathrm{d}t, \notag\\ \left\langle R_2^Y(\overline{u},\overline{m}),w\right\rangle_{X^* \times X} &\mathrel{\vcenter{:}}= \int_0^T \left[ \left\langle G-\partial_t \overline{m},w\right\rangle_{} -\left(\nu \nabla \overline{m}+ \overline{m}H_p[\nabla \overline{u}],\nabla w\right)_{\Omega} \right]\mathrm{d}t, \tag{24} \end{align}\] for all \((v,w) \in Y_0 \times X\). Let \(\mathcal{R}(\cdot,\cdot):X \times Y \to \mathbb{R}_{\geq 0}\) denote the total residual norm functional defined by \[\label{eq:calR-def} \mathcal{R}(\overline{u},\overline{m}) \mathrel{\vcenter{:}}=\left\lVert R_1^X(\overline{u},\overline{m})\right\rVert_{Y^*_0}+\left\lVert R_2^Y(\overline{u},\overline{m})\right\rVert_{X^*}+ \left\lVert\overline{m}(0)-m_0\right\rVert_\Omega,\tag{25}\] for any pair \((\overline{u},\overline{m}) \in X \times Y\). It is clear from the definitions of \(\mathcal{R}(\cdot,\cdot)\) and the weak formulation in 17 that the weak solution \((u,m) \in X \times Y\) satisfies \(\mathcal{R}(u,m) = 0\). Now we present the main result of this section.

The proof is deferred to Section 3.2 below.

3.1 Stability of the HJB and KFP equations↩︎

Here we establish stability results between the dual norms of the residuals defined in 22 and the approximation error norms.

Next, we record some inf-sup stability results that are useful for proving the stability of the HJB and KFP equations.

We omit the proof of Lemma 2, since ?? is easily deduced from the inf-sup stability of the heat equation, see e.g. [14], and from Gronwall’s inequality, while ?? follows from ?? by duality. The next Lemma shows the stability of the HJB and KFP when considered separately.

As will become clear in the proof of Theorem 1, it is important that we do not further bound the second term in the upper bound of ?? using the continuity of \(H_p\) and the boundedness of \(m\).

3.2 Proof of the equivalence result↩︎

4 Stabilized finite element methods↩︎

In this section, we present the numerical framework of the stabilized FEM used to approximate the MFG system 1 .

4.1 Spatial discretization↩︎

To avoid unnecessary technicalities, we shall only consider a fixed, time-independent spatial discretization of \(\Omega\). Let \(\mathcal{T}\) denote a conforming simplicial mesh of the domain \(\Omega\). We adopt the convention that each element \(K\in\mathcal{T}\) is closed. For each element \(K\in\mathcal{T}\), we let \(h_K\) denote the diameter of \(K\) and \(\varrho_K\) denote the largest diameter of an inscribed ball in \(K\). Let \(\theta_\mathcal{T}\mathrel{\vcenter{:}}= \max_{K\in\mathcal{T}} h_K/\rho_K\) denote the shape-regularity parameter of the mesh \(\mathcal{T}\). Let \(\mathcal{V}\) denote the set of all vertices of \(\mathcal{T}\) and let \(\mathcal{V}^I=\mathcal{V}\cap \Omega\) denote the set of interior vertices of \(\mathcal{T}\). Two distinct vertices are called neighbours if they belong to a common element of \(\mathcal{T}\). Let \(\mathcal{E}\) denote the set of edges of the mesh \(\mathcal{T}\), i.e.the set of all closed line segments formed by all pairs of neighbouring vertices. Given a simplex \(K\in\mathcal{T}\) we denote the set of edges contained in \(K\) by \(\mathcal{E}_K\). Let \(\mathcal{F}\) denote the set of all faces of \(\mathcal{T}\), and let \(\mathcal{F}^I\) denote the subset of all interior faces of \(\mathcal{F}\), i.e.all faces that are not contained in \(\partial\Omega\). Note that for \(d=2\), edges and faces of the mesh coincide. For each face \(F\in\mathcal{F}\), let \(h_F\) denote the diameter of \(F\). For each element \(K\in\mathcal{T}\), we let \(\mathcal{F}^I_K\) denote the set of interior faces \(F\in\mathcal{F}^I\) that are contained in \(K\). For each face \(F\in\mathcal{F}\), we let \(\mathcal{T}_{F}\) denote the set of elements that contain \(F\), and we let \(\omega_F \mathrel{\vcenter{:}}= \bigcup_{K\in\mathcal{T}_{F}}K\) denote its associated patch. For each element \(K\in\mathcal{T}\), we define the set \(\widetilde{\mathcal{T}}_K\) of face-sharing neighbouring elements and the associated patch \(\omega_K\) \[\begin{align} \widetilde{\mathcal{T}}_K\mathrel{\vcenter{:}}= \bigcup_{F\in\mathcal{F}_K^{I}}\mathcal{T}_{F}, &&& \omega_K\mathrel{\vcenter{:}}= \bigcup_{K^\prime\in\widetilde{\mathcal{T}}_K} K^\prime = \bigcup_{F\in\mathcal{F}_K^{I}}\omega_F. \end{align}\] Note that \(K\in\widetilde{\mathcal{T}}_K\) for all \(K\in\mathcal{T}\).

Let \(V_h \subset H^1_0(\Omega)\) denote the linear Lagrange finite element space constructed over \(\mathcal{T}\) defined by \[V_h \mathrel{\vcenter{:}}= \{ v_h \in H^1_0(\Omega)~ :~ v_h |_{K} \in \mathcal{P}_1(K) \quad \forall K \in \mathcal{T}\}, \label{eq:fem-space-def}\tag{36}\] where \(\mathcal{P}_1(K)\) denotes the space of linear polynomials on \(K\). Let \(\{\varphi_z \}_{z \in \mathcal{V}^I}\) denote the set of standard Lagrangian basis functions characterised by \(\varphi_z(z) = 1\) and \(\varphi_z(z')=0\) for all \(z' \in \mathcal{V}\backslash \{ z \}\).

4.2 Space-time finite element spaces↩︎

Let \(\mathcal{J} \mathrel{\vcenter{:}}= \{ I_n \}_{n=1}^N\) denote a partition of \(\overline{I} = [0,T]\) into intervals \(I_n \mathrel{\vcenter{:}}= (t_{n-1},t_n)\) with \(0=t_0 \leq t_{n-1} < t_n \leq t_N\) for each \(n \in \{1,...,N\}\). For each \(n \in \{1,...,N\}\), the time-step size is denoted by \(\tau_n \mathrel{\vcenter{:}}= t_n - t_{n-1}\). Let \(\mathbb{V}_{\mathcal{T},\mathcal{J}}\) denote the space of \(V_h\)-valued functions defined on \([0,T]\) that are piecewise constant with respect to \(\mathcal{J}\). In other words, \(\mathbb{V}_{\mathcal{T},\mathcal{J}}\) is defined by \[\mathbb{V}_{\mathcal{T},\mathcal{J}}\mathrel{\vcenter{:}}= \left\{ v _{h,\tau}\in X\,:\, v _{h,\tau}|_{I_n} \in V_h \text{ is constant } \forall n \in \{1,...,N \} \right\}.\] The implicit Euler time stepping method naturally leads to numerical solutions that can be viewed as pairs of functions in \(\mathbb{V}_{\mathcal{T},\mathcal{J}}\). However, given the forward-backward structure of the equations of the MFG system with initial and final time conditions on the different solution components, it will be natural to view the numerical solutions as functions that are defined everywhere on \([0,T]\), with different left- and right- continuity properties between the value function and density, owing to the different time directions of the equations. For a function \(v _{h,\tau}\in \mathbb{V}_{\mathcal{T},\mathcal{J}}\), let \(v _{h,\tau}(t_n^-)\) and \(v _{h,\tau}(t_n^+)\) denote respectively the left- and right-limits of \(v _{h,\tau}\) at time \(t_n\). Also note that \(v _{h,\tau}(t_n^+) = v _{h,\tau}|_{I_{n+1}}\) and \(v _{h,\tau}(t_n^-) = v _{h,\tau}|_{I_n}\) for all \(n \in \{0,...,N-1\}\), and \(v _{h,\tau}(0^+) = v _{h,\tau}|_{I_1}\) and \(v _{h,\tau}(T^-) = v _{h,\tau}|_{I_N}\). We say that \(v _{h,\tau}\) is left-continuous if \(v _{h,\tau}(t_n) = v _{h,\tau}(t^-_n)\) for all \(n \in \{1,...,N \}\) and right-continuous if \(v _{h,\tau}(t_n) = v _{h,\tau}(t^+_n)\) for all \(n \in \{0,...,N-1 \}\). We define the right- and left- continuous finite element spaces by \[\label{eq:signed95FEM95spaces} \begin{align} \mathbb{V}_1&\mathrel{\vcenter{:}}= \{v _{h,\tau}:[0,T] \to V_h,~ v _{h,\tau}|_{(0,T)} \in \mathbb{V}_{\mathcal{T},\mathcal{J}}~:~ v _{h,\tau}\text{ is right-continuous} \}, \\ \mathbb{V}_2&\mathrel{\vcenter{:}}= \{ v _{h,\tau}:[0,T] \to V_h,~ v _{h,\tau}|_{(0,T)} \in \mathbb{V}_{\mathcal{T},\mathcal{J}}~:~ v _{h,\tau}\text{ is left-continuous} \}. \end{align}\tag{37}\] We emphasize that functions in \(\mathbb{V}_i\) are considered to be defined at all points in \([0,T]\), and two functions that agree on up to a null set of \([0,T]\) are not identified. We take this viewpoint in order to handle the initial and final time conditions of the problem in a natural manner. In order to alleviate the notation, the dependence of the spaces on the mesh \(\mathcal{T}\) and time partition \(\mathcal{J}\) will be left implicit. Note that the spaces \(\mathbb{V}_i\) defined above are only generally conforming with respect to the space \(X\), and not \(Y\), since the spaces \(\mathbb{V}_i\) contain functions that may be discontinuous at the timestep points \(t_n \in (0,T)\). Therefore, in order to consider the error between the exact and discrete solutions, we consider the sum space \(\mathbb{V}_i+ Y\) defined by \[\label{eq:space-time-fem-space-defs} \mathbb{V}_i+Y \mathrel{\vcenter{:}}= \left\{ w \in X~:~ \exists (v _{h,\tau},v) \in \mathbb{V}_i\times Y~\mathrm{s.t.}~ w = v _{h,\tau}+ v \right\}.\tag{38}\] Note that the spaces \(\mathbb{V}_i+Y\) are not generally direct sum spaces. For a given function \(w \in \mathbb{V}_i+Y\) and time level \(t_n \in (0,T)\), the temporal jump operators are defined by \[\left( \! \left| w \right| \! \right)_n \mathrel{\vcenter{:}}= w(t_n^-) - w(t_n^+) \quad \forall n\in \{1,\ldots,N-1\}.\] Moreover, the jumps at the initial and final times are defined by \(\left( \! \left| w \right| \! \right)_0 = w(0) - w(0^+)\) and \(\left( \! \left| w \right| \! \right)_N = w(T^-) - w(T)\) respectively. The temporal reconstruction operators \(\mathcal{I}_i : \mathbb{V}_i+ Y \to Y\) are defined by \[\label{eq:Reconstruction-defs} \begin{align} (\mathcal{I}_1 w)(t) &\mathrel{\vcenter{:}}= w(t) - \frac{t-t_{n-1}}{\tau_n} \left( \! \left| w \right| \! \right)_n,~~t \in [t_{n-1},t_n),~ n \in \{1,...,N\},~ w \in \mathbb{V}_1+Y, \\ (\mathcal{I}_2 w)(t) &\mathrel{\vcenter{:}}= w(t) + \frac{t_n-t}{\tau_n} \left( \! \left| w \right| \! \right)_{n-1}, ~~ t \in (t_{n-1},t_n],~ n \in \{1,...,N\},~ w \in \mathbb{V}_2+Y. \end{align}\tag{39}\] It is easy to check that the operators \(\mathcal{I}_i\) are linear, and moreover take values in \(Y\) as a result of the fact that they map into spaces of functions that are continuous and piecewise smooth in time. Also, by the embedding of \(Y \hookrightarrow C([0,T];L^2(\Omega))\), it is clear that \(\mathcal{I}_i w = w\) for any \(w \in Y\), since the vanishing jumps \(\left( \! \left| w \right| \! \right)_n=0\in H^1_0(\Omega)\) vanish for any \(w\in Y\). We equip the spaces \(\mathbb{V}_i+Y\) with an extended \(Y\)-norm defined by \[\label{eq:error-norm-defs} \begin{align} \left\lVert v\right\rVert_{\mathbb{V}_i+Y} \mathrel{\vcenter{:}}= \left\lVert v-\mathcal{I}_i v\right\rVert_X + \left\lVert\mathcal{I}_i v\right\rVert_Y, \quad \forall v \in \mathbb{V}_i+ Y. \end{align}\tag{40}\] The norm defined above is called the extended \(Y\)-norm since it defines a consistent extension of the norm \(\left\lVert\cdot\right\rVert_Y\) of 5 from the space \(Y\) to the sum space \(\mathbb{V}_i+Y\), which is to say that \(\left\lVert v\right\rVert_{\mathbb{V}_i+Y} = \left\lVert v\right\rVert_Y\) for any \(v \in Y\).

4.3 Discretized problem↩︎

As explained in the introduction, we consider here a family of stabilized FEM, which includes for instance the method in [6], based on an implicit Euler discretization in time, and a piecewise affine FEM is used in space. The stabilization is motivated by the need to preserve nonnegativity of the approximate densities, which plays an important role in the uniqueness of numerical solutions and also in the nonnegativity condition appearing in Theorem 1. We now consider abstract stabilization terms of the form \(\mathcal{S}_i(\cdot;\cdot; \cdot): \mathbb{V}_1\times \mathbb{V}_2\times \mathbb{V}_{\mathcal{T},\mathcal{J}}\to \mathbb{R}\), \(i \in \{ 1,2\}\) that are allowed to be nonlinear with respect to first two arguments, but are assumed to be linear with respect to the third argument. Example 1 below details some concrete stabilizations that we have in mind.

The class of FEM that we consider are of the form: find \((u _{h,\tau},m _{h,\tau})\in \mathbb{V}_1\times \mathbb{V}_2\) such that \[\label{eq:FEM-scheme-def} \begin{align} &\int_0^T \left(-\partial_t U _{h,\tau},v _{h,\tau}\right)_{\Omega} + \left(\nu \nabla u _{h,\tau},\nabla v _{h,\tau}\right)_{\Omega} + \left(H[\nabla u _{h,\tau}],v _{h,\tau}\right)_{\Omega}\mathrm{d}t\\ & + \mathcal{S}_1(u _{h,\tau};m _{h,\tau};v _{h,\tau}) = \int_0^T \left\langle F[m _{h,\tau}],v _{h,\tau}\right\rangle_{} \mathrm{d}t\notag \\ &\int_0^T \left(\partial_t M _{h,\tau},v _{h,\tau}\right)_{\Omega} + \left(\nu \nabla m _{h,\tau},\nabla v _{h,\tau}\right)_{\Omega} + \left(m _{h,\tau}H_p[\nabla u _{h,\tau}],\nabla v _{h,\tau}\right)_{\Omega}\mathrm{d}t\\ & + \mathcal{S}_2(u _{h,\tau};m _{h,\tau};v _{h,\tau}) = \int_0^T \left\langle G,v _{h,\tau}\right\rangle_{} \mathrm{d}t, \notag \end{align}\tag{41}\] for all \(v _{h,\tau}\in \mathbb{V}_{\mathcal{T},\mathcal{J}}\), where we use the shorthand notation \(U _{h,\tau}\mathrel{\vcenter{:}}= \mathcal{I}_1 u _{h,\tau}\) and \(M _{h,\tau}\mathrel{\vcenter{:}}= \mathcal{I}_2 m _{h,\tau}\), and such that \(m _{h,\tau}(0) = \Pi_0 m_0\) and \(u _{h,\tau}(T) = \Pi_T S[m _{h,\tau}(T)]\), where \(\Pi_0\) and \(\Pi_T\) are some quasi-interpolation operators from \(L^2(\Omega)\) to the finite element space \(V_h\) that are chosen by the user to approximate the initial and final time conditions. In practice, the stabilizations \(\mathcal{S}_i\) are chosen to ensure good properties of the FEM 41 , such as existence and uniqueness of the numerical solution, and stability, see [6] for further details. However, the only assumptions that we require for the general a posteriori error analysis of Section 5 below is that a solution of 41 exists, and that discrete densities are nonnegative:

1. if \((u _{h,\tau},m _{h,\tau})\in \mathbb{V}_1\times \mathbb{V}_2\) is a solution of 41 then \(m _{h,\tau}\geq 0\) in \(Q_T\).

Observe that if \(m _{h,\tau}\) is nonnegative in \(Q_T\) then so is \(M _{h,\tau}\mathrel{\vcenter{:}}= \mathcal{I}_2m _{h,\tau}\).

5 A posteriori error bounds↩︎

In this section we present the general analysis of a posteriori error estimators for the class of FEM schemes defined above. For the sake of brevity of exposition, we shall assume in the following that the coupling term \(F\) has images in \(L^2(Q_T)\) for all arguments in \(X\), and that \(G\in L^2(Q_T)\). We refer the reader to [29], [30] and references therein for further details on the treatment of PDE with \(H^{-1}\) source terms.

5.1 The estimators↩︎

We start by defining the estimators that appear in the error bound. The a posteriori estimator \(\eta(u _{h,\tau},m _{h,\tau})\) is defined by \[\label{eq:eta-total-def} \eta(u _{h,\tau},m _{h,\tau}) \mathrel{\vcenter{:}}= \sum_{i=1}^2 \left[ \eta_{\mathrm{J},i} + \eta_{R,i} +\eta_{\mathcal{S},i} \right]+ \eta_0 + \eta_T,\tag{42}\] where the \(\eta_{\mathrm{J},i}\) denote the temporal jump estimators, the \(\eta_{R,i}\) denote the PDE residual estimators, the \(\eta_{\mathcal{S},i}\) denote the stabilization estimators, and \(\eta_0\) and \(\eta_T\) denote the initial and final time condition estimators. These terms are defined respectively in 43 , 48 , 50 and 51 below. To alleviate the notation, the dependency of the above estimators on the numerical solution is left implicit.

5.1.0.1 The temporal jump estimators

The temporal jump estimators \(\eta_{\mathrm{J},i}\) measure the lack of conformity in the space \(Y\) of the functions \(u _{h,\tau}\) and \(m _{h,\tau}\), which are generally discontinuous across time-intervals. Let \(\eta_{\mathrm{J},i}\), \(i\in \{1,2\}\) be defined by \[\eta_{\mathrm{J},i}^2 \mathrel{\vcenter{:}}= \sum_{n=1}^N \sum_{K\in \mathcal{T}} \eta_{\mathrm{J},K,n,i}^2, \quad i \in \{1,2\}. \label{eq:etaRecon-global-defs}\tag{43}\] where the local contributions are given by \(\eta_{\mathrm{J},K,n,1} \mathrel{\vcenter{:}}= \left\lVert\nabla (u _{h,\tau}-U _{h,\tau})\right\rVert_{L^2(I_n \times K)}\) and \(\eta_{\mathrm{J},K,n,2} \mathrel{\vcenter{:}}= \left\lVert\nabla(m _{h,\tau}-M _{h,\tau})\right\rVert_{L^2(I_n \times K)}\) for each time-interval \(I_n\), \(n\in\{1,\dots,N\}\), and each element \(K\in \mathcal{T}\). The temporal jump estimators can be evaluated straightforwardly through the formulas \[\begin{align} [\eta_{\mathrm{J},1}]^2=\frac{1}{3}\sum_{n=1}^N \tau_n \left\lVert\nabla \left( \! \left| u _{h,\tau} \right| \! \right)_n \right\rVert_\Omega^2, &&& [\eta_{\mathrm{J},2}]^2=\frac{1}{3}\sum_{n=1}^N \tau_n \big\lVert \nabla \left( \! \left| m _{h,\tau} \right| \! \right)_{n-1} \big\rVert_\Omega^2 . \end{align}\]

5.1.0.2 The residual estimators

For each time-interval \(I_n\), \(n\in\{1,\dots,N\}\), and element \(K\in \mathcal{T}\), we define the local volume residuals by \[\label{eq:volume-residual-defs} \begin{align} r_{K,n,1} &\mathrel{\vcenter{:}}= F[m _{h,\tau}] + \partial_t U _{h,\tau}+ \nu \Delta (u _{h,\tau}|_K) - H[\nabla u _{h,\tau}], \\ r_{K,n,2} &\mathrel{\vcenter{:}}= G - \partial_t M _{h,\tau}+ \nu \Delta (m _{h,\tau}|_K) + \mathrm{div}(m _{h,\tau}H_p[\nabla u _{h,\tau}]) . \end{align}\tag{44}\] Since the functions \(u _{h,\tau}\) and \(m _{h,\tau}\) are piecewise affine on each element of the mesh, the Laplacian terms in the volume residuals above vanish; however we choose to write them explicitly to emphasize the relation of the volume residuals to the strong form of the original PDE. Also, for each interior face \(F \in \mathcal{F}^I\), we define the local spatial jump residuals by \[\tag{45} \begin{align} j_{F,n,1} &\mathrel{\vcenter{:}}= \nu \left\llbracket \nabla u _{h,\tau}\cdot n_F \right\rrbracket_{F}, \tag{46}\\ j_{F,n,2} &\mathrel{\vcenter{:}}= \nu \left\llbracket \nabla m _{h,\tau}\cdot n_F \right\rrbracket_{F} + m _{h,\tau}\left\llbracket H_p[\nabla u _{h,\tau}] \cdot n_F \right\rrbracket_{F}. \tag{47} \end{align}\] For each \(i\in\{1,2\}\), the total residual estimator \(\eta_{R,i}\) is defined by \[\begin{gather} \eta_{R,i}^2 \mathrel{\vcenter{:}}= \sum_{n=1}^N \sum_{K \in \mathcal{T}} \eta_{R,K,n,i}^2, \tag{48} \\ \eta_{R,K,n,i}^2 \mathrel{\vcenter{:}}= h_K^2 \left\lVert r_{K,n,i}\right\rVert^2_{L^2(I_n \times K)} + \sum_{F \in \mathcal{F}_K^I} h_F \left\lVert j_{F,n,i}\right\rVert_{L^2(I_n \times F)}^2. \tag{49} \end{gather}\]

5.1.0.3 The stabilization estimators

The stabilization estimators are defined by \[\begin{align} \eta_{\mathcal{S},i} \mathrel{\vcenter{:}}= \sup_{v _{h,\tau}\in\mathbb{V}_{\mathcal{T},\mathcal{J}}\setminus\{0\}} \frac{\mathcal{S}_i(u _{h,\tau};m _{h,\tau};v _{h,\tau})}{\left\lVert v _{h,\tau}\right\rVert_X} \quad \forall i \in \{ 1,2 \}.\label{eq:etaStab-defs} \end{align}\tag{50}\] Since the stabilizations are linear in the third argument, the stabilization estimators are computable in practice by solving a discrete linear elliptic problem in the finite element space \(V_h\) on each time-interval. We refer the reader also to [23] for further discussion of how the computations can be made more efficient in practice by approximating the stabilization estimators via standard preconditioners for elliptic problems. However, for the general class of stabilizations above, it does not appear to be possible to localize the stabilization estimators across the spatial mesh in general. This motivates the further analysis in Section 6, where it will be seen that the computation of the stabilization estimators can be avoided in practice if the stabilizations have some additional structure.

5.1.0.4 Initial and final time condition estimators

Recall that the FEM scheme 41 imposes the initial condition \(m _{h,\tau}(0)=\Pi_0 m_0\) and final time condition \(u _{h,\tau}(T)=\Pi_T S[m _{h,\tau}(T)]\) for some user-chosen operators \(\Pi_0\) and \(\Pi_T\). The initial- and final-time condition estimators are defined by \[\eta_0 \mathrel{\vcenter{:}}= \left\lVert m_0 -\Pi_0 m_0\right\rVert_\Omega, \quad \eta_T \mathrel{\vcenter{:}}= \left\lVert S[m _{h,\tau}(T)]-\Pi_TS[m _{h,\tau}(T)]\right\rVert_\Omega. \label{eq:eta0-T-defs}\tag{51}\] One could alternatively consider the terms \(\eta_0\) and \(\eta_T\) as data oscillation terms. There is no practical difference in the different point of views, yet we choose to consider these terms as estimators since one can show their efficiency properties, see Remark 3 and Theorem 4 below.

5.1.0.5 Data oscillation terms

Finally, we define the temporal and spatial residual estimators that will appear below in Theorems 2 and Theorems 3 and 4 respectively. We note that it is not generally possible to avoid the dependence of the data oscillation on the numerical solution, owing to the nonlinearity of the problem. To define the temporal data oscillation terms, let \(\Pi_\tau\) denote the temporal \(L^2\)-orthogonal projection onto piecewise constant functions with respect to the time partition. In other words, for each \(v\in L^2(0,T)\), let \(\Pi_\tau v\in L^2(0,T)\) is defined by \(\Pi_\tau v|_{I_n}\mathrel{\vcenter{:}}= \frac{1}{\tau_n}\int_{I_n} v\mathrm{d}t\) for each time-interval \(I_n\), \(n\in\{1,\dots,N\}\). In a slight abuse of notation, we extend \(\Pi_\tau\) without change of notation to functions in general Bochner spaces \(L^2(0,T;W)\) where \(W\) is any Banach space. We now define the temporal data oscillation terms by \[\label{eq:temporal95oscillations} \begin{align} [\mathrm{osc}_{\tau,1}]^2 &\mathrel{\vcenter{:}}= \int_0^T \left\lVert (\mathrm{I}-\Pi_\tau)\left(F[m _{h,\tau}]-H[\nabla u _{h,\tau}]\right)\right\rVert_{H^{-1}(\Omega)}^2\mathrm{d}t, \\ [\mathrm{osc}_{\tau,2}]^2& \mathrel{\vcenter{:}}= \int_0^T \left\lVert(\mathrm{I}-\Pi_\tau)G\right\rVert_{H^{-1}(\Omega)}^2 + \left\lVert m _{h,\tau}(\mathrm{I}-\Pi_\tau)H_p[\nabla u _{h,\tau}] \right\rVert_{\Omega}^2 \mathrm{d}t, \end{align}\tag{52}\] where \(\mathrm{I}\) denotes the identity operator. To define the spatial data oscillation terms, we consider a fixed but arbitrary polynomial degree \(\kappa \geq 0\) and we define the following spatial approximation operators. For each \(K \in \mathcal{T}\) and \(F \in \mathcal{F}^I\), let \(\Pi_{K,\kappa}:L^2(K) \to \mathcal{P}_\kappa(K)\) and \(\Pi_{F,\kappa}:L^2(F) \to \mathcal{P}_\kappa(F)\) denote the \(L^2\)-projections onto polynomials of degree \(\kappa \geq 0\) over \(K\) and \(F\) respectively. For each element \(K \in \mathcal{T}\) and time interval \(n \in \{1,...,N\}\) we define the local data oscillation terms by \[\begin{gather} [\mathrm{osc}_{K,n,\kappa,i}]^2 \mathrel{\vcenter{:}}= \sum_{K' \in \widetilde{\mathcal{T}}_K} \int_{I_n} \Big[ h_{K'}^2 \left\lVert r_{K',n,i}-\Pi_{K',\kappa} r_{K',n,i}\right\rVert_{K'}^2 \\ + \sum_{F \in \mathcal{F}_{K'}^I} h_F \left\lVert j_{F,n,i}-\Pi_{F,\kappa} j_{F,n,i}\right\rVert_F^2 \Big]\mathrm{d}t\quad i \in \{1,2\}. \end{gather}\] For each \(i \in \{1,2\}\), we define the global spatial oscillation terms by \([\mathrm{osc}_{\kappa,i}]^2\mathrel{\vcenter{:}}= \sum_{n=1}^N \sum_{K \in \mathcal{T}} [\mathrm{osc}_{K,n,\kappa,i}]^2\).

5.2 A posteriori error bounds↩︎

The first main result, given in Theorem 2 below, shows the reliability of the estimator defined in 42 , i.e.the estimator bounds from above the extended norm of the error, as defined in 40 , up to a multiplicative constant and temporal data oscillation terms.

The proof of Theorem 2 is deferred to Section 5.3 below.

The following two main results show the efficiency of the estimator, i.e.the estimator is bounded from above by the error, up to spatial data oscillation terms. Theorem 3 below shows the local efficiency of the temporal jump estimator and of the residual estimators.

Theorem 3 shows that the temporal jump and residual estimators are locally efficient over each time interval-element block \(I_n \times K\). The proof of Theorem 3 is given in Section 5.4.

Note also that it is not possible to show local efficiency or computability of the stabilization estimator for the general abstract class of stabilizations considered above, since it is not locally computable. However, we can show the global efficiency of all components of the estimator, as seen in Theorem 4 below.

The consequence of Theorem 4 is that the total error estimator is bounded from above by the approximation error, proving that the estimator does not overestimate the error, up to a multiplicative constant. Theorems 2 and 4 together establish the equivalence between the error and the estimator, up to data oscillation.

5.3 Proof of reliability↩︎

The starting point for the analysis is the equivalence result of Theorem 1, which implies that \(\left\lVert u-u _{h,\tau}\right\rVert_X + \left\lVert m-M _{h,\tau}\right\rVert_Y \lesssim \mathcal{R}(u _{h,\tau},M _{h,\tau})\), where it is recalled that \(\mathcal{R}(u _{h,\tau},M _{h,\tau})\) defined by 25 . In order to connect the residuals at the continuous level with the numerical scheme, it is helpful to write the definition of the FEM 41 in the equivalent form \[\label{eq:FEM-scheme-def-alt} \begin{align} \mathfrak{R}_i(v _{h,\tau}) = \mathcal{S}_i(u _{h,\tau},m _{h,\tau};v _{h,\tau}) \quad\forall v _{h,\tau}\in \mathbb{V}_{\mathcal{T},\mathcal{J}}, \; \forall i\in \{1,2\}, \end{align}\tag{53}\] where the discrete residual functionals \(\mathfrak{R}_i\in X^*\), \(i\in \{1,2\}\), are defined by \[\begin{align} \mathfrak{R}_1(v) &\mathrel{\vcenter{:}}= \int_0^T \left\langle F[m _{h,\tau}],v\right\rangle_{} \mathrm{d}t\tag{54} \\ &- \int_0^T \big[-(\partial_t U _{h,\tau}, v)_\Omega + (\nu \nabla u _{h,\tau},\nabla v)_\Omega + (H[\nabla u _{h,\tau}],v)_\Omega \big] \mathrm{d}t, \notag \\ \mathfrak{R}_2(v)&\mathrel{\vcenter{:}}= \int_0^T \left\langle G,v\right\rangle_{} \mathrm{d}t\tag{55} \\ &-\int_0^T \big[ (\partial_t M _{h,\tau},v)_{\Omega} + \left(\nu \nabla m _{h,\tau},\nabla v\right)_{\Omega} + \left(m _{h,\tau}H_p[\nabla u _{h,\tau}],\nabla v\right)_{\Omega} \big]\mathrm{d}t, \notag \end{align}\] for all \(v\times X\). By introducing the above discrete residual functionals, we can now give a suitable decomposition of the residuals \(R_1^X(u _{h,\tau},M _{h,\tau})\) and \(R_1^Y(u _{h,\tau},M _{h,\tau})\) that appear in \(\mathcal{R}(u _{h,\tau},M _{h,\tau})\), as shown in the following Lemma.

The decompositions of the residuals of Lemma 4 allows us to bound the dual norms \(R_1^X(u _{h,\tau},M _{h,\tau})\) and \(R_2^Y(u _{h,\tau},M _{h,\tau})\). First, we show how the discrete residual functionals \(\mathfrak{R}_i\) that appear as the first terms in the decompositions lead to the residual and stabilization estimators, plus the temporal data oscillation.

We now complete the proof of Theorem 2.

5.4 Proof of efficiency↩︎

We begin by proving local efficiency of the temporal jump and residual estimators.

5.4.1 Proof of Theorem 3↩︎

5.4.2 Proof of Theorem 4↩︎

6 Stabilization-free error bounds↩︎

We conclude the analysis by showing that under additional assumptions on the FEM scheme in 41 , the a posteriori error bound does not require explicit computation of the stabilization estimators. To demonstrate that this is attainable, we consider a general class of stabilization schemes defined by \[\label{eq:Stabfree-defs} \begin{align} \mathcal{S}_1(u _{h,\tau};m _{h,\tau};v _{h,\tau}) &\mathrel{\vcenter{:}}= \int_0^T \left(\partial_t U _{h,\tau},v _{h,\tau}\right)_{\Omega} - \left(\partial_t U _{h,\tau},v _{h,\tau}\right)_{\Omega,h} + \mathcal{S}_{h,1}(u _{h,\tau}; m _{h,\tau}; v _{h,\tau})\mathrm{d}t, \\ \mathcal{S}_2(u _{h,\tau};m _{h,\tau};w _{h,\tau}) &\mathrel{\vcenter{:}}= \int_0^T \left(\partial_t M _{h,\tau},w _{h,\tau}\right)_{\Omega,h} - \left(\partial_t M _{h,\tau},w _{h,\tau}\right)_{\Omega}+\mathcal{S}_{h,2}(u _{h,\tau}; m _{h,\tau}; v _{h,\tau})\mathrm{d}t, \end{align}\tag{69}\] where \(\mathcal{S}_{h,i}:[V_h]^3 \to \mathbb{R}\) denotes a spatial stabilization that is allowed to be nonlinear in its first two arguments, but is linear with respect to its third argument.

The main result of this section, presented in Theorem 7 below, is a stabilization-free a posteriori error bound which is independent of the stabilization estimator and efficient, up to data oscillation terms. The analysis requires conditions on the discretization scheme and stabilization terms. Including the nonnegativity Hypothesis [H:stabilization95main], we require four additional conditions: the time-discretization is larger than the squared mesh size (up to a hidden constant) [H:ht-ratio], the spatial stabilization is affine-preserving [H:affine-preserving] and Lipschitz continuous [H:Lipschitz], and the density approximation is uniformly bounded [H:Mht-infty].

6.1 Mass lumping error bound↩︎

In this section, we show that the component of the stabilization estimator that results from lumping the mass matrix is bounded by the temporal jump estimators. A key ingredient of the analysis is an original explicit formula for the mass lumping error given in Lemma 7 below, which is of independent interest. Recall that, for each edge \(E \in \mathcal{E}\) of the mesh, the vector \(t_E\) denotes a chosen tangent vector to the edge; the orientation of \(t_E\) is of no consequence in the following.

We now prove that the mass lumping error, when applied to the time reconstruction of a function, is bounded up to constants by the temporal jump estimator, under the condition that each time-step size is bounded from below by a constant times the square of the maximum mesh size:

1. There exists a constant \(C_{\mathcal{T}\mathcal{J}} > 0\) such that the discretization parameters satisfy \(\tau_n \geq C_{\mathcal{T}\mathcal{J}} h_K^2\) for all \(n \in \{1,...,N\}\) and \(K \in \mathcal{T}\).

Hypothesis [H:ht-ratio] is usually satisfied in practical computations, and it permits locally refined spatial meshes and long time-steps. This hypothesis is not to be confused with the much more restrictive assumption of a parabolic CFL condition.

6.2 Spatial stabilization↩︎

The analysis of the remaining spatial stabilization makes use of our previous results on affine-preserving spatial stabilization for steady-state MFG [23]. Here we denote the patch of elements containing a given vertex \(z \in \mathcal{V}^I\) by \(\omega_z\). This analysis requires three additional conditions, which are patchwise affine preservation, Lipschitz continuity, and uniform boundedness of the discrete densities.

1. The stabilizations are patchwise affine-preserving: for every interior vertex \(z \in \mathcal{V}^I\) and each \(i \in \{1,2\}\), it holds that \(\mathcal{S}_{h,i}(v_h; \widetilde{v}_h;\varphi_z) = 0\) whenever \(v_h\) and \(\widetilde{v}_h\) are affine functions over \(\omega_z\), i.e.\(v_h|_{\omega_z} \in \mathcal{P}_1(\omega_z)\) and \(\widetilde{v}_h|_{\omega_z} \in \mathcal{P}_1(\omega_z)\).

2. The spatial stabilization scheme \(\mathcal{S}_{h,i}(\cdot,\cdot)\) satisfies for each \(i \in \{1,2\}\) \[\begin{align} \left\lvert\mathcal{S}_{h,i}(v_h;w_h;\varphi_z) - \mathcal{S}_{h,i}(\widetilde{v}_h; \widetilde{w}_h;\varphi_z) \right\rvert& \\ \leq L_{\mathcal{S}} \left\lvert\omega_z \right\rvert_d^{\frac{1}{2}} &\left( \left\lVert\nabla (v_h - \widetilde{v}_h)\right\rVert_{\omega_z} + \left\lVert\nabla (w_h - \widetilde{w}_h)\right\rVert_{\omega_z} \right), \end{align}\] for all \(v_h, w_h, \widetilde{v}_h, \widetilde{w}_h \in V_h\) and \(z \in \mathcal{V}^I\), for some bounded constant \(L_{\mathcal{S}} > 0\) that is independent of the discretization parameters.

3. The density approximation \(m _{h,\tau}\in \mathbb{V}_2\) satisfies \[\left\lVert m _{h,\tau}\right\rVert_{L^\infty(Q_T)} \leq \widetilde{M}_\infty,\] for some bounded constant \(\widetilde{M}_\infty > 0\) that is independent of the discretization parameters.

We refer the reader to [23] for some examples of stabilizations that verify the hypotheses above. Under these conditions, we can bound the spatial component of the stabilization scheme by the jump components of the residual estimators.

6.3 A posteriori bounds↩︎

Our last main result in Theorem 7 below shows that, under the above hypotheses, we can obtain a locally computable estimator in which the stabilization estimators can be removed.

Note that the local efficiency of the residual and temporal jump estimators was already shown above in Theorem 3, and the local efficiency of the initial and final time estimators was shown in Remark 3. Therefore, the estimator on the right-hand side of ?? is both locally computable and locally efficient.

Conclusion↩︎

We obtained the first a posteriori error bounds for a general class of stabilized finite element approximations of time-dependent mean field games. The estimators are reliable and efficient for a very general class of stabilizations. Under some additional structural assumptions on the stabilization and a weak condition on the discretization parameters, we also showed that the stabilization estimators are bounded by the temporal and spatial jump estimators, resulting in locally computable and locally efficient estimators.

References↩︎