Numerical Approximation of Young Measure Solutions to Parabolic Systems of Forward-Backward Type

TL;DR

This paper proves the existence and develops a numerical method for approximating Young measure solutions to complex nonlinear parabolic systems of forward-backward type, which are relevant in atmospheric modeling.

Contribution

It introduces a novel numerical algorithm for these systems and proves its convergence to Young measure solutions, addressing nonstandard growth and non-monotonicity.

Findings

Algorithm converges to Young measure solutions

Handles systems with nonstandard growth rates

Addresses problems beyond monotone operator theory

Abstract

This paper is concerned with the proof of existence and numerical approximation of large-data global-in-time Young measure solutions to initial-boundary-value problems for multidimensional nonlinear parabolic systems of forward-backward type of the form $\partial_t u - \mbox{div}(a(Du)) + Bu = F$, where $B \in \mathbb{R}^{m \times m}$, $Bv \cdot v \geq 0$ for all $v \in \mathbb{R}^m$, $F$ is an $m$-component vector-function defined on a bounded open Lipschitz domain $\Omega \subset \mathbb{R}^n$, and $a$ is a locally Lipschitz mapping of the form $a(A)=K(A)A$, where $K\,:\, \mathbb{R}^{m \times n} \rightarrow \mathbb{R}$. The function $a$ may have a nonstandard growth rate, in the sense that it is permitted to have unequal lower and upper growth rates. Furthermore, $a$ is not assumed to be monotone, nor is it assumed to be the gradient of a potential. Problems of this type arise in mathematical models of the atmospheric boundary layer and fall beyond the scope of monotone operator theory. We develop a numerical algorithm for the approximate solution of problems in this class, and we prove the convergence of the algorithm to a Young measure solution of the system under consideration.

Equations

$${\partial }_{t}u-\mathrm{div}\left(a(Du)\right)+Bu=F,$$

$$p_{i}>\max\bigg{\{}1,\frac{2n}{n+2}\bigg{\}}.$$

$$p:=\underset{i=1,\dots ,m}{\min }{p}_i\text{mbox}andq:=\underset{i=1,\dots ,m}{\max }{p}_i,$$

$$q-p<1.$$

$$F\in L^{p^{\prime }}(Q_T;{\mathbb{R}}^m)\cap L^2(Q_T;{\mathbb{R}}^m),$$

$$Bv\cdot v\ge 0\forall v\in {\mathbb{R}}^m,$$

$$u(0,\cdot )=u_0(\cdot )\in L^2(\Omega ;{\mathbb{R}}^m),$$

$$u{|}_{(0,T]\times \partial \Omega }=0.$$

$$a(A)=K(A)A,$$

$$c_0\sum _{i=1}^m({\mu }_i^2+|A_i{|}^2{)}^{\frac{p_i-2}{2}}\le K(A)\le c_1\sum _{i=1}^m({\mu }_i^2+|A_i{|}^2{)}^{\frac{p_i-2}{2}}\forall A\in {\mathbb{R}}^{m\times n},$$

$${\mu }_i\{\begin{array}{l}\in \mathbb{R}\text{if}{p}_i\ge 2,\\ \ne 0\text{if}1<p_i<2,\end{array}$$

$$\displaystyle K(\partial_{x}u)=\left\{\begin{array}[]{cc}\frac{(|\partial_{x}u_{1}|^{2}+|\partial_{x}u_{2}|^{2})^{\frac{3}{2}}}{|\partial_{x}u_{1}|^{2}+|\partial_{x}u_{2}|^{2}+\partial_{x}u_{3}}&\text{if}\ \partial_{x}u_{3}>0,\\ {}\hfil&{}\hfil\\ \sqrt{|\partial_{x}u_{1}|^{2}+|\partial_{x}u_{2}|^{2}-\partial_{x}u_{3}}&\text{if}\ \partial_{x}u_{3}\leq 0.\end{array}\right.$$

$${\partial }_{t}u-{\partial }_x(K({\partial }_{x}u){\partial }_{x}u)+Bu$$

$$u(t,1)$$

$$u_1(t,0)$$

$$u_2(t,0)$$

$${\partial }_x{u}_3(t,0)$$

$$u(0,x)=u_0(x)\text{for}x\in (0,1).$$

$$B=\left(\begin{array}{ccc}0& 1& 0\\ -1& 0& 0\\ 0& 0& 0\end{array}\right)$$

$$F(t,x)=\left(\begin{array}{c}-V\\ U\\ 0\end{array}\right).$$

$$c_0(|{\xi }_1|+|{\xi }_2|+\sqrt{|{\xi }_3|})\le K(\xi )\le c_1(|{\xi }_1|+|{\xi }_2|+\sqrt{|{\xi }_3|})$$

$$u\in L^{\infty }(0,T;L^2(\Omega ;{\mathbb{R}}^m))\cap L^p(0,T;W_0^{1,p}(\Omega ;{\mathbb{R}}^m))\text{mbox}withD{u}_i\in L^{p_i}(0,T;L^{p_i}(\Omega ;{\mathbb{R}}^{m\times n})),i=1,\dots ,m,$$

$${\partial }_{t}u\in L^{{\hat{q}}^{\prime }}(0,T;W^{-1,{\hat{q}}^{\prime }}(\Omega ;{\mathbb{R}}^m)),$$

$$⟨\nu ,a⟩\in L^{{\hat{q}}^{\prime }}(Q_T;{\mathbb{R}}^{m\times n}),$$

$${\int }_0^T⟨{\partial }_{t}u,\varphi {⟩}_{(W_0^{1,\hat{q}}(\Omega ;{\mathbb{R}}^m),W^{-1,{\hat{q}}^{\prime }}(\Omega ;{\mathbb{R}}^m))}+{\int }_{\Omega }⟨{\nu }_{t,x},a⟩:D\varphi +Bu\cdot \varphi \mathrm{d}x\mathrm{d}t={\int }_0^T{\int }_{\Omega }F\cdot \varphi \mathrm{d}x\mathrm{d}t$$

$$a=\nabla \Phi$$

$$E_r:=\left\{g\in C({\mathbb{R}}^{m\times n};{\mathbb{R}}^m):\underset{|A|\to \infty }{\lim }\frac{g(A)}{1+|A{|}^r}\text{exists}\right\},$$

$$\Vert g{\Vert }_{E_r}=\underset{A\in {\mathbb{R}}^{m\times n}}{\sup }\frac{|g(A)|}{1+|A{|}^r}.$$

$${\tilde{E}}_r:=\left\{g\in C({\mathbb{R}}^{m\times n};{\mathbb{R}}^m):\underset{A\in {\mathbb{R}}^{m\times n}}{\sup }\frac{|g(A)|}{1+|A{|}^r}<\infty \right\},$$

$$E_r^{\prime }\subset \mathcal{M}({\mathbb{R}}^{m\times n};{\mathbb{R}}^m)=C_0({\mathbb{R}}^{m\times n};{\mathbb{R}}^m{)}^{\prime }.$$

$$\underset{j\to \infty }{\lim }|\{x\in \Omega :z_j(x)\notin U\}|=0$$

$$\underset{s\to \infty }{\lim }\underset{k}{\sup }|\{x\in \Omega \cap B_R(0):|z_{j_k}(x)|\ge s\}|=0,$$

Full text

Dedicated to Academician Professor Gradimir Milovanović on the occasion of his 70th birthday.

1. Introduction

The paper is concerned with the numerical approximation of Young measure solutions to initial-boundary-value problems for nonlinear multidimensional parabolic systems of forward-backward type, where the existence of a weak solution cannot in general be guaranteed, because the nonlinearity in the equation is neither monotone, nor globally Lipschitz, nor indeed is it the gradient of a potential. Nonlinear parabolic systems of this type arise in certain mathematical models of the atmospheric boundary layer and, to date, there have been no attempts at the rigorous mathematical analysis of numerical methods for their approximate solution.

The systems of nonlinear parabolic partial differential equations considered here are of the following form:

[TABLE]

with $u:[0,T]\times\overline{\Omega}\rightarrow\mathbb{R}^{m}$, where $\Omega\subset\mathbb{R}^{n}$ is a bounded open Lipschitz domain. Here $Du$ denotes the spatial gradient of $u$, with $(Du)_{ik}=\partial u_{i}/\partial x_{k}$, $i=1,\dots,m$, $k=1,\dots,n$. Throughout we will take real numbers $p_{i}$, for $i=1,\ldots,m$, that satisfy

[TABLE]

We define

[TABLE]

and require that

[TABLE]

We shall also assume that

[TABLE]

with $\frac{1}{p}+\frac{1}{p^{\prime}}=1$, where $Q_{T}:=(0,T)\times\Omega$, and that

[TABLE]

where $B\in\mathbb{R}^{m\times m}$ is a constant matrix. All of our results extend directly to the case where $B\in L^{\infty}(Q_{T};\mathbb{R}^{m\times m})$, without assuming (1.4), however a more general supposition of this kind would contribute no new insight, so for the sake of clarity of the exposition we shall continue to assume that $B\in\mathbb{R}^{m\times m}$ is a constant matrix satisfying (1.4). The system (1.1) is supplemented with the initial condition

[TABLE]

and the homogeneous Dirichlet boundary condition

[TABLE]

The function $a:\mathbb{R}^{m\times n}\rightarrow\mathbb{R}^{m\times n}$ is assumed to be a locally Lipschitz mapping of the form

[TABLE]

where $K:\mathbb{R}^{m\times n}\rightarrow\mathbb{R}$ is continuous, and there exist positive constants $c_{0}$ and $c_{1}$, with $c_{0}\leq c_{1}$, such that

[TABLE]

with $p_{i}>\max\{1,\frac{2n}{n+2}\}$ for all $i=1,\ldots,m$, and $\mu_{i}$ is a constant defined by

[TABLE]

and where $|A_{i}|$ denotes the Euclidean norm of the $i$-th row of the matrix $A$. More generally, $|\cdot|$ will signify the absolute value of a real number, the Euclidean norm of a vector, or the Frobenius norm of a matrix: its precise meaning will be clear from the context. We emphasize that we are making no assumptions here about the monotonicity of the mapping $a$, nor will we assume that $a$ is the gradient of a certain $C^{1}$-function. The results of the paper can be easily adapted to systems in which the form of $a$ is slightly more general: for example, instead of (1.7)–(1.9) one can assume that $|a(x)|\leq c_{0}(1+|x|^{q-1})$ and $a(x)\cdot x\geq c_{1}|x|^{p}$, with $p>\max\{1,\frac{2n}{n+2}\}$ and $p\leq q<p+1$.

In order to motivate the weak structural assumptions on the function $a$ stated above we shall now present some examples.

Example 1.1**.**

In [2] the authors consider a system of equations modelling the behaviour of the atmospheric boundary layer, which served as our original motivation for the present study. In that model, $n=1$, $m=3$, $\Omega=(0,1)$, and the function $K$, referred to as a stability function, is defined as follows:

[TABLE]

The resulting system of PDEs for the vector of unknowns $u(t,x)=(u_{1}(t,x),u_{2}(t,x),u_{3}(t,x))^{\rm T}$ is

[TABLE]

subject to the initial condition

[TABLE]

Here, $u:[0,T]\times[0,1]\rightarrow\mathbb{R}^{3}$ is a vector with components $u=(u_{1},u_{2},u_{3})^{\rm T}$, $B$ is the skew-symmetric matrix

[TABLE]

and $F$ is a forcing term given by the constant vector

[TABLE]

What is immediately seen is that, given the above definition of the function $K$, the product $K(\partial_{x}u)\partial_{x}u$ lacks a coercivity estimate in a Sobolev space on the range $\partial_{x}u_{3}>0$, because the mapping $\xi\in\mathbb{R}^{3}\mapsto a(\xi):=K(\xi)\xi\in\mathbb{R}^{3}$ is not coercive in the range $\xi_{3}>0$; indeed, for $\xi=(1,0,\xi_{3})^{\rm T}$ with $\xi_{3}>0$, $a(\xi)\cdot\xi/|\xi|\rightarrow 1$ as $\xi_{3}\rightarrow+\infty$. Another difficulty arises from the fact that the function $\xi\mapsto a(\xi):=K(\xi)\xi$ is not monotone on $\mathbb{R}^{3}$, i.e., it is not true that $(a(\xi)-a(\eta)):(\xi-\eta)\geq 0$ for all $\xi,\eta\in\mathbb{R}^{3}$, and therefore monotone operator theory cannot be applied in this setting. To see that the mapping $\xi\mapsto a(\xi):=K(\xi)\xi$ is not monotone on $\mathbb{R}^{3}$, take $\xi=(0.035,0,-0.01)^{\rm T}$ and $\eta=(0.05,0,0)^{\rm T}$; then, computing $(a(\xi)-a(\eta))\cdot(\xi-\eta)$ will a give a negative number. Similarly taking $\xi=(-0.2,-0.1,0.2)^{\rm T}$ and $\eta=(-0.1,0,0.5)^{\rm T}$ will result in a negative number. Finally to show that the vector field $a$ is not a potential, consider the contour integral of $a$ around the circle $C$ parametrized by $r(\theta)=(1+\cos(\theta),0,\sin(\theta))$, $\theta\in[0,2\pi)$. One can show that this is nonzero, and so $a$ is not a conservative vector-field, and thus it is not a potential.

Example 1.2**.**

For $n=1$ and $m=3$ consider the function $K$ given by $\xi\in\mathbb{R}^{3}\mapsto K(\xi)=\sqrt{\xi_{1}^{2}+\xi_{2}^{2}+|\xi_{3}|}$. Then we are in the regime $\partial_{x}u_{3}\leq 0$ from (1.13). One has the existence of constants $c_{0}$ and $c_{1}$ such that

[TABLE]

through bounding and using the equivalence of norms on $\mathbb{R}^{3}$. This function $K$ is therefore covered by our assumptions. We are unfortunately unable to say anything about the region from (1.13) where $\partial_{x}u_{3}>0$, as our analysis does not cover the case $p=1$.

As will transpire from the discussion that follows, under the stated hypotheses on the function $a$ we are unable to prove the existence of a weak solution to the problem (1.1)–(1.9) under consideration and have to weaken the notion of solution to be able to show its existence. We shall therefore, instead, show the existence of a Young measure solution in the sense of the next definition, and will then consider the numerical approximation of such Young measure solutions.

Definition 1.3**.**

We say that a pair $(u,\nu)$, where $u$ is a function with

[TABLE]

and

[TABLE]

and $(t,x)\in Q_{T}\mapsto\nu_{t,x}$ is an $\mathbb{R}^{m\times n}$-valued Young measure, such that

[TABLE]

is called a Young measure solution to the problem described by (1.1), (1.5) and (1.6) if

[TABLE]

for all $\phi\in L^{\hat{q}}(0,T;W_{0}^{1,\hat{q}}(\Omega;\mathbb{R}^{m}))$, where $\hat{q}:=\max\{q,2\}$, $\hat{q}^{\prime}=\frac{\hat{q}}{\hat{q}-1}$, and the initial condition is satisfied in the sense that $u(0,x)=u_{0}(x)$ for almost every $x\in\Omega$.

Measure-valued solutions and Young measure solutions to parabolic partial differential equations have been studied by a number of authors. Frehse & Specovius-Neugebauer, in [12], showed the existence of Hölder continuous Young measure solutions to two-dimensional nonmonotone quasilinear problems exhibiting quadratic growth (that is, $p=q=n=2$), under the assumption that their nonlinearity $a$ was expressible as

[TABLE]

for some $C^{1}$ function $\Phi$. Our terminology Young measure solution follows that of Frehse & Specovius-Neugebauer in [12]; in particular, we have consciously avoided referring to the solutions considered here as measure-valued solutions, as the function $u$, whose existence we prove, is still a real-valued function with spatial Sobolev regularity; however, instead of being a standard weak solution, the function $u$ satisfies the PDE in the sense of gradient Young-measures. Young measure solutions to forward-backward parabolic problems were studied by Demoulini in [7]: under the assumptions that $m=1$, $p=q=2$ and that (1.14) is satisfied, a sequence of approximating solutions was constructed by means of minimizing an integral, which involves the convexification of $\Phi$. Using this construction, a family of Young measures was generated displaying an independence property, from which it was possible to deduce a uniqueness result regarding the function $u$ (not regarding the family of Young measures). From this, a weak-strong uniqueness result follows if the classical solution, which is assumed to exist, satisfies additional constraints on its gradient. More recently, in Thanh [20] and some of the references contained therein, in particular [21], possible regularizing terms are discussed for forward-backward parabolic equations, the inclusion of which allows for solutions with higher regularity to be obtained than those considered here; Thanh also obtained results concerning the long-time behaviour of certain Young measure solutions, and (in the case of an equation in one space dimension) the support of the Young measure appearing in the Young measure solutions. However, once again, that work makes the assumption that (1.14) holds, and assumes more restrictive growth conditions on $a$ than those considered here. Even more recently, Kim & Yan in [16] have expanded their earlier work in [15] to show the existence of infinitely many Lipschitz solutions (and under further assumptions, the existence of a unique classical solution) for particular classes of forward-backward parabolic equations. Although their work allows one to consider stronger notions of solution than the Young measure solutions that we shall discuss here, many of their structural assumptions regarding the nonlinearities are incompatible with our setting. Furthermore, their focus is solely on scalar equations, whereas we are concerned with systems here.

From the computational point of view, Carstensen & Roubíček [5] considered the numerical approximation of Young measures arising in certain minimization problems in the calculus of variations. As we lack this variational structure ourselves (in this work we are not assuming (1.14)), neither of the methods in [7] or [5] will be suitable in our setting. In [3] the authors consider numerical experiments for particular forward-backward equations that have a regularizing term added, but only for small values of the regularizing parameter and do not consider numerical schemes for the limiting case. We mention in passing that for the Keller–Segel system, which is a coupled parabolic-elliptic system of partial differential equations, a stochastic interacting particle approximation of global-in-time measure-valued solutions in two space dimensions was considered by Haškovec & Schmeiser in [14].

The aim of this work is to develop and analyze a numerical scheme for the approximation of Young measure solutions of parabolic systems of equations under minimal assumptions on the form taken by the nonlinearity $a$, and under no restriction on the dimension of the domain $\Omega$. We make the point here that the scheme discussed in this article will not allow us to compute the Young measures appearing in the references [7, 12, 5, 20, 3] cited above as taking different approximating sequences to the same function will potentially result in different families of Young measures being generated. Furthermore, without a concept of uniqueness, there is no guarantee that the functions appearing in the Young measure solutions will even be the same.

The paper is structured as follows. In the next section we shall collect miscellaneous results that are used in the subsequent analysis. In Section 3 we will show the existence of a large-data global-in-time Young measure solution to the system (1.1)–(1.9) under consideration by performing a spatial finite element approximation and passing to the limit in the spatial discretization parameter. In Section 4, we shall also discretize with respect to the temporal variable and show the convergence of the fully discrete scheme to a Young measure solution of the problem. We note here, however, that the results from Section 3 and Section 4 only yield to us indirect information about the measure $\nu$ through its action on functions belonging to a particular function space. This is the motivation behind Section 5, in which, stimulated by recent contributions by Fjordholm et al. [11] and [10], in particular, we shall develop a numerical algorithm for the computation of Young measure solutions which allows for a more direct approximation of the Young measure $\nu$, and will prove that the algorithm converges to a Young measure solution of the problem under consideration.

2. Preliminary Definitions and Results

2.1. Function Spaces

We begin by defining, for $r>0$, the function space

[TABLE]

which will be of use to us later. On this space we consider the norm

[TABLE]

Lemma 2.1**.**

The following statements hold true for $r>0$:

1. (1)

The space $E_{r}$ is separable; 2. (2)

The space $E_{r}$ is metrizable.

Proof.

For the first property we refer to Kinderlehrer & Pedregal [17]. The second follows by taking the metric induced by the norm $\|\cdot\|_{E_{r}}$. ∎

Remark 2.2**.**

We remark here that $E_{r}$ is a subspace of the larger space

[TABLE]

which is not separable. For example, under the assumption (1.8), $a\in\tilde{E}_{q-1}$, and for any $\alpha\in(1,\frac{p}{q-1})$ we have that $0<q-1<\frac{p}{\alpha}<p$, and therefore $a\in E_{\frac{p}{\alpha}}$.

As $C_{0}(\mathbb{R}^{m\times n};\mathbb{R}^{m})\subset E_{r}$ for all $r>0$, we have the following nesting of dual spaces for all $r>0$:

[TABLE]

Since the dual of a normed space is a Banach space, we have that $E_{r}^{\prime}$ is a Banach space for all $r>0$, as we can equip $E_{r}$ with the norm induced by its metric. In particular, thanks to the Banach–Alaoglu Theorem, any bounded sequence in $E_{r}^{\prime}$ has a weakly-$\ast$ convergent subsequence.

2.2. Young Measures

We recall here from [1] the so-called Fundamental Theorem for Young Measures, together with some related remarks.

Theorem 2.3**.**

Let $\Omega\subset\mathbb{R}^{n}$ be Lebesgue measurable, let $\mathfrak{K}\subset\mathbb{R}^{m}$ be closed and let $z_{j}:\Omega\rightarrow\mathbb{R}^{m}$, $j=1,2,\ldots,$ be a sequence of Lebesgue measurable functions satisfying

[TABLE]

for any open neighbourhood $U$ of $\mathfrak{K}$ in $\mathbb{R}^{m}$. Then, there exists a subsequence $\{z_{j_{k}}\}$ of $\{z_{j}\}$ and a family of positive measures on $\mathbb{R}^{m}$, $\{\nu_{x}\}_{x\in\Omega}$, depending measurably on $x\in\Omega$, such that:

1. (1)

$\|\nu_{x}\|_{M}:=\int_{\mathbb{R}^{m}}{\rm d}\nu_{x}\leq 1$* for almost every $x\in\Omega$;* 2. (2)

$\mathrm{supp}\ \nu_{x}\subset\mathfrak{K}$* for almost every $x\in\Omega$; and* 3. (3)

$f(z_{j_{k}})\overset{\ast}{\rightharpoonup}\langle\nu_{x},f\rangle=\int_{\mathbb{R}^{m}}f(\lambda){\,{\rm d}}\nu_{x}(\lambda)$* in $L^{\infty}(\Omega)$ for each continuous function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ satisfying $\lim_{|\lambda|\rightarrow\infty}f(\lambda)=0$.*

Suppose further that $\{z_{j_{k}}\}$ satisfies the following tightness condition:

[TABLE]

for every $R>0$. Then, we have that $\|\nu_{x}\|_{M}=1$ for almost every $x\in\Omega$, and, for any measurable subset $A$ of the set $\Omega$,

[TABLE]

in $L^{1}(A)$ for any continuous function $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ such that $\{f(z_{j_{k}})\}$ is sequentially weakly relatively compact in $L^{1}(A)$.

Remark 2.4**.**

If $\Omega$ is bounded and the sequence $\{z_{j}\}$ is bounded in $L^{p}(\Omega;\mathbb{R}^{m})$ for some $p\in(1,\infty)$, then we obtain from the theorem the existence of a family of probability measures $\{\nu_{x}\}_{x\in\Omega}$ and a subsequence $z_{j_{k}}$ such that

[TABLE]

in $L^{r}(\Omega)$, where $f:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is a continuous function satisfying

[TABLE]

with $1<r<\frac{p}{s}$ and $s>0$.

Remark 2.5**.**

The assumption that

[TABLE]

for any open neighbourhood $U$ of $\mathfrak{K}$ is only used to prove that (2) in Theorem 2.3 holds, and is not needed to obtain any of its other conclusions. The proof of this theorem from [1] is still valid if we take $\mathfrak{K}=\mathbb{R}^{m}$, only it then provides no information about the support of the Young measure $\nu$. Therefore, if we are not concerned with the support of the Young measure $\nu$ we may always take $\mathfrak{K}=\mathbb{R}^{m}$ in order to apply this theorem.

2.3. Miscellaneous Results

We recall the following result from Strauss [19] (cf. also Lions & Magenes [18], Lemma 8.1, Ch. 3, Sec. 8.4).

Lemma 2.6**.**

Suppose that $X$ and $Y$ are Banach spaces. Assume that the space $X$ is reflexive and is continuously embedded in the space $Y$​; then,

[TABLE]

where $C_{w}([0,T];X)$ and $C_{w}([0,T];Y)$ denote the spaces of weakly continuous functions from $[0,T]$ into $X$ and $Y$, respectively.

We shall also require the following consequence of the Arzelà–Ascoli theorem.

Lemma 2.7**.**

Let $\{u^{j}\}$ be a sequence of functions with the following properties:

1. (1)

$\{u^{j}\}$* is bounded in $L^{\infty}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$ and*

[TABLE] 2. (2)

$\{\partial_{t}u^{j}\}$* is bounded in $L^{s^{\prime}}\!(0,T;W^{-1,s^{\prime}}\!(\Omega;\mathbb{R}^{m}))$ for some $s\in(1,\infty)$, and*

[TABLE]

Then, $u^{j_{k}},u\in C_{w}([0,T];L^{2}(\Omega;\mathbb{R}^{m})$, and there is a subsequence $\{u^{j_{k}}\}$ such that

[TABLE]

uniformly in $C([0,T])$ for all $w\in L^{2}(\Omega;\mathbb{R}^{m})$ as $k\rightarrow\infty$.

Proof.

We shall apply the Arzelà–Ascoli theorem for sequences of uniformly bounded and equicontinuous functions. For $w\in L^{2}(\Omega;\mathbb{R}^{m})\cap W_{0}^{1,s}(\Omega;\mathbb{R}^{m})$, we let $f^{j}:[0,T]\rightarrow\mathbb{R}$ be the family of functions defined by

[TABLE]

Since $\partial_{t}u^{j}\in L^{s^{\prime}}\!(0,T;W^{-1,s^{\prime}}\!(\Omega;\mathbb{R}^{m}))$, it automatically follows that $u^{j}\in C([0,T];W^{-1,s^{\prime}}\!(\Omega;\mathbb{R}^{m}))$, whereby also $u^{j}\in C_{w}([0,T];W^{-1,s^{\prime}}\!(\Omega;\mathbb{R}^{m}))$. As $u^{j}\in L^{\infty}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$, it then follows from Lemma 2.6 with $X=L^{2}(\Omega;\mathbb{R}^{m})$ and $Y=L^{2}(\Omega;\mathbb{R}^{m})+W^{-1,s^{\prime}}\!(\Omega;\mathbb{R}^{m})=(L^{2}(\Omega;\mathbb{R}^{m})\cap W^{1,s}_{0}(\Omega;\mathbb{R}^{m}))^{\prime}$, that $u^{j}\in C_{w}([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$. We note here that the duality $(L^{2}(\Omega;\mathbb{R}^{m})\cap W^{1,s}_{0}(\Omega;\mathbb{R}^{m}))^{\prime}=L^{2}(\Omega;\mathbb{R}^{m})+W^{-1,s^{\prime}}\!(\Omega;\mathbb{R}^{m})$ is a consequence of the Duality Theorem (cf. Theorem 2.7.1 on p.32 in Bergh & Löfström [4]), because $L^{2}(\Omega;\mathbb{R}^{m})\cap W^{1,s}_{0}(\Omega;\mathbb{R}^{m})$ is dense in both $L^{2}(\Omega;\mathbb{R}^{m})$ and $W^{1,s}_{0}(\Omega;\mathbb{R}^{m})$ (e.g., because $C^{\infty}_{0}(\Omega;\mathbb{R}^{m})$ is dense in both $L^{2}(\Omega;\mathbb{R}^{m})$ and $W^{1,s}_{0}(\Omega;\mathbb{R}^{m})$). Hence, $f^{j}\in C([0,T])$. We further have that

[TABLE]

which gives the boundedness of the sequence $\{f^{j}\}$ in $C([0,T])$ for each $w\in L^{2}(\Omega;\mathbb{R}^{m})\cap W_{0}^{1,s}(\Omega;\mathbb{R}^{m})$. In order to verify equicontinuity of the sequence, note that, for all $w\in L^{2}(\Omega;\mathbb{R}^{m})\cap W_{0}^{1,s}(\Omega;\mathbb{R}^{m})$,

[TABLE]

Thus by the Arzelà–Ascoli theorem there is a subsequence $\{f^{j_{k}}\}$, which converges uniformly as $k\rightarrow\infty$ to some function $f\in C([0,T])$. We note that as this convergence is uniform and the sequence $f^{j_{k}}$ is uniformly bounded we have that

[TABLE]

However, by the assumed weak convergence, we also have that

[TABLE]

allowing us to identify

[TABLE]

Thus we have shown that, for all $w\in L^{2}(\Omega;\mathbb{R}^{m})\cap W_{0}^{1,s}(\Omega;\mathbb{R}^{m})$,

[TABLE]

As $L^{2}(\Omega;\mathbb{R}^{m})\cap W_{0}^{1,s}(\Omega;\mathbb{R}^{m})$ is dense in $L^{2}(\Omega;\mathbb{R}^{m})$, it then follows that, for all $w\in L^{2}(\Omega;\mathbb{R}^{m})$,

[TABLE]

which in particular means that for all $t\in[0,T]$ we have $u^{j_{k}}(t,\cdot)\rightharpoonup u(t,\cdot)$ in the space $L^{2}(\Omega;\mathbb{R}^{m})$. ∎

Finally, we recall the following standard Lebesgue space interpolation result, whose proof is omitted.

Lemma 2.8**.**

Let $u:[0,T]\times\Omega\rightarrow\mathbb{R}^{m}$, where $\Omega\subset\mathbb{R}^{n}$ is open and bounded, satisfy

[TABLE]

for some $p\geq 1$. Then,

[TABLE]

3. Spatial discretization and convergence of the semidiscrete problem

In this section we prove the existence of solutions to the system (1.1)–(1.9). In preparation for the construction of the fully discrete numerical approximation of the class of problems under consideration described in the next section, our proof of existence of large-data global-in-time Young measure solutions will be based on performing a spatial finite element approximation. We shall therefore suppose henceforth that $\Omega$ is a Lipschitz polytope in $\mathbb{R}^{n}$. The proof presented below can be replicated on any bounded Lipschitz domain $\Omega$ by replacing the finite element basis, consisting of continuous piecewise linear basis functions satisfying a homogeneous Dirichlet boundary condition on $\partial\Omega$ by an abstract Galerkin basis. Suppose further that $0<h_{0}\ll\mbox{diam}(\Omega)$ and that $\{\mathcal{T}_{h}\}_{h\in(0,h_{0}]}$ is a shape-regular family of subdivisions $\mathcal{T}_{h}$ of $\overline{\Omega}$ into closed simplexes $\Delta$, where h=\max_{\Delta\in\mathcal{T}_{h}}\mbox{diam(\Delta)}. Consider the finite element space

[TABLE]

By $V^{h}_{m}$ we denote the space of $m$-component vector-valued functions, each of whose components lies in $V^{h}$. We shall assume that $\{\mathcal{T}_{h}\}_{h\in(0,h_{0}]}$ is such that the $L^{2}(\Omega;\mathbb{R}^{m})$ orthogonal projector $P^{h}$ onto $V^{h}_{m}$ is stable in $W^{1,\hat{q}}_{0}(\Omega;\mathbb{R}^{m})$; i.e., there exists a positive constant $C$, independent of $h\in(0,h_{0}]$, such that

[TABLE]

On globally quasiuniform subdivisions (3.1) is a consequence of a global inverse inequality and inequality (7) in [8]; we note, however, that the stability inequality (3.1) is in fact valid under less restrictive assumptions on the subdivision than quasiuniformity (see, for example, [6]).

Let $\{\phi^{h}_{i}\}_{i=1}^{\mathcal{N}(h)}$ be a basis for $V^{h}_{m}$. We shall therefore seek an approximate solution $u_{h}\in V^{h}_{m}$ in the form

[TABLE]

satisfying

[TABLE]

with $u^{h}_{0}\in V^{h}_{m}$ and $u^{h}_{0}\rightarrow u_{0}$ (strongly) in $L^{2}(\Omega;\mathbb{R}^{m})$ as $h\rightarrow 0_{+}$. The system given by (3.2a) is a system of ODEs for the coefficients $\alpha^{h}_{i}$, $i=1,\dots,\mathcal{N}(h)$, and a solution to (3.2a), (3.2b) exists locally by Peano’s Theorem on some interval $[0,T_{h})\subset[0,T]$, thanks to the assumed (local Lipschitz) continuity of $a$. Our existence theorem below is based on proving that $u^{h}$ can be extended to the final time $T>0$ for each $h\in(0,h_{0}]$, and that as $h\rightarrow 0_{+}$ a subsequence of the sequence of approximate solutions converges, in a sense to be made precise, to a Young measure solution of the problem.

Theorem 3.1**.**

Suppose that $\Omega\subset\mathbb{R}^{n}$ is a Lipschitz polytope, let $u_{0}\in L^{2}(\Omega;\mathbb{R}^{m})$, assume that

[TABLE]

satisfy $q-p<1$, let $F\in L^{p^{\prime}}\!(Q_{T};\mathbb{R}^{m})\cap L^{2}(Q_{T};\mathbb{R}^{m})$, and suppose that $a:\mathbb{R}^{m\times n}\rightarrow\mathbb{R}^{m\times n}$ is a locally Lipschitz mapping satisfying the assumptions (1.7)–(1.9). Then, there exists a Young measure solution $(u,\nu)$ of the problem (1.1) with data given by (1.5) and (1.6). Furthermore, there exists a subsequence (not indicated) of solutions $\{u^{h}\}$ to the semidiscrete problem (3.2a), (3.2b) such that

[TABLE]

[TABLE]

[TABLE]

Proof.

As was noted above, (3.2a), (3.2b) is an initial-value problem for a system of ODEs for the coefficients $\alpha^{h}_{i}$, and a solution exists locally by Peano’s Theorem on some interval $[0,T_{h})\subset[0,T]$, thanks to the assumed (local Lipschitz) continuity of $a$. We begin by showing that we can extend the numerical solution $u^{h}\in V^{h}_{m}$, defined on $[0,T_{h})\times\overline{\Omega}$ up to time $T$ for all $h\in(0,h_{0}]$, so that it is defined on the whole of $[0,T]\times\overline{\Omega}$. To this end, we take $v^{h}=u^{h}$ in (3.2a); hence,

[TABLE]

which, by using (1.4) and (1.8) on the left-hand side and Hölder’s inequality on the right-hand side, yields

[TABLE]

In order to bound the second term on the left-hand side of (3.4) from below, we note that

[TABLE]

and the following inequality holds for all $s\in[0,\infty)$ if $1<p_{i}<2$ and $\mu_{i}\neq 0$:

[TABLE]

Using these, we deduce that

[TABLE]

Thus, using Young’s inequality and Poincaré’s inequality in $W^{1,p}_{0}(\Omega;\mathbb{R}^{m})$ in order to absorb the final factor on the right-hand side into the second term on the left-hand side, we arrive at the following energy inequality:

[TABLE]

which provides uniform bounds on the sequence $\{u^{h}\}$, thanks to the strong convergence of $u_{0}^{h}\rightarrow u_{0}$ in $L^{2}(\Omega;\mathbb{R}^{m})$ as $h\rightarrow 0_{+}$, which implies in particular that $\|u_{0}^{h}\|_{L^{2}(\Omega;\mathbb{R}^{m})}$ is uniformly bounded in $h$. Hence, the right-hand side of (3.5) is bounded, independent of $h\in(0,h_{0}]$. This means that $T_{h}$ cannot be the maximal existence time, and $u^{h}$ can be therefore extended from $[0,T_{h})$ beyond $T_{h}$ to the whole of the time interval $[0,T]$ for all $h\in(0,h_{0}]$. Hence,

[TABLE]

and also

[TABLE]

The uniform bound (3.7) implies the existence of a function

[TABLE]

for which, up to a subsequence, we have

[TABLE]

and

[TABLE]

Furthermore, thanks to (1.8) and because the Euclidean norm of the vector $(|Du_{1}^{h}|,\ldots,|Du_{m}^{h}|)^{\rm T}$ is bounded by its 1-norm, we have the following bound:

[TABLE]

Now we seek to bound each of the terms appearing in the last sum in (3.10). Clearly for the case $i=k$ we have an expression of the form $(\mu_{i}^{2}+|Du_{i}^{h}|^{2})^{\frac{p_{i}-1}{2}}$ appearing, which is bounded in $L^{p_{i}^{\prime}}\!(Q_{T})$, uniformly with respect to $h\in(0,h_{0}]$, by the energy estimate (3.7). In the case $i\neq k$ we proceed by using Hölder’s inequality together with (3.6) and (3.7), our objective being to bound all of these terms in $L^{p_{i}^{\prime}}\!(Q_{T})$ if $p_{i}\geq 2$ or in $L^{2}(Q_{T})$ if $1<p_{i}<2$. We begin with the terms for which $p_{i}\geq 2$:

[TABLE]

Now for the terms in which $1<p_{i}<2$, thanks to (3.6), we have

[TABLE]

Thus we have bounded each term in the sum appearing on the right-hand side of (3.10) in $L^{p_{i}^{\prime}}\!(Q_{T})$ (if $p_{i}\geq 2$) or in $L^{2}(Q_{T})$ (if $1<p_{i}<2$), uniformly with respect to $h\in(0,h_{0}]$, and therefore $\{a(Du^{h})\}_{0<h\leq h_{0}}$ is a bounded sequence in $L^{{\hat{q}}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$. More precisely, there exists a positive constant $C$, independent of $h$, such that

[TABLE]

where we recall that $\hat{q}=\max\{q,2\}$, whereby $\hat{q}^{\prime}=\min\{q^{\prime},2\}$. Therefore, there is a subsequence (still indexed only by $h$) and a $\chi\in L^{{\hat{q}}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$ such that

[TABLE]

Using similar estimates to those above we can also show that there exists a positive constant $C$, independent of $h$, such that

[TABLE]

Indeed, for each $\varphi\in L^{\hat{q}}(0,T;W^{1,{\hat{q}}}_{0}(\Omega;\mathbb{R}^{m}))$, by (3.2a), (3.1) and (3.11), we have that

[TABLE]

and then (3.13) follows by noting that the dual space of $L^{\hat{q}}(0,T;W^{1,{\hat{q}}}_{0}(\Omega;\mathbb{R}^{m}))$ is $L^{\hat{q}^{\prime}}\!(0,T;W^{-1,{\hat{q}}^{\prime}}\!(\Omega;\mathbb{R}^{m}))$.

Now let $\psi\in L^{\hat{q}}(0,T;W_{0}^{1,{\hat{q}}}(\Omega;\mathbb{R}^{m}))$ be the strong limit, as $h\rightarrow 0_{+}$ and $l\rightarrow\infty$, of the sequence of functions

[TABLE]

in the $L^{\hat{q}}(0,T;W_{0}^{1,{\hat{q}}}(\Omega;\mathbb{R}^{m}))$ norm, with $\beta^{l}_{i}\in C([0,T])$ and $\phi^{h}_{i}\in V^{h}_{m}$, $i=1,\dots,\mathcal{N}(h)$, $l=1,2,\dots$; such a limit necessarily exists by density. We then have that

[TABLE]

The right-hand side of this inequality converges to zero in the limit as $h\rightarrow 0_{+}$ (and therefore $\mathcal{N}(h)\rightarrow\infty$) and $l\rightarrow\infty$. Thus we have that

[TABLE]

for all $\psi\in L^{\hat{q}}(0,T;W_{0}^{1,{\hat{q}}}(\Omega;\mathbb{R}^{m}))$.

The uniform in $h$ bounds established in (3.13) and (3.7) now allow passage to the limit in the first and third term on the left-hand side of (3.14): by the uniform in $h$ bound on $\|u^{h}\|_{L^{\infty}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))}$, there exists a subsequence (still indexed by $h$) which converges weakly-$\ast$ to $u$ in $L^{\infty}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$. This enables us to also pass to the limit in the third term of the integrand on the left-hand side of (3.14). In addition, thanks to (3.13), one can extract a further subsequence (which we continue to index by $h$) for which, by the uniqueness of the weak limit,

[TABLE]

meaning that

[TABLE]

for all $\psi\in L^{\hat{q}}(0,T;W^{1,{\hat{q}}}_{0}(\Omega;\mathbb{R}^{m}))$, where $\langle\cdot,\cdot\rangle$ is the duality pairing between $W^{-1,{\hat{q}}^{\prime}}\!(\Omega;\mathbb{R}^{m})$ and $W^{1,{\hat{q}}}_{0}(\Omega;\mathbb{R}^{m})$. Hence,

[TABLE]

for all $\psi\in L^{\hat{q}}(0,T;W^{1,{\hat{q}}}_{0}(\Omega;\mathbb{R}^{m}))$. The final convergence result in the statement of the theorem, that $u^{h}\rightarrow u$ in $C_{w}([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$, is a direct consequence of (3.8), (3.15) and Lemma 2.7.

In order to pass to the limit in the second term in (3.16), we shall invoke Theorem 2.3. To this end, we begin by verifying that the subsequence $\{Du^{h_{k}}\}$ satisfies the tightness condition from Theorem 2.3; i.e., that

[TABLE]

for every $R>0$. In order to show this we set

[TABLE]

and apply Chebyshev’s inequality:

[TABLE]

where we have used (3.7) along with the fact that for $r_{1}\geq r_{2}$, we have $|s|^{r_{2}}\leq 1+|s|^{r_{1}}$ to transition from the first line to the second. The resulting bound is independent of $k$ (and $h_{k}$) and $R$, and so by passing $s\rightarrow\infty$ we deduce that (3.17) holds.

By taking $\mathfrak{K}=\mathbb{R}^{m}$ (see Remark 2.5 following the statement of Theorem 2.3), we may apply Theorem 2.3, to deduce the existence of a parameterized family of Young measures $(\nu_{t,x})_{(t,x)\in Q_{T}}$ and a further subsequence (which we still index by $h$) for which

[TABLE]

where $1<\alpha<\frac{p}{q-1}$ (see Remark 2.4, combined with the fact that $a(A)\leq c(1+|A|^{q-1})$, which follows from (1.8) by recalling that $q:=\max_{i=1,\ldots,m}p_{i}$, whence $a\in E_{p/\alpha}$). By comparing (3.12) and (3.18) and noting the uniqueness of weak limits there is then a subsequence of $u^{h}$ (still indexed by $h$) such that

[TABLE]

Here we have used that $\hat{q}^{\prime}=\min\{q^{\prime},2\}=\min\{\frac{q}{q-1},2\}\geq\min\{\frac{p}{q-1},2\}\geq\min\{\alpha,2\}$, with $\alpha$ as above, and noted that (3.18) implies that $a(Du^{h})\rightharpoonup\langle\nu,a\rangle$ in $L^{\min\{\alpha,2\}}(Q_{T};\mathbb{R}^{m\times n})$ because $\alpha\geq{\min\{\alpha,2\}}$.

This allows us to complete the passage to the limit and obtain the existence of a pair $(u,\nu)$ such that

[TABLE]

with

[TABLE]

satisfying

[TABLE]

for all $\phi\in L^{\hat{q}}(0,T;W_{0}^{1,\hat{q}}(\Omega;\mathbb{R}^{m}))$.

It remains to be shown that the initial condition is satisfied. By Lemma 2.7 it follows that we can extract a further subsequence $h_{l}$ along which we have that

[TABLE]

uniformly in $C([0,T])$, for all $\varphi\in L^{2}(\Omega;\mathbb{R}^{m})$. Let us decompose

[TABLE]

The term $I_{l}$ converges to zero as $h_{l}\rightarrow 0_{+}$ by the assumed strong convergence in $L^{2}(\Omega;\mathbb{R}^{m})$ of the discretized initial condition $u_{0}^{h_{l}}$ to $u_{0}$, and the term $III_{l}(t)$ converges to zero, uniformly in $C([0,T])$ as $h_{l}\rightarrow 0_{+}$ by (3.20). Hence, for any $\varepsilon>0$, there exists an $h_{l}^{\ast}=h_{l}^{\ast}(\varepsilon,\varphi)\in(0,h_{0}]$, such that

[TABLE]

Consequently, for any $\varepsilon>0$, there exists an $h_{l}^{\ast}\in(0,h_{0}]$, such that

[TABLE]

With $h_{l}^{\ast}$ fixed, we now pass to the limit $t\rightarrow 0_{+}$ in this inequality. The term $II_{l}(t)$ converges to zero as $t\rightarrow 0_{+}$ because $u^{h_{l}^{\ast}}$, as a solution to the Galerkin equation (3.2a), is continuous and satisfies the discretized initial condition. The term $IV(t)$ converges to as $t\rightarrow 0_{+}$ since by Lemma 2.7 we have that $u\in C_{w}([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$, the space of weakly continuous functions from the interval $[0,T]$ into $L^{2}(\Omega;\mathbb{R}^{m})$. Thus, by passing to the limit $t\rightarrow 0_{+}$ in the last inequality we deduce that, for any $\varphi\in L^{2}(\Omega;\mathbb{R}^{m})$ and any $\varepsilon>0$,

[TABLE]

Taking $\varphi(x)=u_{0}(x)-u(0,x)$ and letting $\varepsilon\rightarrow 0_{+}$ then yields that $u_{0}(x)=u(0,x)$ for almost every $x\in\Omega$. ∎

4. The fully discrete scheme and its convergence analysis

In the course of the proof of Theorem 3.1 Peano’s Theorem was applied in order to show that one can solve the Galerkin equation (3.2a). We shall now employ an implicit Euler discretization scheme to numerically approximate the solution whose existence on $[0,T]$ is guaranteed by Peano’s Theorem. Starting from our discretized initial condition $u_{0}^{h}\in V^{h}_{m}$, we inductively define $u_{i+1}^{h}$, $i=0,\dots,N-1$, as solutions to the following fully discrete problem.

Let $\Delta t=T/N$ where $N\geq 1$, and given $u_{i}^{h}\in V^{h}_{m}$, find $u_{i+1}^{h}\in V^{h}_{m}$ such that

[TABLE]

Here,

[TABLE]

and, letting $t_{i}:=i\Delta t$, we define on $Q_{T}$ the function

[TABLE]

Given a function $u^{h}_{i}$, the fact that there exists a function $u^{h}_{i+1}$ satisfying (4.1) is standard, and can be shown using a simple consequence of Brouwer’s Fixed Point Theorem (cf., for example, Corollary 1.1 on p.279 of [13]). We begin the convergence analysis with the following estimate, which is reminiscent of the energy estimate obtained in the proof of Theorem 3.1, for the semidiscrete problem considered there.

Lemma 4.1**.**

There exists a positive constant $c$, independent of $h$ and $\Delta t$, such that the functions $u_{i}^{h}$, $i=1,\ldots,N$, satisfy, for all $h\in(0,h_{0}]$,

[TABLE]

We note here that $a(\xi):\xi\geq 0$ for all $\xi\in\mathbb{R}^{m\times n}$, and therefore each of the terms appearing on the left-hand side of the inequality (4.2) is nonnegative. In addition, the expression on the right-hand side of (4.2) is further bounded above by a constant that is independent of $h$ and $\Delta t$.

Proof.

We begin by choosing $\phi^{h}=u_{i+1}^{h}$ in (4.1), and noting the following identity:

[TABLE]

Using this in (4.1) gives

[TABLE]

Now we use that $Bv\cdot v\geq 0$ for all $v\in\mathbb{R}^{m}$, multiply by $2\Delta t$ and sum through $i=0,\ldots,M$ for some $M<N$ to get

[TABLE]

By omitting the second and the third term from the left-hand side of the inequality (4.3) noting the independence of its right-hand side of $M$, and then taking the maximum over $M=1,\dots,N-1$, yields the desired bound on the first term on the left-hand side of the inequality (4.2). We then return with that bound to the inequality (4.3) to further estimate its right-hand side from above, whilst omitting the first term from the left-hand side of (4.3); thus we arrive at the desired bound on the second and third term on the left-hand side of the inequality (4.2). The expression on the right-hand side of (4.2) is clearly independent of $\Delta t$. The final assertion in the statement of the lemma is a consequence of the assumed strong convergence of $u_{0}^{h}$ to $u_{0}$ in $L^{2}(\Omega;\mathbb{R}^{m})$. ∎

As the functions $u^{h}_{i}$, $i=0,\dots,N$, are defined only on $\{t_{0}=0,t_{1},\dots,t_{N}=T\}\times\overline{\Omega}$ rather than on the whole of $\overline{Q_{T}}$, we shall next extend them to $\overline{Q_{T}}$. Thus we define the following two functions:

[TABLE]

which is continuous and piecewise linear with respect to $t\in[0,T]$, and

[TABLE]

which is piecewise constant with respect to $t\in[-\Delta t,T]$. With these new notations, (4.1) can be written as

[TABLE]

for all $\phi^{h}\in V^{h}_{m}$. We aim now to obtain uniform bounds on the sequences $\{u_{\Delta t}^{h}\}$ and $\{\tilde{u}_{\Delta t}^{h}\}$, which will allow us to complete our passage to the limit in $\Delta t$, possibly along some subsequence.

Lemma 4.2**.**

There exists a positive constant $c$, independent of $h$ and $\Delta t$, such that, for all $h\in(0,h_{0}]$,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Furthermore, the right-hand sides of these inequalities are bounded above by a positive constant, independent of $h\in(0,h_{0}]$ and $\Delta t$.

Proof.

The first two bounds follow directly from the bounds on the first two terms in (4.2), while the fourth inequality is proved by an argument that is completely analogous to the proof of the bound on the second term appearing on the left-hand side of (3.7). We shall therefore focus our attention on showing the third bound in the statement of the lemma. To this end, we begin by noting that since $a(\xi):\xi\geq 0$ for all $\xi\in\mathbb{R}^{m\times n}$, we have that

[TABLE]

First suppose that $q>2$; then, $\hat{q}=q$ and $\hat{q}^{\prime}=q^{\prime}\in(1,2)$; hence, by Young’s inequality,

[TABLE]

Therefore, by (4.2) and (1.8), we have that

[TABLE]

As $q:=\max\{p_{1},\ldots,p_{m}\}>2$, for any $j\in\{1,\dots,m\}$ such that $p_{j}\geq 2$ we have that $\frac{(p_{j}-2)q}{q-2}<p_{j}$; it thus follows that

[TABLE]

The final term on the right-hand side of this inequality is then bounded analogously as the corresponding term in the proof of Theorem 3.1 and we therefore omit the details (cf. (3.7)):

[TABLE]

Substituting this into the right-hand side of the previous inequality completes the proof of the lemma in the case when $\hat{q}=q>2$.

If $\hat{q}=2$ (and therefore $\hat{q}^{\prime}=2$), then $p_{j}\in(1,2]$ for all $j=1,\ldots,m$. It then follows from (1.8) that

[TABLE]

whereby

[TABLE]

Hence, by Lemma 4.1,

[TABLE]

The final assertion in the statement of the lemma is a consequence of the assumed strong convergence of $u^{h}_{0}$ to $u_{0}$ as $h\rightarrow 0_{+}$. That completes the proof of the lemma. ∎

Next, we prove the required uniform bounds on the sequence $\{u^{h}_{\Delta t}\}$.

Lemma 4.3**.**

There exists a positive constant $c$, independent of $h$ and $\Delta t$, such that, for all $h\in(0,h_{0}]$,

[TABLE]

furthermore, with $\hat{q}:=\max\{q,2\}$,

[TABLE]

and

[TABLE]

Proof.

The first inequality is just a restatement of the bound on the first term in the inequality (4.2) in terms of $u^{h}_{\Delta t}$, and the third inequality is a restatement of the second inequality in terms of $\tilde{u}^{h}_{\Delta t}$. We shall therefore concentrate on the proof of the second stated inequality. The argument proceeds along the same lines as the proof of (3.13). For each $\varphi\in L^{\hat{q}}(0,T;W^{1,\hat{q}}_{0}(\Omega;\mathbb{R}^{m}))$ we have, using the stability in $W^{1,\hat{q}}_{0}(\Omega;\mathbb{R}^{m})$ of the $L^{2}(\Omega;\mathbb{R}^{m})$ orthogonal projector $P^{h}:L^{2}(\Omega;\mathbb{R}^{m})\rightarrow V^{h}_{m}$ (cf. (3.1)), that

[TABLE]

where in the transition to the last line we have used the third and the first bound from Lemma 4.2, together with the assumed strong convergence of $u^{h}_{0}$ to $u_{0}$ in $L^{2}(\Omega;\mathbb{R}^{m})$ to deduce that the constant $c$ is independent of $h$ and $\Delta t$. Dividing through by $\|\varphi\|_{L^{\hat{q}}(0,T;W^{1,{\hat{q}}}(\Omega;\mathbb{R}^{m}))}$, taking the supremum over all $\varphi\in L^{\hat{q}}(0,T;W^{1,{\hat{q}}}_{0}(\Omega;\mathbb{R}^{m}))$ and recalling the definition of the norm of the dual space $L^{\hat{q}^{\prime}}\!(0,T;W^{-1,\hat{q}^{\prime}}\!(\Omega;\mathbb{R}^{m}))=(L^{\hat{q}}(0,T;W^{1,{\hat{q}}}_{0}(\Omega;\mathbb{R}^{m})))^{\prime}$ the stated inequality directly follows. ∎

We shall now discuss the convergence of the sequence of functions $\{u_{\Delta t}^{h}\}$ with $h\in(0,h_{0}]$ fixed.

Lemma 4.4**.**

There exists a subsequence of $\Delta t$, labelled $\Delta t_{k}$, and a function $u^{h}\in C([0,T];V^{h}_{m})$ such that $u_{\Delta t_{k}}^{h}\rightarrow u^{h}$ in $C([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$, and $\partial_{t}u_{\Delta t_{k}}^{h}\rightharpoonup\partial_{t}u^{h}$ in $L^{\hat{q}^{\prime}}\!(0,T;W^{-1,\hat{q}^{\prime}}\!(\Omega;\mathbb{R}^{m}))$ for any (fixed) $h\in(0,h_{0}]$, as $k\rightarrow\infty$.

Proof.

We begin by noting that, for any $\Delta t$ and any $h\in(0,h_{0}]$, the functions $u^{h}_{\Delta t}$ are continuous in both $t$ and $x$. Thanks to the first two bounds in Lemma 4.3 there exists a subsequence of $\Delta t$, labelled $\Delta t_{k}$, and a function $u^{h}\in L^{\infty}(0,T;V^{h}_{m})$ such that $u_{\Delta t_{k}}^{h}\overset{\ast}{\rightharpoonup}u^{h}$ in $L^{\infty}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$, and $\partial_{t}u_{\Delta t_{k}}^{h}\rightharpoonup\partial_{t}u^{h}$ in $L^{\hat{q}^{\prime}}\!(0,T;W^{-1,\hat{q}^{\prime}}\!(\Omega;\mathbb{R}^{m}))$. By Lemma 2.7, $u^{h}\in C_{w}([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$ and $u_{\Delta t_{k}}^{h}\rightarrow u^{h}$ in $C_{w}([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$. As both $u^{h}_{\Delta t}(t,\cdot)$ and $u^{h}(t,\cdot)$ belong to the finite-dimensional space $V^{h}_{m}$ for all $t\in[0,T]$, it follows that $u^{h}\in C([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$ and $u_{\Delta t_{k}}^{h}\rightarrow u^{h}$ in $u^{h}\in C([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$, because $u^{h}\in C_{w}([0,T];V^{h}_{m})=C([0,T];V^{h}_{m})$ and because a sequence converges in $C_{w}([0,T];V^{h}_{m})$ if, and only if, it converges in $C([0,T];V^{h}_{m})$ as $k\rightarrow\infty$ for any (fixed) $h\in(0,h_{0}]$. ∎

Next, we study the convergence of the sequence of functions $\{\tilde{u}_{\Delta t}^{h}\}$ with $h\in(0,h_{0}]$ fixed.

Lemma 4.5**.**

There exists a function $\tilde{u}^{h}\in L^{\infty}(0,T;V^{h}_{m})$ and a subsequence $\Delta t_{k}$ such that

[TABLE]

for any $r\in[1,\infty)$ and fixed $h\in(0,h_{0}]$, as $k\rightarrow\infty$. Furthermore, along this subsequence we have $a(D\tilde{u}_{\Delta t_{k}}^{h})\rightarrow a(D\tilde{u}^{h})$ in $L^{s}(Q_{T};\mathbb{R}^{m\times n})$ for any $s\in[1,\hat{q}^{\prime})$ and $a(D\tilde{u}_{\Delta t_{k}}^{h})\rightharpoonup a(D\tilde{u}^{h})$ in $L^{\hat{q}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$, for any fixed $h\in(0,h_{0}]$, as $k\rightarrow\infty$.

Proof.

Thanks to the first and fourth bound from Lemma 4.2, there exists a $\tilde{u}^{h}\in L^{\infty}(0,T;V^{h}_{m})$ and a subsequence $\Delta t_{k}$ such that $\tilde{u}_{\Delta t_{k}}^{h}\overset{\ast}{\rightharpoonup}\tilde{u}^{h}$ in $L^{\infty}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$ and

[TABLE]

for any fixed $h\in(0,h_{0}]$, as $k\rightarrow\infty$.

By applying Theorem 1 from [9] (with $X=W^{1,p}_{0}(\Omega;\mathbb{R}^{m})$, $B=L^{2}(\Omega;\mathbb{R}^{m})$ and $Y=W^{-1,\hat{q}^{\prime}}\!(\Omega;\mathbb{R}^{m})$ and $r=1$, there), we deduce that $\{\tilde{u}^{h}_{\Delta t_{k}}\}$ is relatively compact in $L^{p}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$, for each fixed $h\in(0,h_{0}]$; therefore, because it is a bounded sequence in $L^{\infty}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$, by interpolation $\{\tilde{u}^{h}_{\Delta t_{k}}\}$ is also relatively compact in $L^{r}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$, for each fixed $h\in(0,h_{0}]$ and all $r\in[1,\infty)$.

Hence, in particular, $\{\tilde{u}^{h}_{\Delta t_{k}}\}$ is relatively compact in $L^{1}(0,T;L^{1}(\Omega;\mathbb{R}^{m}))$, for each fixed $h\in(0,h_{0}]$. As $\tilde{u}^{h}_{\Delta t_{k}}(t,\cdot)\in V^{h}_{m}$ for all $t\in[0,T]$, and $V^{h}_{m}$ is finite-dimensional, for each fixed $h\in(0,h_{0}]$, by norm equivalence in finite-dimensional spaces, it follows that $\{\tilde{u}^{h}_{\Delta t_{k}}\}$ is relatively compact in $L^{1}(0,T;W^{1,1}_{0}(\Omega;\mathbb{R}^{m}))$, for each fixed $h\in(0,h_{0}]$. Therefore $D\tilde{u}^{h}_{\Delta t_{k}}\rightarrow D\tilde{u}^{h}$ in $L^{1}(Q_{T};\mathbb{R}^{m\times n})$ as $k\rightarrow\infty$, for each fixed $h\in(0,h_{0}]$. Thus we can extract a subsequence with respect to $k$ (not indicated) such that $D\tilde{u}^{h}_{\Delta t_{k}}\rightarrow D\tilde{u}^{h}$ a.e. on $Q_{T}$ as $k\rightarrow\infty$, for each fixed $h\in(0,h_{0}]$. Thanks to the continuity of $a$, we then have that $a(D\tilde{u}^{h}_{\Delta t_{k}})\rightarrow a(D\tilde{u}^{h})$ a.e. on $Q_{T}$ as $k\rightarrow\infty$, for each fixed $h\in(0,h_{0}]$. Since by the third inequality from Lemma 4.2 the sequence $\{a(D\tilde{u}^{h}_{\Delta t_{k}})\}$ is weakly compact in $L^{\hat{q}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$ and hence in particular also in $L^{1}(Q_{T};\mathbb{R}^{m\times n})$, for each fixed $h\in(0,h_{0}]$, it follows by Vitali’s convergence theorem that $a(D\tilde{u}^{h}_{\Delta t_{k}})\rightarrow a(D\tilde{u}^{h})$ strongly in $L^{1}(Q_{T};\mathbb{R}^{m\times n})$ as $k\rightarrow\infty$, for each fixed $h\in(0,h_{0}]$, and therefore also strongly in $L^{s}(Q_{T};\mathbb{R}^{m\times n})$ as $k\rightarrow\infty$, for each fixed $h\in(0,h_{0}]$ and any $s\in[1,\hat{q}^{\prime})$. By the uniqueness of the weak limit it then follows that $a(D\tilde{u}_{\Delta t_{k}}^{h})\rightharpoonup a(D\tilde{u}^{h})$ weakly in $L^{\hat{q}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$, for any fixed $h\in(0,h_{0}]$, as $k\rightarrow\infty$. ∎

We now show that the limit functions $u^{h}$ and $\tilde{u}^{h}$ defined above are equal when considered as elements of the space $L^{2}(Q_{T};\mathbb{R}^{m})$.

Lemma 4.6**.**

The limiting functions $u^{h}$ and $\tilde{u}^{h}$ are equal in the space $L^{2}(Q_{T};\mathbb{R}^{m})=L^{2}(0,T;L^{2}(\Omega;\mathbb{R}^{m}))$.

Proof.

We apply the triangle inequality to see that, with $\Delta t_{k}=T/N_{k}$ and $h\in(0,h_{0}]$,

[TABLE]

The first and last terms converge to zero in the limit as $\Delta t_{k}\rightarrow 0_{+}$ thanks to the convergence results stated in Lemma 4.4 and Lemma 4.5 (with $r=2$), respectively, for each fixed $h\in(0,h_{0}]$. We now show that the second term also converges to zero in the limit:

[TABLE]

Thanks to the second inequality in Lemma 4.2, the right-hand side of the last inequality converges to zero in the limit of $\Delta t_{k}\rightarrow 0_{+}$ ($N_{k}\rightarrow\infty$ as $k\rightarrow\infty$). Hence, $u^{h}=\tilde{u}^{h}$, as has been asserted. ∎

Next, we discuss the attainment of the initial condition.

Lemma 4.7**.**

The limiting function $u^{h}$ satisfies the initial condition in the following sense:

[TABLE]

for each fixed $h\in(0,h_{0}]$.

Proof.

Thanks to Lemma 4.4, for any $\epsilon>0$ there is an $\ell\in\mathbb{N}$ such that $\|u^{h}(t,\cdot)-u_{\Delta t_{k_{l}}}^{h}(t,\cdot)\|_{L^{2}(\Omega;\mathbb{R}^{m})}<\epsilon$ for all $l\geq\ell$ and all $t\in[0,T]$. Furthermore, by the triangle inequality:

[TABLE]

for any $l\in\mathbb{N}$. Hence, in particular, for $l=\ell$ (fixed),

[TABLE]

As $u_{\Delta t_{k_{\ell}}}^{h}$ is a continuous function in $t$, passage to the limit in this inequality, using that $\lim_{t\rightarrow 0_{+}}u_{\Delta t_{k_{l}}}^{h}(t,x)=u_{0}^{h}(x)$, implies that, for each $\epsilon>0$,

[TABLE]

Thus we have the asserted result. ∎

Finally, we focus on the convergence of the source term as the time step tends to [math].

Lemma 4.8**.**

We have that $\lim_{\Delta t\rightarrow 0_{+}}\|\tilde{F}_{\Delta t}(\cdot,x)-F(\cdot,x)\|_{L^{r}(0,T)}$ for all $r\in[1,\infty)$ and a.e. $x\in\Omega$.

Proof.

The proof of this result is standard and is therefore omitted. ∎

We are now ready for passage to the limit in (4.4) as $\Delta t\rightarrow 0_{+}$. By Lemma 4.4, Lemma 4.5, and Lemma 4.8, we see that the limit function $u^{h}$ satisfies

[TABLE]

with $u^{h}(0,x)=u^{h}_{0}(x)$ for $x\in\overline{\Omega}$. We can now pass to the limit $h\rightarrow 0_{+}$ by (the proof of) Theorem 3.1, to deduce that as the discretization parameters converge to zero the solution of the fully discrete scheme (4.1) converges to a Young measure solution of the problem (1.1)–(1.9) under consideration.

In preparation for the considerations in the next section, we now discuss continuous dependence of $u^{h}_{\Delta t}$ on the initial data. We start by considering two functions, $u_{0}^{h}\in V^{h}_{m}$ and $v_{0}^{h}\in V^{h}_{m}$, and we let $u_{1}^{h}$ and $v_{1}^{h}$ be respective solutions to (4.1) after a single time step. Subtracting the resulting equations satisfied by $u_{1}^{h}$ from the equation satisfied by $v_{1}^{h}$ and choosing as the test function $\phi^{h}=u_{1}^{h}-v_{1}^{h}$, it follows that

[TABLE]

which after applying Young’s inequality to the right-hand side and absorbing terms, gives

[TABLE]

As the functions $u_{1}^{h}$ and $v_{1}^{h}$ are expressed in terms of a finite basis, we see through the equivalence of norms in finite-dimensional spaces and the first inequality in Lemma 4.3 that $\|Du_{1}^{h}\|_{L^{\infty}(\Omega;\mathbb{R}^{m\times n})}$ and $\|Du_{1}^{h}\|_{L^{\infty}(\Omega;\mathbb{R}^{m\times n})}$ are both bounded by a constant $C=C(h,\|u_{0}\|_{L^{2}(\Omega;\mathbb{R}^{m})},\|F\|_{L^{2}(Q_{T};\mathbb{R}^{m})})$. Therefore we can make use of the local Lipschitz condition satisfied by $a$ and the inverse inequality

[TABLE]

(with $C(h)=Ch^{-1}$ when the family of triangulations $\mathcal{T}^{h}$ is quasiuniform) to deduce, by suppressing in our notation for the constant $C(h,\|a\|_{\rm Lip,loc})$ below its dependence on the (fixed) data $\|u_{0}\|_{L^{2}(\Omega;\mathbb{R}^{m})}$ and $\|F\|_{L^{2}(Q_{T};\mathbb{R}^{m})}$, that

[TABLE]

Therefore by choosing $\Delta t$ sufficiently small for $h\in(0,h_{0}]$ fixed, so that

[TABLE]

we have thus shown the desired continuous dependence on the initial data for a single time step. By iterating this estimate, noting that the constant $C(h,\|a\|_{\rm Lip,loc})$ remains the same regardless of the choice of the time level $i\in\{1,\dots,N\}$, we obtain that, for each $i=1,\ldots,N$, the following estimate holds:

[TABLE]

Recall that $\Delta t$ and $N$ are related via the identity $T=N\Delta t$. Hence,

[TABLE]

To transfer this continuous dependence on the initial data to the limiting functions $u^{h}$ and $v^{h}$ as $\Delta t\rightarrow 0_{+}$, we consider $u^{h}_{\Delta t}$ and $v^{h}_{\Delta t}$, defined, for $i=1,\dots,N$, by

[TABLE]

and, analogously,

[TABLE]

It then follows from (4.5) that

[TABLE]

Let $\Delta t_{k}$ be a sequence such that $u^{h}_{\Delta t_{k}}\rightarrow u^{h}$ and $v^{h}_{\Delta t_{k}}\rightarrow v^{h}$ in $C([0,T];L^{2}(\Omega;\mathbb{R}^{m}))$ and $k\rightarrow\infty$ (cf. Lemma 4.4). Now passing $k\rightarrow\infty$ (whereby $\Delta t_{k}=T/N_{k}\rightarrow 0_{+}$ and $N_{k}\rightarrow\infty$) we see that

[TABLE]

We have therefore shown the following result.

Proposition 4.9**.**

Let all of the assumptions of Theorem 3.1 hold, with the additional assumption that, for $h\in(0,h_{0}]$ fixed, $\Delta t=T/N$ is chosen so that $0<1-C(h,\|a\|_{\rm Lip,loc})\Delta t<1$. Then, there is a unique solution to (3.2a), which depends continuously on the choice of discretized initial data, and which can be constructed numerically by the approximations considered in this subsection.

Remark 4.10**.**

The above proposition in particular implies that the map $S^{h}_{t}$ sending the initial condition $u_{0}^{h}$ to $S^{h}_{t}u_{0}^{h}=u^{h}(t,\cdot)$ is (Lipschitz) continuous from $V^{h}_{m}$ to $V^{h}_{m}$ for all $t\in[0,T]$ and $h\in(0,h_{0}]$ fixed, and by the equivalence of norms in finite-dimensional vector spaces, so is the map sending the initial condition $u_{0}^{h}$ to $Du^{h}(t,\cdot)$. Recall that we are writing the functions $u^{h}$ in terms of a finite element basis

[TABLE]

The above considerations enable us to view the map from $\mathbb{R}^{\mathcal{N}(h)}$ to $\mathbb{R}^{\mathcal{N}(h)}$, which sends the vector $\underline{\alpha}^{h}(0)$ of coefficients of the discretized initial condition in terms of the finite-element basis to the vector $\underline{\alpha}^{h}(t)$ of coefficients of $u^{h}$, in the same finite-element basis, as a continuous function. To see this, let $\underline{\alpha}^{h}$ and $\underline{\beta}^{h}$ be the coefficients of $u^{h}$ and $v^{h}$ respectively, and let

[TABLE]

The (normalized) Gramm matrix $M^{h}=(M_{i,j}^{h})_{i,j=1}^{\mathcal{N}(h)}$ is symmetric positive definite, with $C_{0}I\leq M^{h}\leq C_{1}I$, where $0<C_{0}\leq C_{1}$, in the sense of symmetric positive definite matrices. By expressing

[TABLE]

it follows that,

[TABLE]

[TABLE]

Thus the map $\underline{\alpha}^{h}(0)\in\mathbb{R}^{\mathcal{N}(h)}\mapsto\underline{\alpha}^{h}(t)\in\mathbb{R}^{\mathcal{N}(h)}$ is (Lipschitz) continuous for all $t\in[0,T]$ and $\Delta t$ sufficiently small so as to ensure that $0<1-C(h,\|a\|_{\rm Lip,loc})\Delta t<1$, with $h\in(0,h_{0}]$ fixed.

5. Computation of Young measure solutions

Theorem 3.1 and the results of Section 4 (in particular, Proposition 4.9) give criteria under which Young measure solutions exist as limits of solutions to corresponding semi- or fully-discrete problems. However, the results of Section 3 and Section 4 only give “indirect” information about the Young measure $\nu$ through its action on functions belonging to the space $E_{\frac{p}{\alpha}}$ for $\alpha\in(1,\frac{p}{q-1})$. Indeed, note that for any function $b$ which has the form (1.7)–(1.9), we can copy the calculations used to deduce (3.19) to see that $b(Du^{h})\rightharpoonup\langle\nu,b\rangle$ in $L^{\hat{q}^{\prime}}(Q_{T};\mathbb{R}^{m})$. In this section we shall discuss an algorithm for the numerical approximation of Young measure solutions to systems of the form (1.1), which allows for a more direct approximation of the Young measure $\nu$. The ideas presented here have been inspired by [10] and [11], where the authors develop and analyze numerical schemes for the approximation of measure-valued solutions to systems of hyperbolic conservation laws. We will show that the algorithms developed there can be adapted to systems of forward-backward parabolic equations of the form (1.1) exhibiting Young measure solutions. In particular, we shall demonstrate that, for the class of problems under consideration here, there are Young measure solutions for which the measure $\nu$ may be constructed as a limit of averaged sums of particular Dirac masses, each of whose support is the value taken by an appropriate approximating solution. Throughout this section, we shall be taking subsequences to approximate Young measure solutions, for which we have not shown a uniqueness result: therefore, as we have already noted in the Introduction, different subsequences could potentially converge to different Young measure solutions.

5.1. Preliminary definitions and results

We recall here a selection of definitions and results from Appendix 1 of [10] that are pertinent to the discussion herein, extending them, where necessary.

Definition 5.1**.**

We let $\mathcal{P}(\mathbb{R}^{m})$ denote the set of all probability measures on $\mathbb{R}^{m}$, and define, for $r>0$, the subset $\mathcal{P}^{r}(\mathbb{R}^{m})$ of all probability measures $\mu\in\mathcal{P}(\mathbb{R}^{m})$ for which $\langle\mu,|\xi|^{r}\rangle<\infty$, where $\xi$ denotes the identity function.

We denote by $(\mathfrak{X},\mathcal{F},P)$ a probability space, with $\mathcal{F}$ being a $\sigma$-algebra of sets on the space $\mathfrak{X}$, and $P$ being a probability measure, and we let $u:\mathfrak{X}\times Q_{T}\rightarrow\mathbb{R}^{m}$ be a random field (that is, a jointly measurable function). Given this, we define the law of $u$ as follows:

[TABLE]

for Borel subsets $B\subset\mathbb{R}^{m}$. It is clear that $\mu_{t,x}(\mathbb{R}^{m})=1$ for all $(t,x)\in Q_{T}$. We have the following result.

Lemma 5.2**.**

Let $(\omega;t,x)\in\mathfrak{X}\times Q_{T}\mapsto u(\omega;t,x)\in\mathbb{R}^{m}$ be a random field; then, (5.1) is equivalent to

[TABLE]

for every continuous $r$-component vector function $g$, defined on $\mathbb{R}^{m}$, which is such that $\int_{\mathfrak{X}}|g(u(\omega;t,x))|{\,{\rm d}}P(\omega)$ is finite for a.e. $(t,x)\in Q_{T}$.

Proof.

As the asserted equality is understood componentwise, it suffices to prove it component-by-component. We shall therefore assume in the argument below that $r=1$, i.e., that $g$ is a continuous mapping from $\mathbb{R}^{m}$ into $\mathbb{R}$ such that $\int_{\mathfrak{X}}|g(u(\omega;t,x))|{\,{\rm d}}P(\omega)$ is finite for a.e. $(t,x)\in Q_{T}$. Let us first consider the case when the function $g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ is nonnegative, and take a sequence of simple functions $g_{n}$ defined on $\mathbb{R}^{m}$ by

[TABLE]

where $A^{n}_{k}$, $k=1,\dots,n$, are disjoint measurable subsets of $\mathbb{R}^{m}$, $n$ is a positive integer, and $b^{n}_{k}\in\mathbb{R}$, $k=1,\dots,n$, are such that $g_{n}(\xi)$ increases to $g(\xi)$. This can be achieved by defining

[TABLE]

and

[TABLE]

Observe that $g_{n}(\xi)\leq g(\xi)$ for all $\xi\in\mathbb{R}^{m}$. Further, for a given $\xi\in\mathbb{R}^{m}$ we have that $g_{n}(\xi)=n$ for $g(\xi)\geq n$ whereas for $g(\xi)<n$ we have that

[TABLE]

where $\lfloor y\rfloor$ denotes the greatest integer $\leq y$. As $\lfloor 2^{n+1}g(\xi)\rfloor\geq 2\lfloor 2^{n}g(\xi)\rfloor$, it follows that $g_{n+1}(\xi)\geq g_{n}(\xi)$ for all $\xi\in\mathbb{R}^{m}$. Furthermore, as $\lfloor 2^{n}g(\xi)\rfloor\geq 2^{n}g(\xi)-1$, we have that $g_{n}(\xi)\geq g(\xi)-2^{-n}$ as soon as $n>g(\xi)$. Since $g(\xi)\geq g_{n}(\xi)$ it follows that $\lim_{n\rightarrow\infty}g_{n}(\xi)=g(\xi)$ for all $\xi\in\mathbb{R}^{m}$.

Now, let $\mathcal{F}^{n}_{k}:=\{\omega\in\mathfrak{X}:\ u(\omega;t,x)\in A^{n}_{k}\}$. We observe that

[TABLE]

Hence we have that

[TABLE]

where we have applied the Monotone Convergence Theorem in the first and last lines to exchange the limit with the integral over $\mathbb{R}^{m}$. From here the extension to any, not necessarily nonnegative, function $g:\mathbb{R}^{m}\rightarrow\mathbb{R}$ such that $\int_{\mathfrak{X}}|g(u(\omega;t,x))|{\,{\rm d}}P(\omega)$ is finite for a.e. $(t,x)\in Q_{T}$, is straightforward: we decompose $g$ as $g=g^{+}-g^{-}$, where $g^{+}:=\frac{1}{2}(|g|+g)$ and $g^{-}:=\frac{1}{2}(|g|-g)$ are (respectively) the positive and negative part of $g$, and apply the above argument to $g^{+}$ and $g^{-}$ separately, noting that both are nonnegative on $\mathbb{R}^{m}$. ∎

Remark 5.3**.**

It is easy to extend Lemma 5.2 to functions of the form $g(t,x,u(\omega,t,x))$ where $g$ is a Carathéodory function, i.e., it is continuous in its third variable for a.e. $(t,x)\in Q_{T}$ and measurable in $(t,x)\in Q_{T}$ for each value of the third variable, and $\int_{\mathfrak{X}}|g(t,x,u(\omega,t,x))|\,{\rm d}P(\omega)<\infty$ for a.e. $(t,x)\in Q_{T}$.

We have the following proposition, which shows that the law $\mu=\mu_{t,x}$ defined by (5.1) is in fact a Young measure. This is Proposition 1 in Appendix 1 of [10].

Proposition 5.4**.**

If $u:\mathfrak{X}\times Q_{T}\rightarrow\mathbb{R}^{m}$ is jointly measurable, then the law $\mu=\mu_{t,x}$ (cf. (5.1)) defines a Young measure.

Proof.

The proof is the same as that in [10]. ∎

Finally we consider the following result adapted from [10], found therein as Proposition 2 in Appendix 1.

Proposition 5.5**.**

For every Young measure $\{\nu_{t,x}\}_{(t,x)\in Q_{T}}\in\mathcal{P}^{r}$, $r>0$, on $\mathbb{R}^{m}$ there exists a probability space $(\mathfrak{X},F,P)$ with $P\in\mathcal{P}^{r}$ and a Borel measurable function $u:\mathfrak{X}\times\Omega\rightarrow\mathbb{R}^{m}$ such that $u$ has law given by $\nu$. In particular we may choose the probability space to be the Lebesgue measure on the interval $[0,1)$ with the Borel $\sigma$-algebra.

Proof.

The proof of this result relies on using the characterization given in Lemma 5.2. In [10] a proof is supplied with $g\in C_{0}(\mathbb{R}^{m})$; we extend this result below to cover our particular class of Young measures. We give the proof in the case of scalar-valued functions $u$; the extension to vector-valued functions follows by a component-wise argument.

To extend Proposition 2 in Appendix 1 of [10] from $g\in C_{0}(\mathbb{R})$ to $g\in C(\mathbb{R})$ with $\int_{\mathfrak{X}}|g(u(\omega;t,x))|{\,{\rm d}}P(\omega)<\infty$ for a.e. $(t,x)\in Q_{T}$, let $\varphi\in C^{\infty}_{0}(\mathbb{R})$ be such that $0\leq\varphi\leq 1$, $\varphi(\xi)\equiv 1$ for $|\xi|\leq 1$ and $\varphi(\xi)\equiv 0$ for $|\xi|\geq 2$. Let $g_{k}(x):=g(x)\varphi(x/k)$. Clearly, $g_{k}$ is a sequence of functions in $C_{0}(\mathbb{R})$, $g_{k}(\xi)=g(\xi)$ if $|\xi|\leq k$, and $g_{k}$ converges pointwise to $g$ on $\mathbb{R}$ as $k\rightarrow\infty$. Furthermore, $|g_{k}(x)|\leq|g(x)|$ for all $x\in\mathbb{R}$. Thus, thanks to the Dominated Convergence Theorem, we have that

[TABLE]

We remark that progressing from the left-hand side of the second equality above to the right-hand side of that equality is precisely the result of Proposition 2 in Appendix 1 of [10]. ∎

Finally, we state and prove a standard lemma regarding independent and identically distributed random variables, and their images under measurable functions.

Lemma 5.6**.**

Let $(\mathfrak{X},F,P)$ be a probability space and let $f:\mathfrak{X}\rightarrow\mathbb{R}^{m}$ be a measurable function. If two random variables, $Y_{1}$ and $Y_{2}$, defined on $\mathfrak{X}$ are independent and identically distributed, then the random variables $f(Y_{1})$ and $f(Y_{2})$ are independent and identically distributed on $\mathfrak{X}$.

Proof.

We first prove that the random variables $f(Y_{1})$ and $f(Y_{2})$ are identically distributed. We let $\mathcal{A}\subset\mathbb{R}^{m}$ be a measurable set. It then follows from the measurability of $f$ that

[TABLE]

To show independence, we let $\mathcal{A}_{1}$ and $\mathcal{A}_{2}$ be measurable sets and calculate

[TABLE]

That completes the proof of the lemma. ∎

5.2. An overview of the algorithm

The algorithms discussed below involve generating (at random) a large set of initial data, evolving each of these forward under the solution operator of a semi- or fully-discrete numerical method, and considering the arithmetic average of the resulting functions at a fixed time, before completing various limit passages. The algorithms considered here have been motivated by similar algorithms used in [10] and [11] for the approximation of measure-valued solutions to hyperbolic problems. To this end, let $Y(\Omega,\mathbb{R}^{m})$ denote the set of all Young measures from $\Omega$ to $\mathbb{R}^{m}$.

Algorithm A: Let the initial data for an underlying time-dependent PDE be given as a Young measure $\sigma\in Y(\Omega;\mathbb{R}^{m})$ and let $\Delta$ denote a PDE discretization parameter associated with a certain numerical scheme (in [10] and [11] a finite difference scheme is utilized, and so $\Delta$ would correspond to a vector of grid sizes in the various co-ordinate directions).

Step 1: Let $u_{0}:\omega\in\mathfrak{X}\mapsto u_{0}(\omega;\cdot)\in L^{p}(\Omega;\mathbb{R}^{m})$ be a random field on a probability space $(\mathfrak{X},F,P)$ with law $\sigma_{x}$, meaning that $\sigma_{x}(E)=P(u_{0}(\omega;x)\in E)$ for all Borel sets $E\subset\mathbb{R}^{m}$ and $x\in\Omega$.

Step 2: Evolve the initial random field by applying a suitable numerical scheme, with solution map $\textbf{S}_{t}^{\Delta}$, to the initial data $u_{0}(\omega;\cdot)$ for every $\omega\in\mathfrak{X}$, obtaining an approximate random field $u^{\Delta}(\omega;\cdot,t):=\textbf{S}_{t}^{\Delta}u_{0}(\omega;\cdot)$, $t\in(0,T]$.

Step 3: Define the approximate measure-valued solution $\mu^{\Delta}$ as the law of $u^{\Delta}$ with respect to $P$, that is, for all Borel sets $E\subset\mathbb{R}^{m}$ and $(t,x)\in Q_{T}$,

[TABLE]

Note that we should expect a slight difference compared with [10] and [11] in what we are trying to do, as the Young measures appearing in our context are generated by sequences of gradients. Next, we need a method to approximate the random field $(\omega;\cdot,\cdot)\mapsto u^{\Delta}(\omega;\cdot,\cdot)$, as performing these computations for every $\omega\in\mathfrak{X}$ would be infeasible. This is Algorithm 4.3 in [11].

Algorithm B: Let, as above, $\Delta$ denote a discretization parameter, and let $M\in\mathbb{N}$. Let $\sigma^{\Delta}$ be the initial Young measure.

Step 1: From some probability space $(\mathfrak{X},F,P)$ draw $M$ independent and identically distributed random fields $u_{0}^{\Delta,1},\ldots,u_{0}^{\Delta,M}$ all with the same law $\sigma^{\Delta}$.

Step 2: For each $k\in\{1,\dots,M\}$ and for a fixed $\omega\in\mathfrak{X}$ approximate the solution to the PDE using the solution operator with initial data $u_{0}^{\Delta,k}(\omega)$; denote $u^{\Delta,k}(\omega;\cdot,t):=\textbf{S}_{t}^{\Delta}u_{0}^{\Delta,k}(\omega;\cdot)$.

Step 3: Define the approximate measure-valued solution by

[TABLE]

Results from [11] (see Theorem 4.5, Theorem 5.1 and Corollary 5.4) then guarantee convergence of the sequence of approximate measure-valued solutions to a measure-valued solution as one passes to the limit (diagonally). In this section we discuss in more detail these algorithms and their adaptation to the systems of parabolic PDEs which we are interested in.

5.3. Modifications for systems of forward-backward parabolic PDEs

We now describe the necessary changes to Algorithms A and B described above in order to be able to apply them to the problem (1.1)–(1.9).

Below we formulate the analogues of Algorithm A and Algorithm B, which are needed for our parabolic problem. Since in our case the solution is a Sobolev function, we make here the additional restriction that we only consider atomic initial data: that is, our initial datum is assumed to be given by a function $u_{0}\in L^{2}(\Omega;\mathbb{R}^{m})$, which we view as the atomic Young measure $\delta_{u_{0}}=\delta_{u_{0}(x)}$, $x\in\Omega$.

Algorithm C: Let the initial datum for the problem (1.1)–(1.9) be given as a function $u_{0}\in L^{2}(\Omega;\mathbb{R}^{m})$ and let $h\in(0,h_{0}]$ be the spatial grid size parameter.

Step 1: Let $\upsilon:\mathfrak{X}\rightarrow L^{2}(\Omega;\mathbb{R}^{m})$ be a random field on a probability space $(\mathfrak{X},F,P)$, and discretize $\upsilon$ by a finite element approximation of random fields $\upsilon^{h}\in V^{h}_{m}$, so that $\|\upsilon^{h}(\omega;\cdot)\|_{L^{2}(\Omega;\mathbb{R}^{m})}\leq 1$ for $P$-a.e. $\omega\in\mathfrak{X}$. Discretize $u_{0}$ by a finite element approximation $u_{0}^{h}\in V^{h}_{m}$, and then perturb this discretization by defining

[TABLE]

where $\epsilon\in(0,1]$, and let $\sigma^{h,\epsilon}$ be the law of $u_{0}^{h,\epsilon}$, meaning that $\sigma^{h,\epsilon}_{x}(E)=P(u_{0}^{h,\epsilon}(\omega;x)\in E)$ for $\omega\in\mathfrak{X}$ and $x\in\Omega$.

Step 2: Evolve the initial random field by applying a suitable numerical scheme, with solution map $\textbf{S}_{t}$, to the initial datum $u_{0}^{h,\epsilon}(\omega)$ for every $\omega\in\mathfrak{X}$, obtaining an approximate random field $u^{h,\epsilon}(\omega,t,\cdot):=\textbf{S}_{t}u_{0}^{h,\epsilon}(\omega;\cdot)$, $t\in(0,T]$.

Step 3: Define the approximate Young measure solution $\mu^{h,\epsilon}$ as the law of $u^{h,\epsilon}$ with respect to $P$, that is, for all Borel sets $E\subset\mathbb{R}^{m}$,

[TABLE]

and define, analogously, $\nu^{h,\epsilon}$ to be the law of $Du^{h,\epsilon}$ with respect to $P$.

As before we need a method to approximate the random field $u_{0}^{\epsilon}(\omega;x)$, as well as the measures $\mu^{h,\epsilon}$ and $\nu^{h,\epsilon}$.

Algorithm D: Let the initial datum for the problem (1.1)–(1.9) be given as a function $u_{0}\in L^{2}(\Omega;\mathbb{R}^{m})$, let $h\in(0,h_{0}]$ be the spatial grid size parameter, and let $M\in\mathbb{N}$.

Step 1: From some probability space $(\mathfrak{X},F,P)$ draw $M$ independent and identically distributed random fields $\upsilon^{h,1},\ldots,\upsilon^{h,M}$ such that $\|\upsilon^{h,i}(\omega;\cdot)\|_{L^{2}(\Omega;\mathbb{R}^{m})}\leq 1$, for $P$-a.e. $\omega\in\mathfrak{X}$ and for all $i=1,\ldots,M$, and such that $u_{0}^{h,1,\epsilon},\ldots,u_{0}^{h,M,\epsilon}$ all have the same law $\sigma^{h,\epsilon}$, where

[TABLE]

with $\epsilon\in(0,1]$. We make explicit here the fact that $\|u_{0}^{h,k,\epsilon}(\omega;\cdot)\|_{L^{2}(\Omega;\mathbb{R}^{m})}\leq\|u_{0}^{h}\|_{L^{2}(\Omega;\mathbb{R}^{m})}+\epsilon$ for all $\omega\in\mathfrak{X}$, so that $\|u_{0}^{h,k,\epsilon}(\omega;\cdot)\|_{L^{2}(\Omega;\mathbb{R}^{m})}$ can be assumed to be bounded independent of $k$, $\epsilon$, and $h$ (thanks to the assumed strong convergence of $u_{0}^{h}$ to $u_{0}$ in $L^{2}(\Omega;\mathbb{R}^{m})$ and because $\epsilon\in(0,1]$).

Step 2: For each $k$ and for a fixed $\omega\in\mathfrak{X}$ approximate the solution to the initial boundary-value problem (1.1)–(1.9) under consideration using the solution operator with initial data $u_{0}^{h,k,\epsilon}(\omega)$; denote $u^{h,k,\epsilon}(\omega;t,\cdot)=\textbf{S}_{t}u_{0}^{h,k,\epsilon}(\omega;\cdot,\cdot)$.

Step 3: Define the approximate Young measure solution $(\mu_{t,x}^{h,M,\epsilon},\nu_{t,x}^{h,M,\epsilon})$ by

[TABLE]

5.4. Applying the algorithm to construct Young measure solutions

In what follows we will show that Algorithm D described above converges to a Young measure solution of the system of PDEs (1.1)–(1.9) under consideration as we pass $\Delta t\rightarrow 0_{+}$, $M\rightarrow\infty$, $h,\epsilon\rightarrow 0_{+}$ (in that order).

Remark 5.7**.**

The parameter $\omega$ taken from the probability space $\mathfrak{X}$ is intended to represent the “seed” in the numerical algorithm. Different seeds $\omega\in\mathfrak{X}$ will give rise to different approximate solutions. What will be shown in our analysis of the limit process $\Delta t\rightarrow 0_{+}$, $M\rightarrow\infty$, $h,\epsilon\rightarrow 0_{+}$ is that one can pass to the limit in $M$ (along a subsequence) to obtain convergence to a quantity which is independent of $\omega\in\mathfrak{X}$, for $P$-a.e. $\omega\in\mathfrak{X}$.

Step 1: Passage to the limits $\Delta t\rightarrow 0_{+}$, $M\rightarrow\infty$. We shall use the semidiscrete numerical scheme discussed in Section 3, augmented with the time-stepping procedure discussed in Section 4. By Lemma 4.9 we have that, for any $\omega\in\mathfrak{X}$, any $\epsilon>0$, any positive integer $k$ and any $h\in(0,h_{0}]$ and $\Delta t=T/N$, $N\in\mathbb{N}$, we can construct, via the time-stepping procedure from Section 4, a function $u^{h,k,\epsilon}_{\Delta t}(\omega;\cdot,\cdot)$, which for any $h\in(0,h_{0}]$ and any $\Delta t>0$ (sufficiently small, for a given fixed $h\in(0,h_{0}]$) uniquely solves the fully discrete scheme and depends continuously on the discretized initial condition $u_{0}^{h,k,\epsilon}$. For the rest of the analysis we shall assume for the sake of brevity that the passage to the limit $\Delta t\rightarrow 0_{+}$ has already been made for $\omega\in\mathfrak{X}$ and $h\in(0,h_{0}]$ fixed, by repeating the analysis performed in Section 4, and we shall therefore consider the limiting function $u^{h,k,\epsilon}(\omega;\cdot,\cdot)$ resulting from the passage to the limit $\Delta t\rightarrow 0_{+}$ instead of $u^{h,k,\epsilon}_{\Delta t}(\omega;\cdot,\cdot)$. Passage to the limit $\Delta t\rightarrow 0_{+}$, for reasons that will be explained in due course, must happen before the limit passage $M\rightarrow\infty$. In particular we remark that as we are at this stage considering finitely many perturbations of the discretized initial condition, the existence of a suitable subsequence $\Delta t_{j}\rightarrow 0_{+}$, which works for each $k\in\{1,\ldots,M\}$, is straightforward.

For each initial datum $u_{0}^{h,k,\epsilon}(\omega;\cdot)$ (for fixed $\omega\in\mathfrak{X}$ and $k\in\{1,\ldots,M\}$), we have from (3.7) the following energy estimate:

[TABLE]

Furthermore, (3.11) implies, with $\hat{q}:=\max\{q,2\}$ and $\hat{q}^{\prime}=\frac{\hat{q}}{\hat{q}-1}$, that

[TABLE]

and (3.13) yields

[TABLE]

where $C$ is a positive constant that is independent of $h,k$ and $\epsilon$. This follows from the form of the perturbed initial condition in Algorithm D, which guarantees that the perturbed initial data are bounded uniformly in $h,k$ and $\epsilon$. Summing the equations satisfied by $u^{h,k,\epsilon}(\omega;\cdot,\cdot)$ over $k=1,\ldots,M$ and dividing by $M$ we see that the functions $u^{h,k,\epsilon}$ satisfy, for all $\psi^{h}\in L^{q}(0,T;V^{h}_{m})$, the equality:

[TABLE]

Integration by parts in time then yields, for all $\psi^{h}\in W^{1,1}_{0}(0,T;V^{h}_{m})$,

[TABLE]

We now wish to let $M\rightarrow\infty$ in (5.6). To this end, we need to prove certain convergence results that will enable passage to this limit. We begin by considering the sequence of functions $\{\frac{1}{M}\sum_{k=1}^{M}a(Du^{h,k,\epsilon}(\omega;\cdot,\cdot))\}_{M\geq 1}$, for a fixed $\omega\in\mathfrak{X}$, which is, thanks to (3.11), bounded in the $L^{\hat{q}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$ norm, uniformly in $h\in(0,h_{0}]$ and $\epsilon\in(0,1]$, and as such there is a subsequence $\{M_{n}\}_{n\geq 1}$ and a function $\chi^{h,\epsilon}(\omega;\cdot,\cdot)\in L^{\hat{q}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$ such that $\frac{1}{M_{n}}\sum_{k=1}^{M_{n}}a(Du^{h,k,\epsilon}(\omega;\cdot,\cdot))\rightharpoonup\chi^{h,\epsilon}(\omega;\cdot,\cdot)$ in $L^{\hat{q}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$, as $n\rightarrow\infty$. We aim to identify this function $\chi^{h,\epsilon}(\omega;\cdot,\cdot)$ for $\omega\in\mathfrak{X}$, by identifying the limit in some weaker space and appealing to the uniqueness of weak limits.

Let us consider the initial datum $u^{h,\epsilon}_{0}(\omega;x)=u^{h}_{0}(x)+\epsilon\upsilon^{h}(\omega;x)$ that is independent of the initial data $u_{0}^{h,1,\epsilon},\ldots,u_{0}^{h,M,\epsilon}$, and has the same distribution as these initial data (meaning in particular that the law of $u_{0}^{h,\epsilon}$ is given by $\sigma^{h,\epsilon}$), and satisfies $\|\upsilon^{h}(\omega;\cdot)\|_{L^{2}(\Omega;\mathbb{R}^{m})}\leq 1$ for $P$-almost every $\omega\in\mathfrak{X}$. We evolve this according to the algorithm for each $\omega\in\mathfrak{X}$ to obtain a function $u^{h,\epsilon}(\omega;t,x)$. We denote the law of $u^{h,\epsilon}$ by $\mu^{h,\epsilon}$, and similarly denote the law of $Du^{h,\epsilon}$ by $\nu^{h,\epsilon}$. From Theorem 3.1 we see that the following regularity results hold:

[TABLE]

Analogously to (5.2)–(5.4), the norms of these functions in the respective spaces can be bounded by a positive constant $C$, independent of $h$ and $\epsilon$. Next, we formulate an analogue of Theorem 11 from [10].

Lemma 5.8**.**

Let $\alpha\in(1,\frac{p}{q-1})$ and suppose that $\gamma\in L^{\alpha^{\prime}}\!(Q_{T};E_{\frac{p}{\alpha}})$ (understood as an $m\times n$ matrix-valued function). Along the subsequence $\{M_{n}\}_{n\geq 1}$ we have

[TABLE]

in $L^{2}(\mathfrak{X})$ as $n\rightarrow\infty$. In particular, for each $\gamma\in L^{\alpha^{\prime}}\!(Q_{T};E_{\frac{p}{\alpha}})$ there is a subsequence of $M_{n}$ such that for $P$-a.e. $\omega\in\mathfrak{X}$ the convergence is pointwise.

Proof.

Given a random field $\eta:\mathfrak{X}\rightarrow L^{1}(Q_{T})$ we define its expectation with respect to the probability measure $P$ as

[TABLE]

Similarly to [10] we define, for $k=1,\ldots,M$, $i=1,\dots,m$ and $j=1,\dots,n$ the following quantities:

[TABLE]

and

[TABLE]

We need to show that $G^{1}_{i,j},\ldots,G^{M_{n}}_{i,j}$ are independent and identically distributed.

We begin by remarking that as the initial data $u_{0}^{h,1,\epsilon},\ldots,u_{0}^{h,M_{n},\epsilon}$ are independent and identically distributed (in the variable $\omega\in\mathfrak{X}$), and as the solution operator $\textbf{S}_{t}$ mapping an initial condition $\tilde{u}_{0}$ to $\tilde{u}(t,\cdot)$ is a continuous map (see Lemma 4.9), the functions $u^{h,k,\epsilon}$ are all independent and identically distributed as well (see Lemma 5.6). Furthermore, by the remarks following the statement of Theorem 4.9, the functions $Du^{h,k,\epsilon}$ are also independent and identically distributed.

Recall that the functions $\gamma_{i,j}$ are measurable, they are continuous in their third argument, and satisfy the growth rate coming from the space $E_{\frac{p}{\alpha}}$. This implies that the functions $\gamma_{i,j}(t,x,Du^{h,k,\epsilon}(\omega;t,x))$ are independent and identically distributed with respect to $\omega\in\mathfrak{X}$. From Tonelli’s Theorem it then follows that the functions $G^{k}_{i,j}$ are measurable in $\omega$, and thus by Lemma 5.6 are also independent and identically distributed in $\omega\in\mathfrak{X}$.

We now compute the $L^{2}$ error with respect to the probability measure $P$ of the following quantity:

[TABLE]

and as in [10], using the fact that the random variables $G_{i,j}^{k}$ are independent and identically distributed, we can reduce this to

[TABLE]

Given $\alpha\in(1,\frac{p}{q-1})$ we have that

[TABLE]

Next, we can bound the norm $\|1+|Du^{h,\epsilon}(\omega)|^{\frac{p}{\alpha}}\|_{L^{\alpha}(Q_{T})}$ independently of $h$, $\epsilon$ and $\omega$ by using the energy estimate (3.7) (analogously to (5.3)). Thus, since $P$ is a probability measure, we have that $\mathbb{E}((G_{i,j})^{2})$ is bounded, and so the right-hand side of (5.10) converges to zero as $M_{n}\rightarrow\infty$.

We wish to apply Lemma 5.2 (note also Remark 5.3). In order to do this we must show that

[TABLE]

Note that if we show that

[TABLE]

then the desired bound is true for almost every $(t,x)\in Q_{T}$, which is sufficient for our purposes. We see, using that $\gamma\in L^{\alpha^{\prime}}\!(Q_{T};E_{\frac{p}{\alpha}})$ with $\alpha\in(1,\frac{p}{q-1})$, that

[TABLE]

which is finite by the assumptions on $\gamma$ and by (5.7). Then, by applying Lemma 5.2, we see that:

[TABLE]

Repeating these calculations for all components $\gamma_{i,j}$ gives us the convergence of

[TABLE]

in $L^{2}(\mathfrak{X})$ as $M_{n}\rightarrow\infty$, from which we deduce the existence of a subsequence (not relabelled) converging for $P-$ almost every $\omega\in\mathfrak{X}$ using standard results in measure theory. ∎

What we would like to conclude is that by the above theorem we have, for $P-$almost every $\omega\in\mathfrak{X}$, the weak-star convergence of $\nu^{h,M_{n},\epsilon}_{\omega}=\nu^{h,M_{n},\epsilon}_{\omega,\cdot,\cdot}$ to $\nu^{h,\epsilon}=\nu^{h,\epsilon}_{\cdot,\cdot}$ as $n\rightarrow\infty$. This does not follow immediately, as a-priori there is no reason why the subsequence of $M_{n}$ along which we have pointwise convergence should be independent of the function $\gamma$ (a point not clarified in [10] or [11]). To get around this we note that as $E_{\frac{p}{\alpha}}$ is separable (see [17]), the space $L^{\alpha^{\prime}}\!(Q_{T};E_{\frac{p}{\alpha}})$ is also separable for our choice of $\alpha$. Therefore there is a countable dense subset $\{\gamma_{i}\}_{i=1}^{\infty}$. The desired weak-star convergence then follows from a diagonal argument.

Recall from the calculations performed prior to Lemma 5.8 that, as $M_{n}\rightarrow\infty$ ($n\rightarrow\infty$),

[TABLE]

for every $\psi\in L^{\hat{q}}([0,T];W^{1,\hat{q}}_{0}(\Omega;\mathbb{R}^{m}))$, and for all $\omega\in\mathfrak{X}$. Furthermore, by Lemma 5.8, there is a subsequence of $M_{n}$ such that, as $M_{n}\rightarrow\infty$ ($n\rightarrow\infty$),

[TABLE]

for every $\psi\in L^{\alpha^{\prime}}\!(0,T;W^{1,\alpha^{\prime}}_{0}\!(\Omega;\mathbb{R}^{m}))$, and for $P$-a.e. $\omega\in\mathfrak{X}$.

As both (5.11) and (5.12) hold for $\psi\in L^{\max\{\hat{q},\alpha^{\prime}\}}(0,T;W^{1,\max\{\hat{q},\alpha^{\prime}\}}_{0}(\Omega;\mathbb{R}^{m}))$, by the uniqueness of the weak limit we can identify $\chi^{h,\epsilon}(\omega;\cdot,\cdot)=\langle\nu^{h,\epsilon},a\rangle$ for $P$-a.e. $\omega\in\mathfrak{X}$. As the right-hand side of this equality is independent of $\omega\in\mathfrak{X}$, it follows that $\chi^{h,\epsilon}$ is also independent of $\omega$, and we shall therefore suppress the dependence on $\omega$ in our notation and write $\chi^{h,\epsilon}(t,x)$ instead of $\chi^{h,\epsilon}(\omega;t,x)$, noting that $\chi^{h,\epsilon}=\langle\nu^{h,\epsilon},a\rangle\in L^{\hat{q}^{\prime}}\!(Q_{T};\mathbb{R}^{m\times n})$. Thus, by combining (5.11) and (5.12), which guarantee weak convergence in $L^{\hat{q}^{\prime}}\!(0,T;W^{1,\hat{q}^{\prime}}_{0}\!(\Omega;\mathbb{R}^{m}))$ and $L^{\alpha}(0,T;W^{1,\alpha}_{0}(\Omega;\mathbb{R}^{m}))$, respectively, for $P$-a.e. $\omega\in\mathfrak{X}$, and therefore, by uniqueness of the weak limit, also in the function space $L^{\max\{\hat{q}^{\prime},\alpha\}}(0,T;W^{1,\max\{\hat{q}^{\prime},\alpha\}}_{0}(\Omega;\mathbb{R}^{m}))$, for $P$-a.e. $\omega\in\mathfrak{X}$, we have that

[TABLE]

for all $\psi\in L^{\min\{\hat{q},\alpha^{\prime}\}}(0,T;W^{1,\min\{\hat{q},\alpha^{\prime}\}}_{0}(\Omega;\mathbb{R}^{m}))$, for $P$-a.e. $\omega\in\mathfrak{X}$.

The terms involving the measures $\mu^{h,\epsilon}$ are treated completely analogously. Thus, by passing to the limit $M_{n}\rightarrow\infty$ ($n\rightarrow\infty$) in (5.6), we have that for all $\psi^{h}\in W^{1,1}_{0}(0,T;V^{h}_{m})$ the following identity is satisfied:

[TABLE]

Step 2: Passing $h,\epsilon\rightarrow 0_{+}$. By Lemma 5.8 and noting that $a\in E_{\frac{p}{\alpha}}$, it follows that $\langle\nu^{h,M_{n},\epsilon}_{\omega,\cdot,\cdot},a\rangle\rightharpoonup\langle\nu^{h,\epsilon}_{\cdot,\cdot},a\rangle$ weakly in $L^{\alpha}(Q_{T};\mathbb{R}^{m\times n})$, for $\alpha\in(1,\frac{p}{q-1})$, and for $P$-a.e. $\omega\in\mathfrak{X}$. Thus, by weak lower-semicontinuity of the norm function, we have that $\|\langle\nu^{h,\epsilon}_{\cdot,\cdot},a\rangle\|_{L^{\alpha}(Q_{T};\mathbb{R}^{m\times n})}\leq\liminf_{n\rightarrow\infty}\|\langle\nu^{h,M_{n},\epsilon}_{\omega,\cdot,\cdot},a\rangle\|_{L^{\alpha}(Q_{T};\mathbb{R}^{m\times n})}$ for $P$-a.e. $\omega\in\mathfrak{X}$.

Let $1<\alpha<\frac{p}{q-1}$ and consider the function $\tilde{a}(\xi)=1+|\xi|^{\frac{p}{\alpha}}$; then, $\tilde{a}\in E_{p/\alpha}$. It then follows from Lemma 5.8 that $\langle\nu^{h,M_{l},\epsilon}_{\omega,\cdot,\cdot},\tilde{a}\rangle\rightharpoonup\langle\nu^{h,\epsilon},\tilde{a}\rangle$ in $L^{\alpha}(Q_{T};\mathbb{R}^{m\times n})$ as $l\rightarrow\infty$ for $P$-a.e. $\omega\in\mathfrak{X}$. Then, by weak lower-semicontinuity of the norm function $\|\cdot\|_{L^{\alpha}(Q_{T};\mathbb{R}^{m\times n})}$, we have that

[TABLE]

where $C$ is a positive constant, independent of $h\in(0,h_{0}]$ and $\epsilon\in(0,1]$. Here, in the transition from the left-hand side of the second inequality to its right-hand side we have used the triangle inequality; in the transition from the left-hand side of the third inequality to its right-hand side we have again used the triangle inequality; in the transition from the left-hand side of the fourth inequality to its right-hand side we have used the energy estimate (5.2), followed by recalling that $\|u_{0}^{h,k,\epsilon}(\omega)\|_{L^{2}(\Omega;\mathbb{R}^{m})}\leq\|u_{0}^{h}\|_{L^{2}(\Omega;\mathbb{R}^{m})}+\epsilon$, with $\epsilon\in(0,1]$, and noting that by the assumed strong convergence of $u_{0}^{h}$ to $u_{0}$ in $L^{2}(\Omega;\mathbb{R}^{m})$ one can bound $\|u_{0}^{h}\|_{L^{2}(\Omega;\mathbb{R}^{m})}$ by a constant, independent of $h\in(0,h_{0}]$.

Having shown that $\|\langle\nu^{h,\epsilon},\tilde{a}\rangle\|_{L^{\alpha}(Q_{T})}$ is bounded by a constant, independent of $h$ and $\epsilon$, we now compute

[TABLE]

which we can further bound from above by a positive constant $C$, independent of $h\in(0,h_{0}]$ and $\epsilon\in(0,1]$, using (5.14). Hence, for each $\epsilon\in(0,1]$ there is a subsequence of the sequence of measures $\{\nu^{h,\epsilon}\}$ indexed by $h_{j}$ and a measure $\nu^{\epsilon}\in L^{\alpha}(Q_{T};E_{p/\alpha}^{\prime})$ such that for every $\gamma\in E_{p/\alpha}$ we have that

[TABLE]

in $L^{\alpha}(Q_{T};\mathbb{R}^{m\times n})$ as $j\rightarrow\infty$. In particular this holds for $\gamma=a\in E_{p/\alpha}$, $1<\alpha<\frac{p}{q-1}$, that is:

[TABLE]

in $L^{\alpha}(Q_{T};\mathbb{R}^{m\times n})$ for any such $\alpha>1$. Similar calculations hold for the family of measures $\mu^{h,\epsilon}$.

Taking the final limit in $\epsilon$ is now straightforward, and is done similarly to the above. In particular we note that we can in fact take a diagonal subsequence of $(h,\epsilon)$ and pass to the limit along this subsequence, obtaining the existence of a two measures $\mu$ and $\nu$ satisfying, for all $\psi\in W^{1,1}_{0}(0,T;W^{1,\alpha^{\prime}}_{0}(\Omega;\mathbb{R}^{m}))$ (contained in $W^{1,1}_{0}(0,T;W^{1,\min\{\hat{q},\alpha^{\prime}\}}_{0}(\Omega;\mathbb{R}^{m}))\subset L^{\min\{\hat{q},\alpha^{\prime}\}}(0,T;W^{1,\min\{\hat{q},\alpha^{\prime}\}}_{0}(\Omega;\mathbb{R}^{m}))$):

[TABLE]

Step 3: Interpretation as a Young measure solution. Next, we discuss how the pair of measures $(\mu,\nu)$ satisfying (5.15) can be can be interpreted as a Young measure solution. Recall that $\mu^{h,\epsilon}$ is the law of the random variable $u^{h}=u^{h,\epsilon}(\omega;t,x)$, and that, for each $\omega\in\mathfrak{X}$ and each $t\in[0,T]$, the function $u^{h,\epsilon}(\omega;t,\cdot)\in V^{h}_{m}\subset W^{1,\infty}_{0}(\Omega;\mathbb{R}^{m})$ is a linear combination of finite element basis functions over the domain $\Omega$. Therefore,

[TABLE]

is weakly differentiable with respect to $x\in\Omega$, and thanks to the definition of weak derivative and Fubini’s theorem,

[TABLE]

Using the weak-star convergence of $\mu^{h,\epsilon}$ and $\nu^{h,\epsilon}$ to, respectively, $\mu$ and $\nu$ established in Step 2, in conjunction with the regularity of $u^{h,\epsilon}$ coming from (5.7)–(5.9), and the fact that $\mathfrak{X}$ is a space of finite $P$-measure, we see that by passing $h,\epsilon\rightarrow 0_{+}$ along subsequences, we get the existence of a function

[TABLE]

such that, for $(t,x)\in(0,T]\times\Omega$,

[TABLE]

and

[TABLE]

This comes from the fact that, for example,

[TABLE]

and we have that

[TABLE]

for all $\phi\in L^{1}(\mathfrak{X};L^{\alpha}(0,T;W_{0}^{1,\alpha}(\Omega;\mathbb{R}^{m})))$. However, as $\mathfrak{X}$ has finite $P$-measure, we can choose test functions which are constant in the parameter $\omega$ to see that

[TABLE]

for all $\phi\in L^{\alpha}(0,T;W_{0}^{1,\alpha}(\Omega;\mathbb{R}^{m}))$, which gives the desired representation for $\langle\mu_{t,x},\xi\rangle$ (the corresponding term for $\langle\nu_{t,x},\xi\rangle$ follows similarly).

Now we define the function

[TABLE]

The function $U$ is locally integrable over $Q_{T}$ and it therefore has a well-defined distributional derivative $DU$. If we can demonstrate that $DU(t,x)=\langle\nu_{t,x},\xi\rangle$ in a suitable sense, then we will have shown that the pairing $(U,\nu)$ is a Young measure solution of the system (1.1)–(1.9) under consideration. What must first be shown is that it makes sense to speak of $DU$ as a function, rather than as a distribution. To that end, note that the function $u\in L^{p}(0,T;W^{1,p}_{0}(\Omega;\mathbb{R}^{m}))$ satisfies $u|_{\partial\Omega}=0$, and so we can extend the function $u$ (and hence the function $U$) by zero from $[0,T]\times\overline{\Omega}$ to the whole of $[0,T]\times\mathbb{R}^{n}$ so that the extended function (still denoted by $u$) belongs to $L^{p}(0,T;W^{1,p}(\mathbb{R}^{n};\mathbb{R}^{m}))$. We consider the difference quotient of $U$, defined by

[TABLE]

where $e_{i}$ is the unit vector in the $i$-th co-ordinate direction. We then have that

[TABLE]

where in the transition to the last line we have used that, by Jensen’s inequality,

[TABLE]

Therefore, by a standard characterization of Sobolev functions in terms of difference quotients, $D_{i}U(t,\cdot)\in L^{p_{i}}(\Omega;\mathbb{R}^{m\times n})$, $i=1,\dots,n$, for almost every $t\in[0,T]$. By integrating over $t\in[0,T]$ in the above inequality it follows that $D_{i}U\in L^{p_{i}}(Q_{T};\mathbb{R}^{m\times n})$, $i=1,\dots,n$, and in particular $U\in L^{p}(0,T;W^{1,p}_{0}(\Omega;\mathbb{R}^{m}))$, as required. By an analogous argument,

[TABLE]

thanks to the regularity results

[TABLE]

which follow from the regularity properties of $u^{h,\epsilon}$ stated in (5.7)–(5.9) using the weak lower-semi-continuity of the norm function.

Next, we write

[TABLE]

and let $\varphi\in C^{\infty}_{0}(Q_{T};\mathbb{R}^{m\times n})$ be a test function. We then compute, using Fubini’s Theorem, with $D=D_{x}$ and $\mbox{div}=\mbox{div}_{x}$,

[TABLE]

Hence, $DU(t,x)=V$, and because $V:=\langle\nu,\xi\rangle$, it follows that $DU=\langle\nu,\xi\rangle$, as desired. Thus we have shown that the pair $(U,\nu)$ is a Young measure solution of the system (1.1)–(1.9) under consideration.

Step 4: Attainment of the initial condition. As

[TABLE]

and

[TABLE]

it follows from Lemma 2.7 that there is a subsequence $(h_{j},\epsilon_{j})$ of $(h,\epsilon)$ along which

[TABLE]

uniformly in $C([0,T])$ for all $\varphi\in L^{2}(\Omega;\mathbb{R}^{m})$. Therefore to show that $U(0,\cdot)=u_{0}$ we let $\varphi\in L^{2}(\Omega;\mathbb{R}^{m})$ and compute:

[TABLE]

The term $I_{j}$ converges to zero as $j\rightarrow\infty$ by the uniform convergence in $C([0,T])$ stated in (5.16). The term $II_{j}(t)$ converges to zero as $t\rightarrow 0_{+}$ as it holds that $u^{h_{j},\epsilon_{j}}(\omega;t,x)\rightarrow u^{h_{j},\epsilon_{j}}(\omega,0,x)$ by continuity of $t\in[0,T]\mapsto u^{h_{j},\epsilon_{j}}(\omega;t,x)$, $\omega\in\mathfrak{X}$, $x\in\Omega$, and we can apply the Dominated Convergence Theorem to interchange this limit with the integral over $\mathfrak{X}$ that appears in the definition of $U^{h_{j},\epsilon_{j}}$. The term $III_{j}(t)$ also converges to zero as $t\rightarrow 0_{+}$ as we have that the numerical functions $u^{h_{j},\epsilon_{j}}$ satisfy the initial condition, and we may once again apply the Dominated Convergence Theorem to interchange the limit with the integral over $\mathfrak{X}$. Finally for the term $IV_{j}$ we compute:

[TABLE]

Integrating this over $\Omega$ we get

[TABLE]

which converges to zero as $j\rightarrow\infty$ by the assumed strong convergence of our discretized initial condition $u_{0}^{h}$ to $u_{0}$ in $L^{2}(\Omega;\mathbb{R}^{m})$. In summary then, $\lim_{j\rightarrow\infty}I_{j}=0$, $\lim_{t\rightarrow 0_{+}}II_{j}(t)=0$ for all $j\geq 1$, $\lim_{t\rightarrow 0_{+}}III_{j}(t)=0$ for all $j\geq 1$, and $\lim_{j\rightarrow\infty}IV_{j}=0$. Passing to the limit $t\rightarrow 0_{+}$ in (5.17) with $j\geq 1$ kept fixed and then with $j\rightarrow\infty$ we deduce that

[TABLE]

Thus we have that $U(0,x)=u_{0}(x)$ for almost every $x\in\Omega$, showing that the initial condition is satisfied.

In summary, we have that that the pair $(U,\nu)$ is a Young measure solution to the model problem satisfying the initial condition $U(0,\cdot)=u_{0}(\cdot)$. For clarity we state the following theorem which collects the most relevant properties proved in Steps 1–4 above.

Theorem 5.9**.**

Algorithms C and D converge to a Young measure solution of the problem (1.1)–(1.9) in the event that the initial Young measure is a Dirac mass, $\delta_{u_{0}(x)}$, concentrated on an initial datum $u_{0}\in L^{2}(\Omega;\mathbb{R}^{m})$.

Remark 5.10**.**

If we have a (nonatomic) measure-valued initial datum $\sigma$, Steps 1–4 above provide an outline of an algorithm for proving the existence of a Young measure solution, assuming that all of the $M$ “perturbations” to $\sigma$ can be chosen to be bounded, independent of $M$. This requires being able to find a random variable $u_{0}$ whose law is given by the initial Young measure $\sigma$. The existence of such a random variable is guaranteed by Proposition 5.5.

Remark 5.11**.**

In [10] and [11], in their analysis of semidiscrete schemes the authors are able to interchange the order in which the limits of the spatial discretization parameter and the parameter $\epsilon$ tending to zero and $M\rightarrow\infty$ are taken. For us here the passages to the limits in $h\rightarrow 0_{+}$, $\epsilon\rightarrow 0_{+}$ and $\Delta t\rightarrow 0_{+}$ are all interchangeable, but the passage in $M$ is problematic from this point of view. The reason for this is that, while the limits obtained in the proof of Theorem 5.8 are independent of $h$, computing these limits in the first place required us to use the independence and the “identicalness” of the distributions of the functions $u^{h,k,\epsilon}$. This required that the solution operator $\textbf{S}_{t}$ was measurable. However, for the model problem under consideration here, we were only able to show continuity of this operator on the functions $u^{h,k,\epsilon}$ for each fixed $h$, and were unable to transfer this property to the limiting function as we pass $h\rightarrow 0_{+}$, meaning that within the context of the algorithms, we must take a limit passage in $M$ before the limit passage in $h$. If instead one could show that the operator $\textbf{S}_{t}$ was measurable when considered as a mapping from the nondiscretized initial condition $u_{0}\in L^{2}(\Omega;\mathbb{R}^{m})$, we would then be able to interchange the limits in any way we desired to, as this property would be sufficient for the independence and “identicalness” of the distributions to carry over to the limiting functions as we take $h\rightarrow 0_{+}$.

Acknowledgement. We are grateful to Professor Mike Cullen (UK Met Office, Exeter) for stimulating discussions and for suggesting to us the problem studied in this paper. This work was supported by the UK Engineering and Physical Sciences Research Council [EP/L015811/1].