Numerical Approximation of Young Measure Solutions to Parabolic Systems of Forward-Backward Type
Miles Caddick Endre Süli


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





Mathematical Institute, University of Oxford, Andrew Wiles Bldg., Woodstock Road, Oxford OX2 6GG, UK
miles.caddick@maths.ox.ac.uk
Mathematical Institute, University of Oxford, Andrew Wiles Bldg., Woodstock Road, Oxford OX2 6GG, UK
endre.suli@maths.ox.ac.uk



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 pa­rabolic systems of forward-backward type of the form ∂_t u - (a(Du)) + Bu = F, where B ∈ℝ^m × m, Bv · v ≥ 0 for all v ∈ℝ^m, F is an m-component vector-function defined on a bounded open Lipschitz domain Ω⊂ℝ^n, and a is a locally Lipschitz mapping of the form a(A)=K(A)A, where K : ℝ^m × n→ℝ. 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.
 



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




[
    Endre Süli
    December 30, 2023
=====================





§ 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:


    ∂_tu-div(a(Du))+Bu=F,


with u:[0,T]×Ω→m, where Ω⊂n is a bounded open Lipschitz domain. Here Du denotes the spatial gradient of u, with (Du)_ik= ∂ u_i/∂ x_k, i=1,…,m, k=1,…,n. Throughout we will take real numbers p_i, for i=1,…,m, that satisfy


    p_i>max{1, 2n/n+2}.


We define


    p:=min_i=1,…,m p_i          q:=max_i=1,…,m p_i,


and require that


    q-p<1.


We shall also assume that


    F∈ L^p'(Q_T;m)∩ L^2(Q_T;m),


with 1/p + 1/p'=1, where Q_T:=(0,T) ×Ω, and that


    Bv· v≥ 0    ∀  v ∈ℝ^m,


where B∈m× m is a constant matrix. All of our results extend directly to the case where B ∈ L^∞(Q_T;m× m), without assuming (<ref>), 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 ∈m× m is a constant matrix satisfying (<ref>). The system (<ref>) is supplemented with the initial condition


    u(0,·)=u_0(·)∈ L^2(Ω;m),


and the homogeneous Dirichlet boundary condition


    u|_(0,T]×∂Ω=0.


The function a:ℝ^m × n→ℝ^m × n is assumed to be a locally Lipschitz mapping of the form


    a(A)=K(A)A,


where K:m × n→ is continuous, and there exist positive constants c_0 and c_1, with c_0 ≤  c_1, such that


    c_0∑_i=1^m(μ_i^2+|A_i|^2)^p_i-2/2≤ K(A)
    ≤ c_1∑_i=1^m(μ_i^2+|A_i|^2)^p_i-2/2    ∀  A ∈ℝ^m × n,


with p_i>max{1,2n/n+2} for all i=1,…,m, and μ_i is a constant defined by


    μ_i
    ∈  if p_i≥2,
    ≠0    if 1<p_i<2,


and where |A_i| denotes the Euclidean norm of the i-th row of the matrix A. More generally, |·| 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 (<ref>)–(<ref>) one can assume that |a(x)| ≤ c_0(1+|x|^q-1) and a(x)· x ≥ c_1|x|^p, with p>max{1,2n/n+2} and p≤ q < p+1.

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


In <cit.> 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, Ω = (0,1), and the function K, referred to as a stability function,  is defined as follows:


    K(∂_xu)={[ (|∂_xu_1|^2+|∂_xu_2|^2)^3/2/|∂_xu_1|^2+|∂_xu_2|^2+∂_xu_3                                             if ∂_xu_3>0,;                                                                                                                  ;                          √(|∂_xu_1|^2+|∂_xu_2|^2-∂_xu_3)                                             if ∂_xu_3≤0. ].


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))^ T is


    ∂_tu-∂_x(K(∂_xu)∂_xu)+ Bu   =F      for (t,x)∈(0,T]×(0,1),
    
    u(t,1)   =f(t)   for t∈(0,T],
    
    u_1(t,0)   =g_1(t)  for t∈(0,T],
    
    u_2(t,0)   =g_2(t)  for t∈(0,T],
    ∂_xu_3(t,0)   =h(t)   for t∈(0,T],


subject to the initial condition


    u(0,x)=u_0(x)    for x∈(0,1).


Here, u:[0,T]× [0,1]→3 is a vector with components u=(u_1,u_2,u_3)^ T, B is the skew-symmetric matrix


    B=[  0  1  0; -1  0  0;  0  0  0 ]


and F is a forcing term given by the constant vector


    F(t,x)=[ -V;  U;  0 ].


What is immediately seen is that, given the above definition of the function K, the product K(∂_xu)∂_xu lacks a coercivity estimate in a Sobolev space on the range ∂_xu_3>0, because the mapping ξ∈3↦ a(ξ) := K(ξ)ξ∈3 is not coercive in the range ξ_3>0; indeed, for ξ=(1,0,ξ_3)^ T with ξ_3>0, a(ξ)·ξ/|ξ| → 1 as ξ_3 → +∞. Another difficulty arises from the fact that the function ξ↦ a(ξ):=K(ξ)ξ is not monotone on 3, i.e., it is not true that (a(ξ)-a(η)):(ξ - η) ≥ 0 for all ξ, η∈3, and therefore monotone operator theory cannot be applied in this setting.
To see that the mapping ξ↦ a(ξ):=K(ξ)ξ is not monotone on 3, take ξ=(0.035,0,-0.01)^ T and η=(0.05,0,0)^ T; then, computing (a(ξ)-a(η))·(ξ-η) will a give a negative number. Similarly taking ξ=(-0.2, -0.1, 0.2)^ T and η=(-0.1,0,0.5)^ 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(θ)=(1+cos(θ),0,sin(θ)), θ∈ [0,2π). One can show that this is nonzero, and so a is not a conservative vector-field, and thus it is not a potential.





For n=1 and m=3 consider the function K given by ξ∈3↦ K(ξ)=√(ξ_1^2+ξ_2^2+|ξ_3|). Then we are in the regime ∂_x u_3 ≤ 0 from (<ref>). One has the existence of constants c_0 and c_1 such that


    c_0(|ξ_1|+|ξ_2|+√(|ξ_3|))≤ K(ξ)≤ c_1(|ξ_1|+|ξ_2|+√(|ξ_3|))


through bounding and using the equivalence of norms on 3. This function K is therefore covered by our assumptions.
We are unfortunately unable to say anything about the region from (<ref>) where ∂_xu_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 (<ref>)–(<ref>) 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.



We say that a pair (u,ν), where u is a function with


    u ∈ L^∞(0,T;L^2(Ω;m))∩ L^p(0,T;W^1,p_0(Ω;m))     Du_i ∈ L^p_i(0,T;L^p_i(Ω;m× n)),   i=1,…,m,


and

    ∂_tu∈ L^q̂'(0,T;W^-1,q̂'(Ω;m)),


and (t,x) ∈ Q_T ↦ν_t,x is an m × n-valued Young measure, such that


    ⟨ν,a⟩∈ L^q̂'(Q_T;m× n),


is called a Young measure solution to the problem described by (<ref>), (<ref>) and (<ref>) if


    ∫_0^T⟨∂_tu,ϕ⟩_(W_0^1,q̂(Ω;m),W^-1,q̂'(Ω;m))+∫_Ω⟨ν_t,x,a⟩ : Dϕ + Bu·ϕxt=∫_0^T∫_Ω F·ϕxt


for all ϕ∈ L^q̂(0,T;W_0^1,q̂(Ω;m)), where q̂:=max{q,2}, q̂'=q̂/q̂-1, and the initial condition is satisfied in the sense that u(0,x)=u_0(x) for almost every x∈Ω.



 The proofs and the main results of this paper can be easily adapted to the case when, instead of assumptions
(<ref>)–(<ref>), the nonlinearity
a satisfies a(A)· A≥ c_1|A|^p,  |a(A)|≤ c_2|A|^q-1+c_3, the function a is locally Lipschitz continuous, (<ref>) holds, and q-p<1, max{1, 2n/n+2} < p ≤  q < ∞.



Measure-valued solutions and Young measure solutions to parabolic partial differential equations have been studied by a number of authors. Frehse & Specovius-Neugebauer, in <cit.>, 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


    a=∇Φ


for some C^1 function Φ. Our terminology Young measure solution follows that of Frehse & Specovius-Neugebauer in <cit.>; 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 <cit.>: under the assumptions that m=1, p=q=2 and that (<ref>) is satisfied,  a sequence of approximating solutions was constructed by means of minimizing an integral, which involves the convexification of Φ. 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 <cit.> and some of the references contained therein, in particular <cit.>, 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 (<ref>) holds, and assumes more restrictive growth conditions on a than those considered here. Even more recently, Kim & Yan in <cit.> have expanded their earlier work in <cit.> 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 <cit.> 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 (<ref>)), neither of the methods in <cit.> or <cit.> will be suitable in our setting. In <cit.> 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 <cit.>.

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 Ω. 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 <cit.> 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 <ref> we will show the existence of a large-data global-in-time Young measure solution to the system (<ref>)–(<ref>) under consideration by performing a spatial finite element approximation and passing to the limit in the spatial discretization parameter. In Section <ref>, 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 <ref> and Section <ref> only yield to us indirect information about the measure ν through its action on functions belonging to a particular function space. This is the motivation behind Section <ref>, in which, stimulated by recent contributions by Fjordholm et al. <cit.> and <cit.>, 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 ν, and will prove that the algorithm converges to a Young measure solution of the problem under consideration.




§ PRELIMINARY DEFINITIONS AND RESULTS




 §.§ Function Spaces


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


    E_r:={g∈ C(m × n;m): lim_|A|→∞g(A)/1+|A|^r exists},


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


    g_E_r=sup_A∈m × n|g(A)|/1+|A|^r.




The following statements hold true for r > 0:



  * The space E_r is separable;

  * The space E_r is metrizable.






For the first property we refer to Kinderlehrer & Pedregal <cit.>. The second follows by taking the metric induced by the norm ·_E_r.





We remark here that E_r is a subspace of the larger space


    Ẽ_r:={g∈ C(m × n;m):sup_A∈m × n|g(A)|/1+|A|^r<∞},


which is not separable. For example, under the assumption (<ref>), a ∈Ẽ_q-1, and for any
α∈ (1,p/q-1) we have that 0<q-1 < p/α<p, and therefore a ∈ E_p/α.



As C_0(m × n;m)⊂ E_r for all r > 0, we have the following nesting of dual spaces for all r > 0:


    E_r'⊂ℳ(m × n;m)= C_0(m × n;m)'.


Since the dual of a normed space is a Banach space, we have that E_r' 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' has a weakly-∗ convergent subsequence.




 §.§ Young Measures


We recall here from <cit.> the so-called Fundamental Theorem for Young Measures, together with some related remarks.



Let Ω⊂n be Lebesgue measurable, let 𝔎⊂m be closed and let
z_j:Ω→m, j=1,2,…, be a sequence of Lebesgue measurable functions satisfying


    lim_j→∞|{x∈Ω: z_j(x)∉ U}|=0


for any open neighbourhood U of 𝔎 in m. Then, there exists a subsequence {z_j_k} of {z_j} and a family of positive measures on m, {ν_x}_x∈Ω, depending measurably on x∈Ω, such that:



  * ν_x_M:=∫_m dν_x≤1 for almost every x∈Ω;

  * supp ν_x⊂𝔎 for almost every x∈Ω; and

  * f(z_j_k)∗⇀⟨ν_x,f⟩=∫_mf(λ)ν_x(λ) in L^∞(Ω)
for each continuous function f:m→ satisfying lim_|λ|→∞f(λ)=0.


Suppose further that {z_j_k} satisfies the following tightness condition:


    lim_s→∞sup_k|{x∈Ω∩ B_R(0): |z_j_k(x)|≥ s}|=0,


for every R>0. Then, we have that ν_x_M=1 for almost every x∈Ω, and, for any measurable subset A of the set Ω,


    f(z_j_k)⇀⟨ν_x,f⟩


in L^1(A) for any continuous function f:m→ such that {f(z_j_k)} is sequentially weakly relatively compact in L^1(A).





If Ω is bounded and the sequence {z_j} is bounded in L^p(Ω;m) for some p∈(1,∞), then we obtain from the theorem the existence of a family of probability measures {ν_x}_x ∈Ω and a subsequence z_j_k such that


    f(z_j_k)⇀⟨ν,f⟩


in L^r(Ω), where f:m→ is a continuous function satisfying


    |f(λ)|≤ c(1+|λ|^s)


with 1<r<p/s and s>0.





The assumption that


    lim_j→∞|{x∈Ω: z_j(x)∉ U}|=0


for any open neighbourhood U of 𝔎 is only used to prove that (2) in Theorem <ref> holds, and is not needed to obtain any of its other conclusions. The proof of this theorem from <cit.> is still valid if we take 𝔎=m, only it then provides no information about the support of the Young measure ν. Therefore, if we are not concerned with the support of the Young measure ν we may always take 𝔎=m in order to apply this theorem.






 §.§ Miscellaneous Results

We recall the following result from Strauss <cit.> (cf. also Lions & Magenes <cit.>, Lemma 8.1, Ch. 3, Sec. 8.4).


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

    L^∞(0, T;X) ∩ C_w([0, T];Y) = C_w([0, T];X),


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.



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



  * {u^j} is bounded in L^∞(0,T; L^2(Ω;m))
and


    u^j∗⇀ u in L^∞(0,T; L^2(Ω;m));



  * {∂_tu^j} is bounded in L^s'(0,T;W^-1,s'(Ω;m)) for some s ∈ (1,∞),
and


    ∂_tu^j⇀∂_tu in L^s'(0,T;W^-1,s'(Ω;m)).




Then, u^j_k, u ∈ C_w([0,T];L^2(Ω;m), and there is a subsequence {u^j_k} such that


    ∫_Ω u^j_k(t,x)· w(x)x →∫_Ω u(t,x)· w(x)x


uniformly in C([0,T]) for all w∈ L^2(Ω;m) as k →∞.



 We shall apply the Arzelà–Ascoli theorem for sequences of uniformly bounded and equicontinuous functions. For w ∈ L^2(Ω;m)∩ W_0^1,s(Ω;m), we let f^j:[0,T]→ be the family of functions defined by


    f^j(t)=∫_Ω u^j(t,x)· w(x)x.


Since ∂_tu^j ∈ L^s'(0,T;W^-1,s'(Ω;m)), it automatically follows that u^j ∈ C([0,T];W^-1,s'(Ω;m)), whereby also u^j ∈ C_w([0,T];W^-1,s'(Ω;m)). As u^j ∈  L^∞(0,T; L^2(Ω;m)), it then follows from Lemma <ref> with X = L^2(Ω;m) and Y = L^2(Ω;m) + W^-1,s'(Ω;m) = (L^2(Ω;m) ∩ W^1,s_0(Ω;m))', that u^j ∈  C_w([0,T]; L^2(Ω;m)). We note here that the duality (L^2(Ω;m) ∩ W^1,s_0(Ω;m))' = L^2(Ω;m) + W^-1,s'(Ω;m) is a consequence of the Duality Theorem (cf. Theorem 2.7.1 on p.32 in Bergh & Löfström <cit.>), because L^2(Ω;m) ∩ W^1,s_0(Ω;m) is dense in both L^2(Ω;m) and W^1,s_0(Ω;m) (e.g., because C^∞_0(Ω;m) is dense in both L^2(Ω;m) and W^1,s_0(Ω;m)). Hence, f^j ∈ C([0,T]).  We further have that


    |f^j(t)|≤∫_Ω |u^j(t,x)||w(x)|x
    ≤u^j_L^∞(0,T;L^2(Ω;m))w_L^2(Ω;m)≤ Cw_L^2(Ω;m),


which gives the boundedness of the sequence {f^j} in C([0,T]) for each w∈ L^2(Ω;m)∩ W_0^1,s(Ω;m). In order to verify equicontinuity of the sequence, note that, for all w ∈ L^2(Ω;m)∩ W_0^1,s(Ω;m),


    |f^j(t+h)-f^j(t)|=|∫_Ω(u^j(t+h,x)-u^j(t,x))· w(x)x|
    =|∫_t^t+h∫_Ω∂_tu^j(s,x)· w(x)xs|
    ≤∂_tu^j_L^s'(0,T;W^-1,s'(Ω;m))|h|^1/sw_W^1,s(Ω;m)≤ C|h|^1/sw_W^1,s(Ω;m).


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


    lim_k→∞∫_0^Tφ(t) ∫_Ω u^j_k(t,x)· w(x)xt
    =lim_k→∞∫_0^T φ(t) f^j_k(t)t=∫_0^T φ(t) f(t)t    ∀ φ∈ C^∞([0,T]).


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


    lim_k→∞∫_0^Tφ(t)∫_Ω u^j_k(t,x)· w(x)xt=∫_0^Tφ(t) ∫_Ω u(t,x)· w(x)xt     ∀ φ∈ C^∞([0,T]),


allowing us to identify


    f(t)=∫_Ω u(t,x)· w(x)x     ∀ t ∈ [0,T].


Thus we have shown that, for all w ∈ L^2(Ω;m)∩ W_0^1,s(Ω;m),


    ∫_Ω u^j_k(t,x)· w(x)x→∫_Ω u(t,x)· w(x)x,    .


As L^2(Ω;m)∩ W_0^1,s(Ω;m) is dense in L^2(Ω;m), it then follows that, for all w ∈ L^2(Ω;m),


    ∫_Ω u^j_k(t,x)· w(x)x→∫_Ω u(t,x)· w(x)x,    ,


which in particular means that for all t∈[0,T] we have u^j_k(t,·)⇀ u(t,·) in the space L^2(Ω;m).


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


Let u:[0,T]×Ω→m, where Ω⊂n is open and bounded, satisfy


    u∈ L^p(0,T;W_0^1,p(Ω;m))∩ L^∞(0,T;L^2(Ω;m))


for some p≥1. Then,


    u∈ L^p(n+2)/n(Q_T;m).








§ SPATIAL DISCRETIZATION AND CONVERGENCE OF THE SEMIDISCRETE PROBLEM


In this section we prove the existence of solutions to the system (<ref>)–(<ref>). 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 Ω is a Lipschitz polytope in n.
The proof presented below can be replicated on any bounded Lipschitz domain Ω by replacing the finite element basis, consisting of continuous piecewise linear basis functions satisfying a homogeneous Dirichlet boundary condition on ∂Ω by an abstract Galerkin basis. Suppose further that 0<h_0≪(Ω) and that {𝒯_h}_h ∈ (0,h_0] is a shape-regular family of subdivisions 𝒯_h of Ω into closed simplexes Δ, where h=max_Δ∈𝒯_h. Consider the finite element space


    V^h:={v^h∈ W^1,∞_0(Ω) : v^h|_Δ is affine for all Δ∈𝒯_h}.


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 {𝒯_h}_h ∈ (0,h_0] is such that the L^2(Ω;m) orthogonal projector P^h onto V^h_m is stable in W^1,q̂_0(Ω;m); i.e., there exists a positive constant C, independent of h ∈ (0,h_0], such that


    D(P^hφ) _L^q̂(Ω;m× n)≤ C D φ_L^q̂(Ω;m × n)    ∀ φ∈ W^1,q̂_0(Ω;m),  .


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

Let {ϕ^h_i}_i=1^𝒩(h) be a basis for V^h_m. We shall therefore seek an approximate solution u_h ∈ V^h_m in the form


    u^h(t,x)=∑_i=1^𝒩(h)α^h_i(t)ϕ^h_i(x),     t ∈ [0,T],    x ∈Ω,


satisfying



    ( ∂_tu^h(t),v^h) + ( a(Du^h(t)),Dv^h)+ (Bu^h,v^h)= ( F,v^h)
        ∀  v^h ∈ V^h_m,


for all t∈[0,T], and


    u^h(0)=u^h_0,



with u^h_0∈ V^h_m and u^h_0→ u_0 (strongly) in L^2(Ω;m) as h → 0_+. The system given by (<ref>) is a system of ODEs for the coefficients α^h_i, i=1,…, 𝒩(h), and a solution to (<ref>), (<ref>) exists locally by Peano's Theorem on some interval [0,T_h) ⊂ [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 ∈ (0,h_0], and that as h → 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.



Suppose that Ω⊂n is a Lipschitz polytope, let u_0∈ L^2(Ω;m), assume that


    p:=min_i=1,…, mp_i>max{2n/n+2,1}     q:=max_i=1,…,mp_i


satisfy q-p<1, let F∈ L^p'(Q_T;m)∩ L^2(Q_T;m), and suppose that a:m × n→m × n is a locally Lipschitz mapping satisfying the assumptions (<ref>)–(<ref>). Then, there exists a Young measure solution (u,ν) of the problem (<ref>) with data given by (<ref>) and (<ref>). Furthermore, there exists a subsequence (not indicated) of solutions {u^h} to the semidiscrete problem (<ref>), (<ref>) such that


    u^h∗⇀u in L^∞(0,T;L^2(Ω;m)),    
    Du^h_i ⇀ Du_i  in L^p_i(0,T; L^p_i(Ω;m× n)),   i=1,…,m.



    a(Du^h)⇀⟨ν,a⟩ in L^q̂'(Q_T;m × n),
        ∂_tu^h⇀∂_tu in L^q̂'(0,T;W^-1,q̂'(Ω;m)),
      ,



    u^h → u  in C_w([0,T];L^2(Ω;m)).





As was noted above, (<ref>), (<ref>) is an initial-value problem for a system of ODEs for the coefficients α^h_i, and a solution exists locally by Peano's Theorem on some interval [0,T_h) ⊂ [0,T], thanks to the assumed (local Lipschitz) continuity of a. We begin by showing that we can extend the numerical solution u^h ∈ V^h_m, defined on [0,T_h) ×Ω up to time T for all h ∈ (0,h_0], so that it is defined on the whole of [0,T]×Ω. To this end, we take v^h = u^h in (<ref>); hence,


    ∫_0^τ∫_Ω∂_tu^h· u^h+K(Du^h)|Du^h|^2
    +Bu^h· u^hxt=∫_0^τ∫_Ω F· u^hxt    ∀ τ∈ (0,T_h),


which, by using (<ref>) and (<ref>) on the left-hand side and Hölder's inequality on the right-hand side, yields


    u^h(τ,·)_L^2(Ω;m)^2+∫_0^τ∫_Ω[∑_i=1^m(μ_i^2+|Du_i^h|^2)^p_i-2/2]|Du^h|^2xt
           ≤ c [u_0^h_L^2(Ω;m)^2
    + (∫_0^TF^p'_L^p'(Ω;m)t)^1/p'(∫_0^τu^h_L^p(Ω;m)^p t)^1/p]    ∀ τ∈ (0,T_h).


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


    |Du_i^h|^p_i≤ (μ_i^2+|Du_i^h|^2)^p_i-2/2 |Du^h|^2    ,


and the following inequality holds for all s ∈ [0,∞) if 1<p_i <2  and μ_i ≠ 0:


    s^p_i≤{[     2^1-p_i/2(μ^2_i + s^2)^p_i -2/2s^2                                      ,; (μ_i^2 + s^2)^p_i -2/2s^2  + |μ_i|^p_i                                      . ].


Using these, we deduce that


    u^h(τ,·)_L^2(Ω;m)^2+∫_0^τ∫_Ω∑_i=1^m
    |Du_i^h|^p_ixt
           ≤ c [1 + u_0^h_L^2(Ω;m)^2
    + (∫_0^TF^p'_L^p'(Ω;m)t)^1/p'(∫_0^τu^h_L^p(Ω;m)^p t)^1/p]    ∀ τ∈ (0,T_h).


Thus, using Young's inequality and Poincaré's inequality in W^1,p_0(Ω;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:


    sup_t∈[0,T_h)u^h(t,·)_L^2(Ω;m)^2+∫_0^T_h∫_Ω∑_i=1^m|Du_i^h|^p_i dxt≤ c(1+u_0^h_L^2(Ω;m)^2+F_L^p'(Q_T;m)^p'),


which provides uniform bounds on the sequence {u^h}, thanks to the strong convergence of u_0^h→ u_0 in L^2(Ω;m) as h → 0_+, which implies in particular that u_0^h_L^2(Ω;m) is uniformly bounded in h. Hence, the right-hand side of (<ref>) is bounded, independent of h ∈ (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 ∈ (0,h_0]. Hence,


    sup_t∈[0,T)u^h(t,·)_L^2(Ω;m)^2+∫_0^T∫_Ω[∑_i=1^m(μ_i^2+|Du_i^h|^2)^p_i-2/2]|Du^h|^2xt
    ≤ c(1 + u_0_L^2(Ω;m)^2 +F_L^p'(Q_T;m)^p')


and also


    sup_t∈[0,T)u^h(t,·)_L^2(Ω;m)^2+∫_0^T∫_Ω∑_i=1^m|Du_i^h|^p_i dxt≤ c(1+u_0_L^2(Ω;m)^2+F_L^p'(Q_T;m)^p').


The uniform bound (<ref>) implies the existence of a function


    u ∈ L^∞(0,T;L^2(Ω;m))∩ L^p(0,T;W^1,p_0(Ω;m))      Du_i ∈ L^p_i(0,T;L^p_i(Ω;m× n)),   i=1,…,m,


for which, up to a subsequence, we have


    u^h∗⇀u in L^∞(0,T;L^2(Ω;m) )


and


    Du^h_i ⇀ Du_i  in L^p_i(0,T; L^p_i(Ω;m× n)),   i=1,…,m.




Furthermore, thanks to (<ref>) and because the Euclidean norm of the vector (|Du_1^h|,…,|Du_m^h|)^ T is bounded by its 1-norm, we have the following bound:


    |a(Du^h)|   ≤ c_1 ∑_i=1^m(μ_i^2+|Du_i^h|^2)^p_i-2/2(∑_k=1^m|Du_k^h|^2)^1/2≤ c_1 ∑_i=1^m∑_k=1^m(μ_i^2+|Du_i^h|^2)^p_i-2/2|Du_k^h|.


Now we seek to bound each of the terms appearing in the last sum in (<ref>). Clearly for the case i=k we have an expression of the form (μ_i^2+|Du_i^h|^2)^p_i-1/2 appearing, which is bounded in L^p_i'(Q_T), uniformly with respect to h ∈ (0,h_0], by the energy estimate (<ref>). In the case i≠ k we proceed by using Hölder's inequality together with (<ref>) and (<ref>), our objective being to bound all of these terms in L^p_i'(Q_T) if p_i≥ 2 or in L^2(Q_T) if 1<p_i<2. We begin with the terms for which p_i≥2:


    ∫_0^T∫_Ω[(μ_i^2+|Du_i^h|^2)^p_i-2/2|Du_k^h|]^p_i'xt
    =∫_0^T∫_Ω (μ_i^2+|Du_i^h|^2)^p_i'(p_i-2)/4|Du_k^h|^p_i' (μ_i^2+|Du_i^h|^2)^p_i'(p_i-2)/4xt
           ≤(∫_0^T∫_Ω(μ_i^2+|Du_i^h|^2)^p_i-2/2|Du_k^h|^2x  dt)^p_i'/2 (∫_0^T∫_Ω(μ_i^2+|Du_i^h|^2)^p_i/2x  dt)^1-p_i'/2
           ≤ c(1+u_0_L^2(Ω;m)^2+F_L^p'(Q_T;m)^p')^p_i'/2(1 + u_0_L^2(Ω;m)^2
    + F^p'_L^p'(Q_T;m))^1-p_i'/2.


Now for the terms in which 1<p_i<2, thanks to (<ref>), we have


    ∫_0^T∫_Ω   [(μ_i^2+|Du_i^h|^2)^p_i-2/2|Du_k^h|]^2xt
       =∫_0^T∫_Ω (μ_i^2+|Du_i^h|^2)^p_i-2/2(μ_i^2+|Du_i^h|^2)^p_i-2/2|Du_k^h|^2xt
       ≤|μ_i|^p_i-2∫_0^T∫_Ω (μ_i^2+|Du_i^h|^2)^p_i-2/2|Du_k^h|^2xt
       ≤ c|μ_i|^p_i-2(1+u_0_L^2(Ω;m)^2+F_L^p'(Q_T;m)^p').


Thus we have bounded each term in the sum appearing on the right-hand side of (<ref>) in L^p_i'(Q_T) (if p_i≥2) or in L^2(Q_T) (if 1<p_i<2), uniformly with respect to h ∈ (0,h_0], and therefore {a(Du^h)}_0<h≤ h_0 is a bounded sequence in L^q̂'(Q_T;m × n). More precisely, there exists a positive constant C, independent of h, such that


    a(Du^h)_L^q̂'(Q_T;m × n)≤ C    ∀  h ∈ (0,h_0],


where we recall that q̂=max{q,2}, whereby q̂'= min{q',2}. Therefore, there is a subsequence (still indexed only by h) and a χ∈ L^q̂'(Q_T;m × n) such that


    a(Du^h)⇀χ    .



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


    ∂_tu^h_L^q̂'(0,T; W^-1,q̂'(Ω;m))≤ C    ∀  h ∈ (0,h_0].


Indeed, for each φ∈ L^q̂(0,T;W^1,q̂_0(Ω;m)), by (<ref>), (<ref>) and (<ref>), we have that


    ∫_0^T (∂_tu^h, φ )  t = ∫_0^T ( ∂_tu^h, P^hφ )  t = - ∫_0^T ( a(Du^h(t)),DP^hφ)  t + ∫_0^T ( F - Bu^h ,P^hφ )  t
           ≤ Ca(Du^h)_L^q̂'(Q_T;m × n)Dφ_L^q̂(Q_T;m× n) + (F_L^2(Q_T;m) + |B| u^h_L^2(Q_T;m)) φ_L^2(Q_T;m)
           ≤ C (Dφ_L^q̂(Q_T;m× n) + φ_L^q̂(Q_T;m))  ≤  C φ_L^q̂(0,T;W^1,q̂(Ω;m)),


and then (<ref>) follows by noting that the dual space of L^q̂(0,T;W^1,q̂_0(Ω;m)) is L^q̂'(0,T;W^-1,q̂'(Ω;m)).

Now let ψ∈ L^q̂(0,T; W_0^1,q̂(Ω;m)) be the strong limit, as h → 0_+ and l →∞, of the sequence of functions


    ψ^𝒩(h),l(t,x)=∑_i=1^𝒩(h)β^l_i(t)ϕ^h_i(x)


in the L^q̂(0,T; W_0^1,q̂(Ω;m)) norm, with β^l_i∈ C([0,T]) and ϕ^h_i∈ V^h_m, i=1,…, 𝒩(h), l=1,2,…; such a limit necessarily exists by density. We then have that


    ∫_0^T∫_Ω[∂_tu^h·ψ+a(Du^h) : Dψ+Bu^h·ψ-F·ψ]xt
       = ∫_0^T∫_Ω∂_tu^h·ψ^𝒩(h),l+a(Du^h) :  Dψ^𝒩(h),l+Bu^h·ψ^𝒩(h),l-F·ψ^𝒩(h),lxt
           +∫_0^T∫_Ω[∂_tu^h· (ψ-ψ^𝒩(h),l)+a(Du^h) : (Dψ-Dψ^𝒩(h),l)
    +Bu^h·(ψ-ψ^𝒩(h),l)-F·(ψ-ψ^𝒩(h),l)]xt
       ≤∂_tu^h_L^q̂'(0,T;W^-1,q̂'(Ω;m))ψ-ψ^𝒩(h),l_L^q̂(0,T;W_0^1,q̂(Ω;m))
           +a(Du^h)_L^q̂'(0,T;L^q̂'(Ω;m × n))Dψ-Dψ^𝒩(h),l_L^q̂(0,T;L^q̂(Ω;m × n))
           +F_L^p'(Q_T;m)ψ-ψ^𝒩(h),l_L^q̂(0,T;W_0^1,q̂(Ω;m))
    +cu^h_L^q̂'(Q_T;m)ψ-ψ^𝒩(h),l_L^q̂(0,T;W_0^1,q̂(Ω;m))
       ≤ c(u_0_L^2(Ω;m),μ_i, T, |Ω|,n,F) (ψ-ψ^𝒩(h),l_L^q̂(0,T;W_0^1,q̂(Ω;m))
    +Dψ-Dψ^𝒩(h),l_L^q̂(0,T;L^q̂(Ω;m × n))).


The right-hand side of this inequality converges to zero in the limit as h→0_+ (and therefore 𝒩(h)→∞) and l →∞. Thus we have that


    lim_h → 0_+∫_0^T∫_Ω∂_tu^h·ψ+a(Du^h) : Dψ + Bu^h·ψxt
    =∫_0^T∫_Ω F·ψxt,


for all ψ∈ L^q̂(0,T; W_0^1,q̂(Ω;m)).


The uniform in h bounds established in (<ref>) and (<ref>) now allow passage to the limit in the first and third term on the left-hand side of (<ref>): by the uniform in h bound on u^h_L^∞(0,T;L^2(Ω;m)), there exists a subsequence (still indexed by h) which converges weakly-∗ to u in L^∞(0,T;L^2(Ω;m)). This enables us to also pass to the limit in the third term of the integrand on the left-hand side of (<ref>). In addition, thanks to (<ref>), one can extract a further subsequence (which we continue to index by h) for which, by the uniqueness of the weak limit,


    ∂_tu^h⇀∂_tu  in   L^q̂'(0,T;W^-1,q̂'(Ω;m)),


meaning that


    lim_h→0_+∫_0^T (∂_tu^h, ψ) t=∫_0^T ⟨∂_tu,ψ⟩t,


for all ψ∈ L^q̂(0,T;W^1,q̂_0(Ω;m)),
where ⟨·,·⟩ is the duality pairing between
W^-1,q̂'(Ω;m) and W^1,q̂_0(Ω;m). Hence,


    ∫_0^T⟨∂_tu , ψ⟩ t +lim_h → 0_+∫_0^T∫_Ωa(Du^h) : Dψ
    + ∫_0^T∫_Ω Bu·ψxt
    =∫_0^T∫_Ω F·ψxt,


for all ψ∈ L^q̂(0,T;W^1,q̂_0(Ω;m)). The final convergence result in the statement of the theorem, that u^h → u in C_w([0,T];L^2(Ω;m)), is a direct consequence of (<ref>), (<ref>) and Lemma <ref>.


In order to pass to the limit in the second term in (<ref>), we shall invoke Theorem <ref>. To this end, we begin by verifying that the subsequence {Du^h_k} satisfies the tightness condition from Theorem <ref>; i.e., that


    lim_s→∞sup_k|{(t,x)∈([0,T]×Ω)∩ Q_R(0): |Du^h_k(t,x)|≥ s}|=0,


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


    A_k,s,R:={(t,x)∈([0,T]×Ω)∩ Q_R(0): |Du^h_k(t,x)|≥ s},


and apply Chebyshev's inequality:


    |A_k,s,R|   ≤1/s^p∫_A_k,s,R|Du^h_k|^px
       ≤c(1+ ||u_0||_L^2(Ω;m)^2 + F^p'_L^p'(Q_T;m))/s^p,


where we have used (<ref>) along with the fact that for r_1≥ r_2, we have |s|^r_2≤ 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→∞ we deduce that (<ref>) holds.


By taking 𝔎=m (see Remark <ref> following the statement of Theorem <ref>), we may apply Theorem <ref>, to deduce the existence of a parameterized family of Young measures (ν_t,x)_(t,x) ∈ Q_T and a further subsequence (which we still index by h) for which


    a(Du^h)⇀⟨ν,a⟩    in L^α(Q_T;m × n),


where 1<α<p/q-1 (see Remark <ref>, combined with the fact that a(A)≤ c(1+|A|^q-1), which follows from (<ref>) by recalling  that q:=max_i=1,…,m p_i, whence a ∈ E_p/α). By comparing
(<ref>) and (<ref>) and noting the uniqueness of weak limits there is then a subsequence of u^h (still indexed by h) such that


    a(Du^h)⇀⟨ν,a⟩    in L^q̂'(Q_T;m).


Here we have used that q̂' = min{q',2} = min{q/q-1,2}≥min{p/q-1,2}≥min{α,2}, with α as above, and noted that (<ref>) implies that a(Du^h)⇀⟨ν ,a⟩ in L^min{α,2}(Q_T;m × n) because α≥min{α,2}.

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


    u ∈ L^∞(0,T;L^2(Ω;m))∩ L^p(0,T;W^1,p_0(Ω;m)),     Du_i  ∈ L^p_i(0,T;L^p_i(Ω;m× n)),   i=1,…,m,


with


    ∂_tu∈ L^q̂'(0,T; W^-1,q̂'(Ω;m)),     u ∈ C_w([0,T];L^2(Ω;m))    ⟨ν ,a⟩∈ L^q̂'(Q_T;m× n)


satisfying


    ∫_0^T⟨∂_tu,ϕ⟩+∫_Ω⟨ν_t,x,a⟩ : Dϕ + Bu·ϕxt=∫_0^T∫_Ω F·ϕxt


for all ϕ∈ L^q̂(0,T; W_0^1,q̂(Ω;m)).

It remains to be shown that the initial condition is satisfied. By Lemma <ref> it follows that we can extract a further subsequence h_l along which we have that


    ∫_Ω u^h_l(t,x)·φ(x)x→∫_Ω u(t,x)·φ(x)x,


uniformly in C([0,T]), for all φ∈ L^2(Ω;m). Let us decompose


    ∫_Ω   (u_0(x)-u(0,x))·φ(x)x
       =∫_Ω(u_0(x)-u^h_l(0,x))·φ(x)x+∫_Ω (u^h_l(0,x)-u^h_l(t,x))·φ(x)x
           +∫_Ω (u^h_l(t,x)-u(t,x))·φ(x)x+∫_Ω (u(t,x)-u(0,x))·φ(x)x
    =: I_l + II_l(t) + III_l(t) + IV(t).


The term I_l converges to zero as h_l→0_+ by the assumed strong convergence in L^2(Ω;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→0_+ by (<ref>). Hence, for any ε>0, there exists an h_l^∗=h_l^∗(ε,φ) ∈ (0,h_0], such that


    |∫_Ω(u_0(x)-u^h_l^∗(0,x))·φ(x)x|<1/4ε    max_t ∈ [0,T]|∫_Ω (u^h_l^∗(t,x)-u(t,x))·φ(x)x| < 1/4ε.


Consequently, for any ε>0, there exists an h_l^∗∈ (0,h_0], such that


    |∫_Ω (u_0(x)-u(0,x))·φ(x)x|
            < 1/2ε +|∫_Ω (u^h_l^∗(0,x)-u^h_l^∗(t,x))·φ(x)x|
    +|∫_Ω (u(t,x)-u(0,x))·φ(x)x|.


With h_l^∗ fixed, we now pass to the limit t→0_+  in this inequality. The term II_l(t) converges to zero as t→0_+ because u^h_l^∗, as a solution to the Galerkin equation (<ref>), is continuous and satisfies the discretized initial condition. The term IV(t) converges to as t → 0_+ since by Lemma <ref> we have that u∈ C_w([0,T];L^2(Ω;m)), the space of weakly continuous functions from the interval [0,T] into L^2(Ω;m).
Thus, by passing to the limit t → 0_+ in the last inequality we deduce that, for any φ∈ L^2(Ω;m) and any ε>0,


    |∫_Ω (u_0(x)-u(0,x))·φ(x)x| < 1/2ε.


Taking φ(x) = u_0(x)-u(0,x) and letting ε→ 0_+ then yields that u_0(x)=u(0,x) for almost every x∈Ω.









§ THE FULLY DISCRETE SCHEME AND ITS CONVERGENCE ANALYSIS



In the course of the proof of Theorem <ref> Peano's Theorem was applied in order to show that one can solve the Galerkin equation (<ref>). 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∈ V^h_m, we inductively define u_i+1^h, i=0,…, N-1, as solutions to the following fully discrete problem.

Let Δ t = T/N where N ≥ 1, and given u_i^h∈ V^h_m, find u_i+1^h ∈ V^h_m such that


    ∫_Ωu_i+1^h-u_i^h/Δ t·ϕ^h+a(Du_i+1^h) : Dϕ^h+Bu_i+1^h·ϕ^hx=∫_Ω F_i+1·ϕ^hx    {[ ,; . ].


Here,


    F_i(x):=1/Δ t∫_t_i-1^t_iF(t,x)t,


and, letting t_i:=iΔ t, we define on Q_T the function


    F̃_Δ t(t,x):=F_i(x)  for t∈(t_i-1,t_i],   i=1,…,N,     x ∈Ω.



Given a function u^h_i, the fact that there exists a function u^h_i+1 satisfying (<ref>) 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 <cit.>). We begin the convergence analysis with the following estimate, which is reminiscent of the energy estimate obtained in the proof of Theorem <ref>, for the semidiscrete problem considered there.




There exists a positive constant c, independent of h and Δ t,
such that the functions u_i^h, i=1,…,N, satisfy, for all h ∈ (0,h_0],


    max_i=1,…,Nu_i^h_L^2(Ω;m)^2+∑_i=1^N∫_Ω |u_i^h-u_i-1^h|^2x + Δ t∑_i=1^N∫_Ω a(Du_i^h) : Du_i^hx
    ≤ c(u_0^h_L^2(Ω;m)^2+TF_L^2(Q_T;m)^2).


We note here that a(ξ) : ξ≥ 0 for all ξ∈m× n, and therefore each of the terms appearing on the left-hand side of the inequality (<ref>) is nonnegative. In addition, the expression on the right-hand side of (<ref>) is further bounded above by a constant that is independent of h and Δ t.


 We begin by choosing ϕ^h=u_i+1^h in (<ref>), and noting the following identity:


    (a-b)· a/Δ t=|a|^2-|b|^2/2Δ t+|a-b|^2/2Δ t    ∀  a,b ∈m.


Using this in (<ref>) gives


    ∫_Ω|u_i+1^h|^2-|u_i^h|^2/2Δ t+|u_i+1^h-u_i^h|^2/2Δ t+a(Du_i+1^h) : Du_i+1^h + Bu_i+1^h· u_i+1^hx=∫_Ω F_i+1· u_i+1^hx.


Now we use that Bv· v≥0 for all v ∈m, multiply by 2Δ t and sum through i=0,…,M for some M<N to get


    ∫_Ω |u_M+1^h|^2x   +∑_i=0^M∫_Ω|u_i+1^h-u_i^h|^2x+2Δ t∑_i=0^M∫_Ω a(Du_i+1^h) : Du_i+1^hx
       ≤∫_Ω |u_0^h|^2x+2∑_i=0^M∫_Ω(∫_t_i^t_i+1F(t,x)t)· u^h_i+1x
       ≤∫_Ω |u_0^h|^2x+ 2TF_L^2(Q_T;m)^2+ 1/2max_i=1,…,Nu_i^h_L^2(Ω;m)^2 .


By omitting the second and the third term from the left-hand side of the inequality (<ref>) noting the independence of its right-hand side of M, and then taking the maximum over M=1,…,N-1, yields the desired bound on the first term on the left-hand side of the inequality (<ref>). We then return with that bound to the inequality (<ref>) to further estimate its right-hand side from above, whilst omitting the first term from the left-hand side of (<ref>);
thus we arrive at the desired bound on the second and third term on  the left-hand side of the inequality (<ref>). The expression on the right-hand side of (<ref>) is clearly independent of Δ 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(Ω;m).





As the functions u^h_i, i=0,…, N, are defined only on {t_0=0,t_1, …, t_N=T}×Ω rather than on the whole of Q_T, we shall next extend them to Q_T. Thus we define the following two functions:


    u_Δ t^h(t):=t-t_i-1/Δ tu_i^h+t_i-t/Δ tu_i-1^h   for t∈[t_i-1,t_i],     i=1,…,N,


which is continuous and piecewise linear with respect to t ∈ [0,T], and


    ũ_Δ t^h(t):=u_i^h   for t∈ (t_i-1,t_i],   i = 1,…, N,    and  ũ_Δ t^h(t):=u_0^h   for t ∈ [-Δ t,0],


which is piecewise constant with respect to t ∈ [-Δ t,T]. With these new notations, (<ref>) can be written as


    ∫_Ω∂_tu_Δ t^h·ϕ^h+a(Dũ_Δ t^h) : Dϕ^h+Bũ_Δ t^h·ϕ^hx=∫_ΩF̃_Δ t·ϕ^hx   ,


for all  ϕ^h∈ V^h_m. We aim now to obtain uniform bounds on the sequences {u_Δ t^h} and {ũ_Δ t^h}, which will allow us to complete our passage to the limit in Δ t, possibly along some subsequence.




There exists a positive constant c, independent of h and Δ t, such that, for all h ∈ (0,h_0],


    ũ^h_Δ t^2_L^∞(0,T;L^2(Ω;m))≤ c (u_0^h_L^2(Ω;m)^2+F_L^2(Q_T;m)^2),



    ũ^h_Δ t-ũ^h_Δ t(·-Δ t)^2_L^2(0,T;L^2(Ω;m))≤ c (u_0^h_L^2(Ω;m)^2+F_L^2(Q_T;m)^2),



    a(Dũ^h_Δ t)_L^q̂'(Q_T;m× n)^q̂'≤ c(1+u_0^h_L^2(Ω;m)^2+F_L^2(Q_T;m)^2 +F_L^p'(Q_T;m)^p'),  ,



    ∑_i=1^m D(ũ^h_Δ t)_i_L^p_i(0,T;L^p_i(Ω;n))^p_i≤ c(1+u_0^h_L^2(Ω;m)^2+F_L^p'(Q_T;m)^p').


Furthermore, the right-hand sides of these inequalities are bounded above by a positive constant, independent of h ∈ (0,h_0] and
Δ t.





The first two bounds follow directly from the bounds on the first two terms in (<ref>), 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 (<ref>). 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(ξ):ξ≥ 0 for all ξ∈m × n, we have that


    |a(Du_i^h)|^q' = |a(Du_i^h) : a(Du_i^h)|^q'/2 = |(a(Du_i^h) : Du_i^h )  K(Du_i^h)|^q'/2 = (a(Du_i^h) : Dũ^h_Δ t)^q'/2  |K(Du_i^h)|^q'/2.



First suppose that q>2; then, q̂=q and q̂'=q' ∈ (1,2); hence, by Young's inequality,


    |a(Du_i^h)|^q̂' = |a(Du_i^h)|^q'≤q'/2 a(Du_i^h) : Du_i^h + 2-q'/2 |K(Du_i^h)|^q'/2-q'.


Therefore, by (<ref>) and (<ref>), we have that


    ∫_0^T∫_Ω |a(Dũ^h_Δ t)|^q̂'xt = ∫_0^T∫_Ω |a(Dũ^h_Δ t)|^q'xt=Δ t∑_i=1^N∫_Ω |a(Du_i^h)|^q'x
           ≤q'/2Δ t∑_i=1^N∫_Ω a(Du_i^h) : Du_i^hx +2-q'/2Δ t∑_i=1^N∫_Ω |K(Du_i^h)|^q'/2-q'x
           ≤ c(q)(u_0^h_L^2(Ω;m)^2+T F_L^2(Ω;m)^2) + c(q) Δ t∑_i=1^N∫_Ω(c_1∑_j=1^m(|μ_j|^2 + |Du_i,j^h|^2)^p_j-2/2)^q'/2-q'x
           ≤ c(q) (u_0^h_L^2(Ω;m)^2+TF_L^2(Q_T;m)^2)
    +c(c_1,q,m)Δ t∑_i=1^N∫_Ω[∑_j=1^m |μ_j|^(p_j-2)q/q-2 + ∑_p_j ≥ 2^m |Du_i,j^h|^(p_j-2)q/q-2] x.


As q:=max{p_1,…,p_m} >2, for any j ∈{1,…,m} such that p_j ≥ 2 we have that  (p_j-2)q/q-2<p_j; it thus follows that


    ∫_0^T∫_Ω |a(Dũ^h_Δ t)|^q̂'xt≤ c(q)(u_0^h_L^2(Ω;m)^2+TF_L^2(Q_T;m)^2)
             +c(c_1, q,m)   T   |Ω| ∑_j=1^m |μ_j|^(p_j-2)q/q-2
    + c(q,m,c_1,p_1,…,p_m,|Ω|) Δ t∑_i=1^N ∫_Ω(∑_j=1^m|Du_i,j^h|^p_j) x.

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


    Δ t∑_i=1^N ∫_Ω(∑_j=1^m|Du_i,j^h|^p_j dx)≤ c(1+u_0^h_L^2(Ω;m)^2+F_L^p'(Q_T;m)^p').


Substituting this into the right-hand side of the previous inequality completes the proof of the lemma in the case when
q̂ = q>2.


If q̂=2 (and therefore q̂'=2), then p_j ∈ (1,2] for all j=1,…,m. It then follows from (<ref>) that


    0 ≤ K(A) ≤ c_1∑_i=1^m |μ_i|^p_i-2    ∀  A ∈ℝ^m × n,


whereby


    |a(Du_i^h)|^q̂' = |a(Du_i^h)|^2= (a(Du_i^h) : Dũ^h_Δ t)  K(Du_i^h) ≤
    c_1(∑_i=1^m |μ_i|^p_i-2) (a(Du_i^h) : Dũ^h_Δ t).


Hence, by Lemma <ref>,


    ∫_0^T∫_Ω |a(Dũ^h_Δ t)|^q̂'xt=Δ t∑_i=1^N∫_Ω |a(Du_i^h)|^q'x
           ≤ c_1(∑_i=1^m |μ_i|^p_i-2) Δ t∑_i=1^N∫_Ω a(Du_i^h) : Du_i^hx
           ≤ c(c_1,p_1,…,p_m, μ_1,…, μ_m, m, T)
    (u_0^h_L^2(Ω;m)^2+F_L^2(Ω;m)^2).


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 → 0_+. That completes the proof of the lemma.



Next, we prove the required uniform bounds on the sequence {u^h_Δ t}.



There exists a positive constant c, independent of h and Δ t, such that,  for all h ∈ (0,h_0],


    u^h_Δ t^2_L^∞(0,T;L^2(Ω;m))≤ c (u_0^h_L^2(Ω;m)^2+F_L^2(Q_T;m)^2);


furthermore, with  q̂:=max{q,2},


    ∂_t u^h_Δ t_L^q̂'(0,T;W^-1,q̂'(Ω;m))≤ c


and


    (Δ t)^-1ũ^h_Δ t-ũ^h_Δ t(·-Δ t)_L^q̂'(0,T;W^-1,q̂'(Ω;m))≤ c.





The first inequality is just a restatement of the bound on the first term in the inequality (<ref>) in terms of
u^h_Δ t, and the third inequality is a restatement of the second inequality in terms of ũ^h_Δ t. We shall therefore concentrate on the proof of the second stated inequality. The argument proceeds along the same lines as the proof of (<ref>). For each φ∈ L^q̂(0,T;W^1,q̂_0(Ω;m)) we have, using the stability in W^1,q̂_0(Ω;m) of the L^2(Ω;m) orthogonal projector P^h:L^2(Ω;m)→ V^h_m (cf. (<ref>)), that


    ∫_0^T (∂_t u^h_Δ t, φ)  t    = ∫_0^T (∂_t u^h_Δ t, P^hφ)  t
     = ∑_i=0^N-1∫_t_i^t_i+1∫_Ωu_i+1^h - u_i^h/Δ t· P^hφ x  t
       = -  ∑_i=0^N-1∫_t_i^t_i+1∫_Ω a(Du^h_i+1) · DP^hφ x  t +
    ∑_i=0^N-1∫_t_i^t_i+1 (F_i+1 - Bu_i+1^h)· P^hφ x  t
       = -  ∫_0^T ( a(Dũ^h_Δ t(t)),DP^hφ)  t + ∫_0^T ( F̃_Δ t - Bũ^h_Δ t ,P^hφ )  t
       ≤a(Dũ^h_Δ t)_L^q̂'(Q_T;m × n)Dφ_L^q̂(Q_T;m× n)
            + (F̃_Δ t_L^2(Q_T;m) + |B| ũ^h_Δ t_L^2(Q_T;m)) φ_L^2(Q_T;m)
       ≤ c (Dφ_L^q̂(Q_T;m× n) + φ_L^q̂(Q_T;m))  ≤  c φ_L^q̂(0,T;W^1,q̂(Ω;m)),


where in the transition to the last line we have used the third and the first bound from Lemma <ref>, together with the assumed strong convergence of u^h_0 to u_0 in L^2(Ω;m) to deduce that the constant c is independent of h and Δ t. Dividing through by φ_L^q̂(0,T;W^1,q̂(Ω;m)), taking the supremum over all φ∈ L^q̂(0,T;W^1,q̂_0(Ω;m)) and recalling the definition of the norm of the dual space L^q̂'(0,T;W^-1,q̂'(Ω;m))=(L^q̂(0,T;W^1,q̂_0(Ω;m)))' the stated inequality directly follows.



We shall now discuss the convergence of the sequence of functions {u_Δ t^h} with h ∈ (0,h_0] fixed.



There exists a subsequence of Δ t, labelled Δ t_k, and a function u^h∈ C([0,T];V^h_m) such that u_Δ t_k^h→ u^h in C([0,T];L^2(Ω;m)), and ∂_tu_Δ t_k^h⇀∂_tu^h in L^q̂'(0,T;W^-1,q̂'(Ω;m)) for any (fixed) h ∈ (0,h_0], as k →∞.





We begin by noting that, for any Δ t and any h ∈ (0,h_0], the functions u^h_Δ t are continuous in both t and x. Thanks to the first two bounds in Lemma <ref> there exists a subsequence of Δ t, labelled Δ t_k, and a function u^h∈ L^∞(0,T;V^h_m) such that u_Δ t_k^h ∗⇀ u^h in L^∞(0,T;L^2(Ω;m)), and ∂_tu_Δ t_k^h⇀∂_tu^h in L^q̂'(0,T;W^-1,q̂'(Ω;m)). By Lemma <ref>, u^h∈ C_w([0,T];L^2(Ω;m)) and u_Δ t_k^h → u^h in C_w([0,T];L^2(Ω;m)). As both u^h_Δ t(t,·) and u^h(t,·) belong to the finite-dimensional space V^h_m for all t ∈ [0,T], it follows that u^h∈ C([0,T];L^2(Ω;m)) and u_Δ t_k^h → u^h in u^h∈ C([0,T];L^2(Ω;m)), because u^h∈ 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 →∞ for any (fixed) h ∈ (0,h_0].



Next, we study the convergence of the sequence of functions {ũ_Δ t^h} with h ∈ (0,h_0] fixed.


There exists a function ũ^h ∈ L^∞(0,T;V^h_m) and a subsequence Δ t_k such that


    ũ_Δ t_k^h ⇀ũ^h
         _j=1^m L^p_j(0,T;W^1,p_j_0(Ω;m))
             ũ_Δ t_k^h →ũ^h
        L^r(0,T;L^2(Ω;m)),


for any r ∈ [1,∞) and fixed h ∈ (0,h_0], as k→∞. Furthermore, along this subsequence we have a(Dũ_Δ t_k^h)→ a(Dũ^h) in L^s(Q_T;m× n) for any s ∈ [1,q̂') and a(Dũ_Δ t_k^h)⇀ a(Dũ^h) in L^q̂'(Q_T;m× n), for any fixed h ∈ (0,h_0], as k→∞.




 Thanks to the first and fourth bound from Lemma <ref>, there exists a
ũ^h ∈ L^∞(0,T;V^h_m) and a subsequence Δ t_k such that ũ_Δ t_k^h ∗⇀ũ^h in L^∞(0,T;L^2(Ω;m)) and


    ũ_Δ t_k^h ⇀ũ^h
           _j=1^m L^p_j(0,T;W^1,p_j_0(Ω;m)),


for any fixed h ∈ (0,h_0], as k→∞.

By applying Theorem 1 from <cit.> (with X=W^1,p_0(Ω;m), B=L^2(Ω;m) and Y=W^-1,q̂'(Ω;m) and r=1, there), we deduce that {ũ^h_Δ t_k} is relatively compact in L^p(0,T;L^2(Ω;m)), for each fixed h ∈ (0,h_0]; therefore, because it is a bounded sequence in L^∞(0,T;L^2(Ω;m)), by interpolation  {ũ^h_Δ t_k} is also relatively compact in L^r(0,T;L^2(Ω;m)), for each fixed h ∈ (0,h_0] and all r ∈ [1,∞).

Hence, in particular, {ũ^h_Δ t_k} is relatively compact in L^1(0,T;L^1(Ω;m)), for each fixed h ∈ (0,h_0]. As ũ^h_Δ t_k(t,·) ∈ V^h_m for all t ∈ [0,T], and V^h_m is finite-dimensional, for each fixed h ∈ (0,h_0], by norm equivalence in finite-dimensional spaces, it follows that {ũ^h_Δ t_k} is relatively compact in L^1(0,T;W^1,1_0(Ω;m)), for each fixed h ∈ (0,h_0]. Therefore Dũ^h_Δ t_k→ Dũ^h in L^1(Q_T;m × n) as k →∞, for each fixed h ∈ (0,h_0]. Thus we can extract a subsequence with respect to k (not indicated) such that Dũ^h_Δ t_k→ Dũ^h a.e. on Q_T as k →∞, for each fixed h ∈ (0,h_0]. Thanks to the continuity of a, we then have that a(Dũ^h_Δ t_k) → a(Dũ^h) a.e. on Q_T as k →∞, for each fixed h ∈ (0,h_0]. Since by the third inequality from Lemma <ref> the sequence {a(Dũ^h_Δ t_k)} is weakly compact in L^q̂'(Q_T;m × n) and hence in particular also in L^1(Q_T;m × n), for each fixed h ∈ (0,h_0], it follows by Vitali's convergence theorem that a(Dũ^h_Δ t_k) → a(Dũ^h) strongly in L^1(Q_T;m × n) as k →∞, for each fixed h ∈ (0,h_0], and therefore also strongly in  L^s(Q_T;m × n) as k →∞, for each fixed h ∈ (0,h_0] and any s ∈ [1,q̂'). By the uniqueness of the weak limit it then follows that a(Dũ_Δ t_k^h)⇀ a(Dũ^h) weakly in L^q̂'(Q_T;m× n), for any fixed h ∈ (0,h_0], as k→∞.



We now show that the limit functions u^h and ũ^h defined above are equal when considered as elements of the space L^2(Q_T;m).


 The limiting functions u^h and ũ^h are equal in the space L^2(Q_T;m)=L^2(0,T;L^2(Ω;m)).




 We apply the triangle inequality to see that, with Δ t_k = T/N_k and h ∈ (0,h_0],


    u^h-ũ^h_L^2(0,T; L^2(Ω;m))   ≤u^h-u_Δ t_k^h_L^2(0,T;L^2(Ω;m))+u_Δ t_k^h-ũ_Δ t_k^h_L^2(0,T;L^2(Ω;m))
             +ũ_Δ t_k^h-ũ^h_L^2(0,T;L^2(Ω;m)).


The first and last terms converge to zero in the limit as Δ t_k→0_+  thanks to the convergence results stated in Lemma <ref> and Lemma <ref> (with r=2), respectively, for each fixed h ∈ (0,h_0]. We now show that the second term also converges to zero in the limit:


    u_Δ t_k^h-ũ_Δ t_k^h_L^2(0,T;L^2(Ω;m))^2   =∫_0^T∫_Ω|u_Δ t_k^h-ũ_Δ t_k^h|^2xt
       =∑_i=1^N_k∫_Ω∫_t_i-1^t_i|t-t_i-1-Δ t_k/Δ t_ku_i^h+t_i-t/Δ t_ku_i-1^h|^2tx
       =Δ t_k/3∑_i=1^N_k∫_Ω |u_i^h-u_i-1^h|^2x
       =Δ t_k/3ũ^h_Δ t_k-ũ^h_Δ t_k(·-Δ t)^2_L^2(0,T;L^2(Ω;m)).


Thanks to the second inequality in Lemma <ref>, the right-hand side of the last inequality converges to zero in the limit of Δ t_k → 0_+ (N_k →∞ as k →∞). Hence, u^h=ũ^h, as has been asserted.



Next, we discuss the attainment of the initial condition.


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


    lim_t→0_+u^h(t,·)-u_0^h(·)_L^2(Ω;m)=0,


for each fixed h ∈ (0,h_0].



Thanks to Lemma <ref>, for any ϵ>0 there is an ℓ∈ℕ such that u^h(t,·)-u_Δ t_k_l^h(t,·)_L^2(Ω;m)<ϵ for all l≥ℓ and all t ∈ [0,T]. Furthermore, by the triangle inequality:


    u^h(t,·)-u_0^h(·)_L^2(Ω;m)≤u^h(t,·)-u_Δ t_k_l^h(t,·)_L^2(Ω;m)+u_Δ t_k_l^h(t,·)-u_0^h(·)_L^2(Ω;m)


for any l ∈ℕ. Hence, in particular, for l = ℓ (fixed),


    u^h(t,·)-u_0^h(·)_L^2(Ω;m)≤u^h(t,·)-u_Δ t_k_ℓ^h(t,·)_L^2(Ω;m) + ϵ.


As u_Δ t_k_ℓ^h is a continuous function in t, passage to the limit in this inequality, using that lim_t→0_+u_Δ t_k_l^h(t,x)=u_0^h(x), implies that, for each ϵ>0,


    0 ≤lim sup_t → 0_+u^h(t,·)-u_0^h(·)_L^2(Ω;m)≤ϵ.


Thus we have the asserted result.




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


 We have that lim_Δ t→0_+F̃_Δ t(·,x)-F(·,x)_L^r(0,T) for all r ∈ [1,∞) and a.e. x ∈Ω.





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




We are now ready for passage to the limit in (<ref>) as Δ t → 0_+. By Lemma <ref>, Lemma <ref>, and Lemma <ref>, we see that the limit function u^h satisfies


    ∫_0^T∫_Ω∂_tu^h·Φ^h+a(Du^h) : DΦ^h+Bu^h·Φ^hxt=∫_0^T∫_Ω F·Φ^hxt     ∀ Φ^h ∈ V^h_m,   ∀ h ∈ (0,h_0],


with u^h(0,x) = u^h_0(x) for x ∈Ω. We can now pass to the limit h → 0_+ by (the proof of) Theorem <ref>, to deduce that as the discretization parameters converge to zero the solution of the fully discrete scheme (<ref>) converges to a Young measure solution of the problem (<ref>)–(<ref>) under consideration.

In preparation for the considerations in the next section, we now discuss continuous dependence of u^h_Δ t on the initial data. We start by considering two functions, u_0^h ∈ V^h_m and v_0^h∈ V^h_m, and we let u_1^h and v_1^h be respective solutions to (<ref>) 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 ϕ^h=u_1^h-v_1^h, it follows that


    ∫_Ω |u_1^h-v_1^h|^2 +Δ t(a(Du_1^h)-a(Dv_1^h))·(Du_1^h-Dv_1^h)x ≤∫_Ω(u_0^h-v_0^h)·(u_1^h-v_1^h)x,


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


    u_1^h-v_1^h_L^2(Ω;m)^2+2Δ t∫_Ω(a(Du_1^h)-a(Dv_1^h))·(Du_1^h-Dv_1^h)x≤u_0^h-v_0^h_L^2(Ω;m)^2.


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 <ref> that Du_1^h_L^∞(Ω;m × n) and Du_1^h_L^∞(Ω;m × n) are both bounded by a constant C=C(h,u_0_L^2(Ω;m),F_L^2(Q_T;m)). Therefore we can make use of the local Lipschitz condition satisfied by a and the inverse inequality


    Du_1^h-Dv_1^h_L^2(Ω;m)≤ C(h) u_1^h-v_1^h_L^2(Ω;m),


(with C(h) = Ch^-1 when the family of triangulations 𝒯^h is quasiuniform) to deduce, by suppressing in our notation for the constant C(h,a_ Lip,loc) below its dependence on the (fixed) data u_0_L^2(Ω;m) and F_L^2(Q_T;m), that


    u_1^h-v_1^h_L^2(Ω;m)^2(1-C(h,a_ Lip,loc)Δ t)≤u_0^h-v_0^h_L^2(Ω;m)^2.


Therefore by choosing Δ t sufficiently small for h ∈ (0,h_0] fixed, so that


    0<1-C(h,a_ Lip,loc)Δ t<1,


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_ Lip,loc) remains the same regardless of the choice of the time level i ∈{1,…, N}, we obtain that, for each i=1,…,N, the following estimate holds: 


    u^h_i-v^h_i_L^2(Ω;m)^2(1-C(h,a_ Lip,loc)Δ t)^i≤u_0^h-v_0^h_L^2(Ω;m)^2.


Recall that Δ t and N are related via the identity T=NΔ t. Hence,


    (1-C(h,a_ Lip,loc)Δ t)^N=(1-TC(h,a_ Lip,loc)/N)^N → e^-TC(h, a_ Lip,loc)    .


To transfer this continuous dependence on the initial data to the limiting functions u^h and v^h as Δ t → 0_+, we consider u^h_Δ t and v^h_Δ t, defined, for i=1,…,N, by


    u^h_Δ t(t,x)=t-t_i-1/Δ tu^h_i(x)+t_i-t/Δ tu^h_i-1(x)  for t∈[t_i-1,t_i] and x ∈Ω,


and, analogously,


    v^h_Δ t(t,x)=t-t_i-1/Δ tv^h_i(x)+t_i-t/Δ tv^h_i-1(x)  for t∈[t_i-1,t_i]  and x ∈Ω.


It then follows from (<ref>) that


    sup_t∈[0,T]u^h_Δ t-v^h_Δ t_L^2(Ω;m)^2(1-C(h,a_ Lip,loc)Δ t)^N
    ≤u_0^h-v_0^h_L^2(Ω;m)^2.


Let Δ t_k be a sequence such that u^h_Δ t_k→ u^h and v^h_Δ t_k→ v^h in C([0,T];L^2(Ω;m)) and k→∞ (cf. Lemma <ref>). Now passing k→∞ (whereby Δ t_k = T/N_k → 0_+ and N_k →∞) we see that


    u^h-v^h_C([0,T];L^2(Ω;m))≤ e^1/2TC(h,a_ Lip,loc)u_0^h-v_0^h_L^2(Ω;m).


We have therefore shown the following result.



Let all of the assumptions of Theorem <ref> hold, with the additional assumption that, for h ∈ (0,h_0] fixed, Δ t = T/N is chosen so that 0 < 1 - C(h,a_ Lip,loc)Δ t < 1. Then, there is a unique solution to (<ref>), which depends continuously on the choice of discretized initial data, and which can be constructed numerically by the approximations considered in this subsection.



The above proposition in particular implies that the map S^h_t sending the initial condition u_0^h to S^h_tu_0^h=u^h(t,·) is (Lipschitz) continuous from V^h_m to V^h_m for all t ∈ [0,T] and h ∈ (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,·). Recall that we are writing the functions u^h in terms of a finite element basis


    u^h(t,x)=∑_i=1^𝒩(h)α_i^h(t)ϕ_i^h(x).


The above considerations enable us to view the map from 𝒩(h) to 𝒩(h), which sends the vector α^h(0) of coefficients of the discretized initial condition in terms of the finite-element basis to the vector α^h(t) of coefficients of u^h, in the same finite-element basis, as a continuous function. To see this, let α^h and β^h be the coefficients of u^h and v^h respectively, and let


    M_i,j^h:= (ϕ^h_i,ϕ^h_j)/ϕ^h_i_L^2(Ω;m)ϕ^h_j_L^2(Ω;m),      i,j=1,…,𝒩(h).


The (normalized) Gramm matrix M^h=(M_i,j^h)_i,j=1^𝒩(h) is symmetric positive definite, with
C_0 I ≤ M^h ≤  C_1 I, where 0<C_0 ≤ C_1, in the sense of symmetric positive definite matrices. By expressing


    u^h(t,x)-v^h(t,x)=∑_i=1^𝒩(h)(α^h_i(t)-β^h_i(t))ϕ^h_i(x),


it follows that,


    C_0min_i=1,…,𝒩(h)ϕ^h_i_L^2(Ω;m)max_t∈[0,T]   (∑_i=1^𝒩(h)|α^h_i(t)-β^h_i(t)|^2)^1/2
       ≤ C_0max_t∈[0,T](∑_i=1^𝒩(h)|α^h_i(t)-β^h_i(t)|^2∫_Ω |ϕ^h_i(x)|^2x)^1/2
       ≤max_t∈[0,T](∑_i,j=1^𝒩(h)(α^h_i(t)-β^h_i(t))(α^h_j(t)-β^h_j(t))∫_Ωϕ^h_i(x) ·ϕ^h_j(x)x)^1/2
       =max_t∈[0,T](∫_Ω∑_i,j=1^𝒩(h)(α^h_i(t)-β^h_i(t))ϕ^h_i(x)(α^h_j(t)-β^h_j(t))ϕ^h_j(x) x)^1/2


    =max_t∈[0,T](∫_Ω|∑_i=1^𝒩(h)(α^h_i(t)-β^h_i(t) ϕ^h_i(x)|^2x)^1/2
        = u^h-v^h_C([0,T];L^2(Ω;m))≤ C(T,h,a_ Lip,loc)u_0^h-v_0^h_L^2(Ω;m)
       ≤ C(T,h,a_ Lip,loc)  C_1 max_i=1,…,𝒩(h)ϕ^h_i_L^2(Ω;m)(∑_i=1^𝒩(h)|α^h_i(0)-β^h_i(0)|^2)^1/2.


Thus the map α^h(0) ∈𝒩(h)↦α^h(t)∈𝒩(h) is (Lipschitz) continuous for all t∈[0,T] and Δ t sufficiently small so as to ensure that 0 < 1 - C(h,a_ Lip,loc)Δ t < 1, with h ∈ (0,h_0] fixed.





§ COMPUTATION OF YOUNG MEASURE SOLUTIONS


Theorem <ref> and the results of Section <ref> (in particular, Proposition <ref>) give criteria under which Young measure solutions exist as limits of solutions to corresponding semi- or fully-discrete problems. However, the results of Section <ref> and Section <ref> only give “indirect" information about the Young measure ν through its action on functions belonging to the space E_p/α for α∈(1,p/q-1). Indeed, note that for any function b which has the form (<ref>)–(<ref>), we can copy the calculations used to deduce (<ref>) to see that b(Du^h)⇀⟨ν,b⟩ in L^q̂'(Q_T;m). In this section we shall discuss an algorithm for the numerical approximation of Young measure solutions to systems of the form (<ref>), which allows for a more direct approximation of the Young measure ν. The ideas presented here have been inspired by <cit.> and <cit.>, 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 (<ref>) 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 ν 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.








 §.§ Preliminary definitions and results


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



We let 𝒫(m) denote the set of all probability measures on m, and define, for r>0, the subset 𝒫^r(m) of all probability measures μ∈𝒫(m) for which ⟨μ,|ξ|^r⟩<∞, where ξ denotes the identity function.



We denote by (𝔛,ℱ,P) a probability space, with ℱ being a σ-algebra of sets on the space 𝔛, and P being a probability measure, and we let u:𝔛× Q_T→m be a random field (that is, a jointly measurable function). Given this, we define the law of u as follows:


    μ_t,x(B):=P({ω∈𝔛: u(w;t,x)∈ B}),


for Borel subsets B⊂m. It is clear that μ_t,x(m)=1 for all (t,x) ∈ Q_T. We have the following result.



Let (ω;t,x) ∈𝔛× Q_T ↦ u(ω;t,x) ∈m be a random field; then, (<ref>) is equivalent to


    ⟨μ_t,x,g⟩=∫_𝔛 g(u(ω;t,x))P(ω)    ,


for every continuous r-component vector function g, defined on m, which is such that ∫_𝔛 |g(u(ω;t,x))|P(ω) is finite for a.e. (t,x) ∈ Q_T.




 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 m into 
such that ∫_𝔛 |g(u(ω;t,x))|P(ω) is finite for a.e. (t,x) ∈ Q_T. Let us first consider the case when the function g:m→ is nonnegative, and take a sequence of simple functions g_n defined on m by


    g_n(ξ)=∑_k=1^n2^n+1 b^n_kχ_A^n_k(ξ),    ξ∈m,   n=1,2,…,


where A^n_k, k=1,…, n, are disjoint measurable subsets of m, n is a positive integer, and b^n_k ∈, k=1,…,n, are such that g_n(ξ) increases to g(ξ). This can be achieved by defining


    A^n_k:= g^-1([k-1/2^n,k/2^n)),      1 ≤ k ≤ n 2^n,


and


    A_n,n2^n+1:=g^-1([n,∞)).


Observe that g_n(ξ) ≤ g(ξ) for all ξ∈m. Further, for a given ξ∈m we have that g_n(ξ) = n for g(ξ)≥ n whereas for g(ξ)<n we have that


    g_n(ξ)= 1/2^n⌊ 2^n g(ξ)⌋,


where ⌊ y ⌋ denotes the greatest integer ≤ y. As ⌊ 2^n+1 g(ξ)⌋≥ 2 ⌊ 2^n g(ξ)⌋, it follows that g_n+1(ξ) ≥ g_n(ξ) for all ξ∈m. Furthermore, as ⌊ 2^n g(ξ)⌋≥ 2^n g(ξ) - 1, we have that g_n(ξ) ≥ g(ξ) - 2^-n as soon as n>g(ξ). Since g(ξ) ≥ g_n(ξ) it follows that lim_n →∞ g_n(ξ) = g(ξ) for all ξ∈m.



Now, let ℱ^n_k:={ω∈𝔛: u(ω;t,x)∈ A^n_k}. We observe that


    χ_ℱ^n_k(ω)
    =
    1    
    
    0     = 
    1    
    
    0     = χ_A^n_k(u(ω;t,x)).


Hence we have that


    ⟨μ_t,x,g⟩   = ∫_m g(ξ)μ_t,x(ξ) =lim_n→∞∫_mg_n(ξ)μ_t,x(ξ) =lim_n→∞∑_k=1^n2^n+1b^n_kμ_t,x(A^n_k)
       =lim_n→∞∑_k=1^n2^n+1 b^n_kP(ℱ^n_k)
    =lim_n→∞∫_𝔛∑_k=1^n2^n+1b^n_kχ_ℱ^n_k(ω)P(ω)
       =lim_n→∞∫_𝔛∑_k=1^n2^n+1b^n_kχ_A^n_k(u(ω;t,x))P(ω)
       =lim_n→∞∫_𝔛 g_n(u(ω;t,x))P(ω) =∫_𝔛 g(u(ω;t,x))P(ω),


where we have applied the Monotone Convergence Theorem in the first and last lines to exchange the limit with the integral over m. From here the extension to any, not necessarily nonnegative, function g:m→ such that ∫_𝔛 |g(u(ω;t,x))|P(ω) is finite for a.e. (t,x) ∈ Q_T, is straightforward: we decompose g as g=g^+ - g^-, where g^+:=1/2(|g| + g) and g^-:=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 m.




 
It is easy to extend Lemma <ref> to functions of the form g(t,x,u(ω,t,x)) where g is a Carathéodory function, i.e., it is continuous in its third variable for a.e. (t,x) ∈ Q_T and measurable in (t,x)∈ Q_T for each value of the third variable, and ∫_𝔛 |g(t,x,u(ω,t,x))| P(ω) < ∞ for a.e. (t,x) ∈ Q_T.


We have the following proposition, which shows that the law μ=μ_t,x defined by (<ref>) is in fact a Young measure. This is Proposition 1 in Appendix 1 of <cit.>.



If u:𝔛× Q_T→m is jointly measurable, then the law μ=μ_t,x (cf. (<ref>)) defines a Young measure.





The proof is the same as that in <cit.>.



Finally we consider the following result adapted from <cit.>, found therein as Proposition 2 in Appendix 1.



For every Young measure {ν_t,x}_(t,x)∈ Q_T∈𝒫^r, r>0, on m there exists a probability space (𝔛,F,P) with P∈𝒫^r and a Borel measurable function u:𝔛×Ω→m such that u has law given by ν. In particular we may choose the probability space to be the Lebesgue measure on the interval [0,1) with the Borel σ-algebra.





The proof of this result relies on using the characterization given in Lemma <ref>. In <cit.> a proof is supplied with g∈ C_0(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 <cit.> from g∈ C_0() to g∈ C() with ∫_𝔛 |g(u(ω;t,x))|P(ω)<∞ for a.e. (t,x) ∈ Q_T, let φ∈ C^∞_0() be such that 0 ≤φ≤ 1, φ(ξ) ≡ 1 for |ξ| ≤ 1 and φ(ξ) ≡ 0 for |ξ|≥ 2. Let g_k(x):= g(x) φ(x/k). Clearly, g_k is a sequence of functions in C_0(), g_k(ξ)=g(ξ) if |ξ|≤ k, and g_k converges pointwise to g on  as k →∞. Furthermore, |g_k(x)| ≤ |g(x)| for all x ∈. Thus, thanks to the Dominated Convergence Theorem, we have that


    ∫_𝔛 g(u(ω;x))P(ω)   =lim_k→∞∫_𝔛 g_k(u(ω;x))P(ω)
    =lim_k→∞∫_g_k(ξ)ν_x(ξ)
    =∫_g(ξ)ν_x(ξ).


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 <cit.>.



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



Let (𝔛,F,P) be a probability space and let f:𝔛→m be a measurable function. If two random variables, Y_1 and Y_2, defined on 𝔛 are independent and identically distributed, then the random variables f(Y_1) and f(Y_2) are independent and identically distributed on 𝔛.





We first prove that the random variables f(Y_1) and f(Y_2) are identically distributed. We let 𝒜⊂m be a measurable set. It then follows from the measurability of f that


    P(f(Y_1(ω))∈𝒜)   =P(Y_1(ω)∈ f^-1(𝒜))
    =P(Y_2(ω) ∈ f^-1(𝒜))
    =P(f(Y_2(ω))∈𝒜)    ∀ ω∈𝔛.


To show independence, we let 𝒜_1 and 𝒜_2 be measurable sets and calculate


    P((f(Y_1(ω))∈𝒜_1)∩(f(Y_2(ω))∈𝒜_2))   =P((Y_1(ω)∈ f^-1(𝒜_1))∩(Y_2(ω)∈ f^-1(𝒜_2)))
       =P(Y_1(ω)∈ f^-1(𝒜_1)) P(Y_2(ω)∈ f^-1(𝒜_2))
       =P(f(Y_1(ω))∈𝒜_1) P(f(Y_2(ω))∈𝒜_2)    ∀ ω∈𝔛.


That completes the proof of the lemma.





 §.§ 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 <cit.> and <cit.> for the approximation of measure-valued solutions to hyperbolic problems. To this end, let Y(Ω,m) denote the set of all Young measures from Ω to m.



Algorithm A: Let the initial data for an underlying time-dependent PDE be given as a Young measure σ∈ Y(Ω;m) and let Δ denote a PDE discretization parameter associated with a certain numerical scheme (in <cit.> and <cit.> a finite difference scheme is utilized, and so Δ would correspond to a vector of grid sizes in the various co-ordinate directions).



Step 1: Let u_0: ω∈𝔛↦ u_0(ω;·) ∈ L^p(Ω;m) be a random field on a probability space (𝔛,F,P) with law σ_x, meaning that σ_x(E)=P(u_0(ω;x)∈ E) for all Borel sets E ⊂m and x ∈Ω.



Step 2: Evolve the initial random field by applying a suitable numerical scheme, with solution map S_t^Δ, to the initial data u_0(ω;·) for every ω∈𝔛, obtaining an approximate random field u^Δ(ω;·,t):=S_t^Δu_0(ω;·), t ∈ (0,T].



Step 3: Define the approximate measure-valued solution μ^Δ as the law of u^Δ with respect to P, that is, for all Borel sets E⊂m and (t,x) ∈ Q_T,


    μ_t,x^Δ(E)=P(u^Δ (ω;t,x)∈ E).




Note that we should expect a slight difference compared with <cit.> and <cit.> 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 (ω;·,·) ↦ u^Δ(ω;·,·), as performing these computations for every ω∈𝔛 would be infeasible. This is Algorithm 4.3 in <cit.>.



Algorithm B: Let, as above, Δ denote a discretization parameter, and let M∈ℕ. Let σ^Δ be the initial Young measure.



Step 1: From some probability space (𝔛,F,P) draw M independent and identically distributed random fields u_0^Δ,1,…,u_0^Δ,M all with the same law σ^Δ.



Step 2: For each k ∈{1,…,M} and for a fixed ω∈𝔛 approximate the solution to the PDE using the solution operator with initial data u_0^Δ,k(ω); denote u^Δ,k(ω;·,t):=S_t^Δu_0^Δ,k(ω;·).



Step 3: Define the approximate measure-valued solution by


    μ_t,x^Δ,M:=1/M∑_k=1^Mδ_u^Δ,k(·;t,x).




Results from <cit.> (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.




 §.§ 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 (<ref>)–(<ref>).


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 ∈ L^2(Ω;m), which we view as the atomic Young measure
δ_u_0=δ_u_0(x), x ∈Ω.



Algorithm C: Let the initial datum for the problem (<ref>)–(<ref>) be given as a function u_0∈ L^2(Ω;m) and let h ∈ (0,h_0] be the spatial grid size parameter.



Step 1: Let υ:𝔛→ L^2(Ω;m) be a random field on a probability space (𝔛,F,P), and discretize υ by a finite element approximation of random fields υ^h ∈ V^h_m, so that υ^h(ω;·)_L^2(Ω;m)≤1 for P-a.e. ω∈𝔛. Discretize u_0 by a finite element approximation u_0^h ∈ V^h_m, and then perturb this discretization by defining


    u_0^h,ϵ(ω;x):=u_0^h(x)+ϵυ^h(ω;x),    ω∈𝔛,   x ∈Ω,


where ϵ∈ (0,1], and let σ^h,ϵ be the law of u_0^h,ϵ, meaning that σ^h,ϵ_x(E)=P(u_0^h,ϵ(ω;x)∈ E) for ω∈𝔛 and x ∈Ω.



Step 2: Evolve the initial random field by applying a suitable numerical scheme, with solution map S_t, to the initial datum u_0^h,ϵ(ω) for every ω∈𝔛, obtaining an approximate random field u^h,ϵ(ω,t,·):=S_tu_0^h,ϵ(ω;·), t ∈ (0,T].



Step 3: Define the approximate Young measure solution μ^h,ϵ as the law of u^h,ϵ with respect to P, that is, for all Borel sets E⊂m, 


    μ_t,x^h,ϵ(E)=P(u^h,ϵ(ω;t,x)∈ E),    ω∈𝔛,   (t,x) ∈ Q_T,


and define, analogously, ν^h,ϵ to be the law of Du^h,ϵ with respect to P.



As before we need a method to approximate the random field u_0^ϵ(ω;x), as well as the measures μ^h,ϵ and ν^h,ϵ.



Algorithm D: Let the initial datum for the problem (<ref>)–(<ref>) be given as a function u_0∈ L^2(Ω;m), let h ∈ (0,h_0] be the spatial grid size parameter, and let M∈ℕ.



Step 1: From some probability space (𝔛, F, P) draw M independent and identically distributed random fields υ^h,1,…,υ^h,M such that υ^h,i(ω;·)_L^2(Ω;m)≤1, for P-a.e. ω∈𝔛 and for all i=1,…,M, and such that u_0^h,1,ϵ,…,u_0^h,M,ϵ all have the same law σ^h,ϵ, where


    u_0^h,k,ϵ(ω;x)=u_0^h(x)+ϵυ^h,k(ω;x),    ω∈𝔛,    x ∈Ω,


with ϵ∈ (0,1]. We make explicit here the fact that u_0^h,k,ϵ(ω;·)_L^2(Ω;m)≤u_0^h_L^2(Ω;m)+ϵ for all ω∈𝔛, so that u_0^h,k,ϵ(ω;·)_L^2(Ω;m) can be assumed to be bounded independent of k, ϵ, and h (thanks to the assumed strong convergence of u_0^h to u_0 in L^2(Ω;m) and because ϵ∈ (0,1]).



Step 2: For each k and for a fixed ω∈𝔛 approximate the solution to the initial boundary-value problem (<ref>)–(<ref>) under consideration using the solution operator with initial data u_0^h,k,ϵ(ω); denote u^h,k,ϵ(ω;t,·)=S_tu_0^h,k,ϵ(ω;·,·).



Step 3: Define the approximate Young measure solution (μ_t,x^h,M,ϵ,ν_t,x^h,M,ϵ)
by


    μ_t,x^h,M,ϵ:=1/M∑_k=1^Mδ_u^h,k,ϵ(·;t,x)        ν_t,x^h,M,ϵ:=1/M∑_k=1^Mδ_Du^h,k,ϵ(·;t,x).





 §.§ 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 (<ref>)–(<ref>) under consideration as we pass Δ t → 0_+, M→∞, h,ϵ→0_+ (in that order).



The parameter ω taken from the probability space 𝔛 is intended to represent the “seed" in the numerical algorithm. Different seeds ω∈𝔛 will give rise to different approximate solutions. What will be shown in our analysis of the limit process Δ t → 0_+,
M →∞, h, ϵ→ 0_+ is that one can pass to the limit in M (along a subsequence) to obtain convergence to a quantity which is independent of ω∈𝔛, for P-a.e. ω∈𝔛.





Step 1: Passage to the limits Δ t → 0_+, M→∞. We shall use the semidiscrete numerical scheme discussed in Section <ref>, augmented with the time-stepping procedure discussed in Section <ref>. By Lemma <ref> we have that, for any ω∈𝔛, any ϵ>0, any positive integer k and any h ∈ (0,h_0] and Δ t = T/N, N ∈ℕ, we can construct, via the time-stepping procedure from Section <ref>, a function u^h,k,ϵ_Δ t(ω;·,·), which for any h ∈ (0,h_0] and any Δ t>0 (sufficiently small, for a given fixed h ∈ (0,h_0]) uniquely solves the fully discrete scheme and depends continuously on the discretized initial condition u_0^h,k,ϵ. For the rest of the analysis we shall assume for the sake of brevity that the passage to the limit Δ t→0_+ has already been made for ω∈𝔛 and h ∈ (0,h_0] fixed,
by repeating the analysis performed in Section <ref>, and we shall therefore consider the limiting function u^h,k,ϵ(ω;·,·) resulting from the passage to the limit Δ t → 0_+ instead of u^h,k,ϵ_Δ t(ω;·,·). Passage to the limit Δ t → 0_+, for reasons that will be explained in due course, must happen before the limit passage M→∞. 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 Δ t_j→0_+, which works for each k∈{1,…,M}, is straightforward.

For each initial datum u_0^h,k,ϵ(ω;·) (for fixed ω∈𝔛 and k∈{1,…,M}), we have from
(<ref>) the following energy estimate:


    u^h,k,ϵ(ω)_L^∞(0,T;L^2(Ω;m))^2   +∑_i=1^mDu_i^h,k,ϵ(ω)_L^p_i(0,T;W_0^1,p_i(Ω;m× n))^p_i
       ≤ c(1+ u_0^h,k,ϵ(ω)_L^2(Ω;m)^2+F_L^p'(Q_T;m)^p').


Furthermore, (<ref>) implies, with q̂:=max{q,2}  and q̂'=q̂/q̂-1, that


    a(Du^h,k,ϵ)(ω)_L^q̂'(Q_T;m× n)≤ C


and (<ref>) yields


    ∂_tu^h,k,ϵ(ω)_L^q̂'(0,T;W^-1,q̂'(Ω;m))≤ C,


where C is a positive constant that is independent of h,k and ϵ. 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 ϵ. Summing the equations satisfied by u^h,k,ϵ(ω;·,·) over k=1,…,M and dividing by M we see that the functions u^h,k,ϵ satisfy, for all ψ^h∈ L^q(0,T;V^h_m), the equality:


    ∫_0^T∫_Ω1/M∑_k=1^M ∂_tu^h,k,ϵ(ω;t,x)·ψ^h(t,x)+
    1/M∑_k=1^M a(Du^h,k,ϵ(ω;t,x)): Dψ^h(t,x)xt
            +∫_0^T∫_Ω1/M∑_k=1^MBu^h,k,ϵ(ω;t,x)·ψ^h(t,x)xt =∫_0^T∫_Ω F·ψ^hxt.



Integration by parts in time then yields, for all ψ^h∈ W^1,1_0(0,T;V^h_m),


    ∫_0^T∫_Ω1/M∑_k=1^M u^h,k,ϵ(ω;t,x)· (∂_tψ^h(t,x)-B^ Tψ^h(t,x))xt
            -∫_0^T∫_Ω1/M∑_k=1^M a(Du^h,k,ϵ(ω;t,x)): Dψ^h(t,x)-F(t,x)·ψ^h(t,x)xt=0.



We now wish to let M→∞ in (<ref>). 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 {1/M∑_k=1^M a(Du^h,k,ϵ(ω;·,·))}_M ≥ 1, for a fixed ω∈𝔛, which is,  thanks to (<ref>), bounded in the L^q̂'(Q_T;m × n) norm, uniformly in h ∈ (0,h_0] and ϵ∈ (0,1], and as such there is a subsequence {M_n}_n ≥ 1 and a function χ^h,ϵ(ω;·,·)∈ L^q̂'(Q_T;m × n) such that 1/M_n∑_k=1^M_n a(Du^h,k,ϵ(ω;·,·)) ⇀χ^h,ϵ(ω;·,·) in L^q̂'(Q_T;m × n), as n →∞. We aim to identify this function χ^h,ϵ(ω;·,·) for ω∈𝔛, by identifying the limit in some weaker space and appealing to the uniqueness of weak limits.

Let us consider the initial datum u^h,ϵ_0(ω;x)=u^h_0(x) + ϵυ^h(ω;x) that is independent
of the initial data u_0^h,1,ϵ,…,u_0^h,M,ϵ, and has the same distribution as these initial data (meaning in particular that the law of u_0^h,ϵ is given by σ^h,ϵ), and satisfies υ^h(ω;·)_L^2(Ω;m)≤ 1 for P-almost every ω∈𝔛. We evolve this according to the algorithm for each ω∈𝔛 to obtain a function u^h,ϵ(ω;t,x). We denote the law of u^h,ϵ by μ^h,ϵ, and similarly denote the law of Du^h,ϵ by ν^h,ϵ. From Theorem <ref> we see that the following regularity results hold:


    u^h,ϵ   ∈ L^∞(𝔛;_i=1^mL^p_i(0,T;W_0^1,p_i(Ω;m)))∩ L^∞(𝔛×(0,T);L^2(Ω;m)),
    ∂_t u^h,ϵ   ∈ L^∞(𝔛;L^q̂'(0,T;W^-1,q̂'(Ω;m))),
    a(Du^h,ϵ)   ∈ L^∞(𝔛;L^q̂'(Q_T;m× n)).


Analogously to (<ref>)–(<ref>), the norms of these functions in the respective spaces can be bounded by a positive constant C, independent of h and ϵ. Next, we formulate an analogue of Theorem 11 from <cit.>.



Let α∈(1,p/q-1) and suppose that γ∈ L^α'(Q_T;E_p/α) (understood as an m× n matrix-valued function). Along the subsequence {M_n}_n ≥ 1 we have


    ∫_0^T∫_Ω⟨ν_ω;t,x^h,M_n,ϵ,γ(t,x,·)⟩xt
    →∫_0^T∫_Ω⟨ν_t,x^h,ϵ,γ(t,x,·)⟩ dxt,


in L^2(𝔛) as n→∞. In particular, for each γ∈ L^α'(Q_T;E_p/α) there is a subsequence of M_n such that for P-a.e. ω∈𝔛 the convergence is pointwise.





Given a random field η:𝔛→ L^1(Q_T) we define its expectation with respect to the probability measure P as 


    𝔼(η):=∫_𝔛η(ω)P(ω).


Similarly to <cit.> we define, for
k=1,…,M, i=1,…, m and j=1,…,n the following quantities:


    G_i,j(ω)=∫_0^T∫_Ωγ_i,j(t,x,Du^h,ϵ(ω;t,x))xt,


and


    G^k_i,j(ω)=∫_0^T∫_Ωγ_i,j(t,x,Du^h,k,ϵ(ω;t,x))xt.


We need to show that G^1_i,j,…,G^M_n_i,j are independent and identically distributed.

We begin by remarking that as the initial data u_0^h,1,ϵ,…,u_0^h,M_n,ϵ are independent and identically distributed (in the variable ω∈𝔛), and as the solution operator S_t mapping an initial condition ũ_0 to ũ(t,·) is a continuous map (see Lemma <ref>), the functions u^h,k,ϵ are all independent and identically distributed as well (see Lemma <ref>). Furthermore, by the remarks following the statement of Theorem <ref>, the functions Du^h,k,ϵ are also independent and identically distributed.

Recall that the functions γ_i,j are measurable, they are continuous in their third argument, and satisfy the growth rate coming from the space E_p/α. This implies that the functions γ_i,j(t,x,Du^h,k,ϵ(ω;t,x)) are independent and identically distributed with respect to ω∈𝔛. From Tonelli's Theorem it then follows that the functions G^k_i,j are measurable in ω, and thus by Lemma <ref> are also independent and identically distributed in ω∈𝔛.

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


    (𝔼(G_i,j(ω))-1/M_n∑_k=1^M_n G^k_i,j(ω))^2,


and as in <cit.>, using the fact that the random variables G_i,j^k are independent and identically distributed, we can reduce this to


    𝔼((𝔼(G_i,j(ω))-1/M_n∑_k=1^M_n G^k_i,j(ω))^2)=1/M_n(𝔼((G_i,j)^2)-(𝔼(G_i,j))^2) ≤1/M_n𝔼((G_i,j)^2).


Given α∈(1,p/q-1) we have that


    𝔼((G_i,j)^2)   =∫_𝔛(∫_0^T∫_Ωγ_i,j(t,x,Du^h,ϵ(ω;t,x))xt)^2P(ω)
       ≤∫_𝔛γ_i,j_L^α'(Q_T;E_p/α)^21+|Du^h,ϵ(ω)|^p/α_L^α(Q_T)^2P(ω).



Next, we can bound the norm 1+|Du^h,ϵ(ω)|^p/α_L^α(Q_T) independently of h, ϵ and ω by using the energy estimate (<ref>) (analogously to (<ref>)). Thus, since P is a probability measure, we have that 𝔼((G_i,j)^2) is bounded, and so the right-hand side of (<ref>) converges to zero as M_n→∞.

We wish to apply Lemma <ref> (note also Remark <ref>). In order to do this we must show that


    ∫_𝔛 |γ_i,j(t,x,Du^h,ϵ(ω;t,x))|P(ω)<∞    .


Note that if we show that


    ∫_𝔛∫_0^T∫_Ω |γ_i,j(t,x,Du^h,ϵ(ω;t,x))|xtP(ω)<∞,


then the desired bound is true for almost every (t,x)∈ Q_T, which is sufficient for our purposes. We see, using that γ∈ L^α'(Q_T;E_p/α) with α∈(1,p/q-1), that



    ∫_𝔛∫_0^T∫_Ω    |γ_i,j(t,x,Du^h,ϵ(ω;t,x))|xtP(ω)
       = ∫_𝔛∫_0^T∫_Ω |γ_i,j(t,x,Du^h,ϵ(ω;t,x))|1+|Du^h,ϵ(ω;t,x)|^p/α/1+|Du^h,ϵ(ω;t,x)|^p/αxtP(ω)
       ≤∫_𝔛∫_0^T∫_Ωγ_i,j(t,x,·)_E_p/α(1+|Du^h,ϵ(ω;t,x)|^p/α)xtP(ω)
       ≤∫_𝔛γ_i,j_L^α'(Q_T;E_p/α)1+|Du^h,ϵ(ω;·,·)|^p/α_L^α(Q_T)P(ω)
       ≤γ_i,j_L^α'(Q_T;E_p/α)∫_𝔛 |Q_T|^1/α+Du^h,ϵ(ω;·,·)_L^p(Q_T;m × n)^p/αP(ω)
       ≤γ_i,j_L^α'(Q_T;E_p/α)(|Q_T|^1/α+Du^h,ϵ_L^∞(𝔛;L^p(Q_T;m × n))^p/α),


which is finite by the assumptions on γ and by (<ref>). Then, by applying Lemma <ref>, we see that: 

    𝔼(G_i,j(ω))   =∫_𝔛∫_0^T∫_Ωγ_i,j(t,x,Du^h,ϵ(ω;t,x))xtP(ω)
       =∫_0^T∫_Ω(∫_𝔛γ_i,j(t,x,Du^h,ϵ(ω;t,x))P(ω))xt
       =∫_0^T∫_Ω⟨ν^h,ϵ_i,j;t,x,γ_i,j(t,x,·)⟩xt.



Repeating these calculations for all components γ_i,j gives us the convergence of 


    ∫_0^T∫_Ω⟨ν_ω;t,x^h,M_n,ϵ,γ(t,x,·)⟩ dxt→∫_0^T∫_Ω⟨ν_t,x^h,ϵ,γ(t,x,·)⟩xt,


in L^2(𝔛) as M_n→∞, from which we deduce the existence of a subsequence (not relabelled) converging for P- almost every ω∈𝔛 using standard results in measure theory.



What we would like to conclude is that by the above theorem we have, for P-almost every ω∈𝔛, the weak-star convergence of ν^h,M_n,ϵ_ω=ν^h,M_n,ϵ_ω,·,· to ν^h,ϵ=ν^h,ϵ_·,· as n →∞. 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 γ (a point not clarified in <cit.> or <cit.>). To get around this we note that as E_p/α is separable (see <cit.>), the space L^α'(Q_T;E_p/α) is also separable for our choice of α. Therefore there is a countable dense subset {γ_i}_i=1^∞. The desired weak-star convergence then follows from a diagonal argument.

Recall from the calculations performed prior to Lemma <ref> that, as M_n →∞ (n →∞),


    ∫_0^T∫_Ω1/M_n∑_k=1^M_na(Du^h,k,ϵ(ω;t,x)) : Dψ(t,x)xt→∫_0^T∫_Ωχ^h,ϵ(ω;t,x) : Dψ(t,x)xt


for every ψ∈ L^q̂([0,T];W^1,q̂_0(Ω;m)), and for all ω∈𝔛. Furthermore, by Lemma <ref>, there is a subsequence of M_n such that, as M_n →∞ (n →∞),


    ∫_0^T∫_Ω1/M_n∑_k=1^M_na(Du^h,k,ϵ(ω;t,x)) : Dψ(t,x)xt→∫_0^T∫_Ω⟨ν^h,ϵ_t,x,a⟩ : Dψ(t,x)xt


for every ψ∈ L^α'(0,T;W^1,α'_0(Ω;m)), and for P-a.e. ω∈𝔛.

As both (<ref>) and (<ref>) hold for ψ∈ L^max{q̂,α'}(0,T;W^1,max{q̂,α'}_0(Ω;m)), by the uniqueness of the weak limit we can identify χ^h,ϵ(ω;·,·) =⟨ν^h,ϵ,a⟩ for P-a.e. ω∈𝔛. As the right-hand side of this equality is independent of ω∈𝔛, it follows that χ^h,ϵ is also independent of ω, and we shall therefore suppress the dependence on ω in our notation and write χ^h,ϵ(t,x) instead of χ^h,ϵ(ω;t,x), noting that χ^h,ϵ = ⟨ν^h,ϵ,a⟩∈ L^q̂'(Q_T;m × n). Thus, by combining (<ref>) and (<ref>), which guarantee weak convergence in L^q̂'(0,T;W^1,q̂'_0(Ω;m)) and L^α(0,T;W^1,α_0(Ω;m)), respectively, for P-a.e. ω∈𝔛, and therefore, by uniqueness of the weak limit, also in the function space L^max{q̂',α}(0,T;W^1,max{q̂',α}_0(Ω;m)), for P-a.e. ω∈𝔛, we have that


    ∫_0^T∫_Ω1/M_n∑_k=1^M_na(Du^h,k,ϵ(ω;t,x)) : Dψ(t,x)xt→∫_0^T∫_Ω⟨ν^h,ϵ_t,x,a⟩ : Dψ(t,x)xt


for all ψ∈ L^min{q̂,α'}(0,T;W^1,min{q̂,α'}_0(Ω;m)), for P-a.e. ω∈𝔛.

The terms involving the measures μ^h,ϵ are treated completely analogously. Thus, by passing to the limit M_n →∞ (n →∞) in (<ref>), we have that for all ψ^h∈ W^1,1_0(0,T;V^h_m) the following identity is satisfied:


    ∫_0^T∫_Ω⟨μ^h,ϵ_t,x,ξ⟩·(∂_tψ^h(t,x)-B^ Tψ^h(t,x))-⟨ν^h,ϵ_t,x,a⟩ : Dψ^h(t,x)+F(t,x)·ψ^h(t,x)xt=0.





Step 2: Passing h,ϵ→0_+.
By Lemma <ref> and noting that a ∈ E_p/α, it follows that ⟨ν^h,M_n,ϵ_ω,·,·,a ⟩⇀⟨ν^h,ϵ_·,·,a ⟩
weakly in L^α(Q_T;m× n), for α∈ (1,p/q-1), and for P-a.e. ω∈𝔛. Thus, by weak lower-semicontinuity of the norm function, we have that ⟨ν^h,ϵ_·,·,a⟩_L^α(Q_T;m× n)≤lim inf_n→∞⟨ν^h,M_n,ϵ_ω,·,·,a⟩_L^α(Q_T;m× n) for P-a.e. ω∈𝔛.

Let 1<α<p/q-1 and consider the function ã(ξ)=1+|ξ|^p/α; then, ã∈ E_p/α. It then follows from Lemma <ref> that ⟨ν^h,M_l,ϵ_ω,·,·,ã⟩⇀⟨ν^h,ϵ,ã⟩ in L^α(Q_T;m× n) as l→∞ for P-a.e. ω∈𝔛. Then, by weak lower-semicontinuity of the norm function ·_L^α(Q_T;m× n), we have that


    ⟨ν^h,ϵ,ã⟩_L^α(Q_T;m × n)   ≤lim inf_l→∞⟨ν^h,M_l,ϵ_ω,·,·,ã⟩_L^α(Q_T;m × n)
        = lim inf_l→∞1/M_l∑_k=1^M_l (1 + |Du^h,k,ϵ(ω;·,·)|^p/α) _L^α(Q_T)
       ≤lim inf_l→∞1/M_l∑_k=1^M_l1 + |Du^h,k,ϵ(ω;·,·)|^p/α_L^α(Q_T)
       ≤lim inf_l→∞ |Q_T|^1/α + 1/M_l∑_k=1^M_lDu^h,k,ϵ(ω;·,·)_L^p(Q_T;m× n)^p/α
       ≤lim inf_l→∞ |Q_T|^1/α + 1/M_l∑_k=1^M_l[c(1+ u_0^h,k,ϵ(ω)_L^2(Ω;m)^2+F_L^p'(Q_T;m)^p')]^1/α
       ≤lim inf_l→∞ |Q_T|^1/α + 1/M_l∑_k=1^M_l[c(1+ 2u_0^h_L^2(Ω;m)^2+ 2ϵ^2 + F_L^p'(Q_T;m)^p')]^1/α
        =  |Q_T|^1/α + c (3+ 2u_0^h_L^2(Ω;m)^2 + F_L^p'(Q_T;m)^p')^1/α≤ C,


where C is a positive constant, independent of h ∈ (0,h_0] and ϵ∈ (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 (<ref>), followed by recalling that u_0^h,k,ϵ(ω)_L^2(Ω;m)≤u_0^h_L^2(Ω;m) + ϵ, with ϵ∈ (0,1], and noting that by the assumed strong convergence of u_0^h to u_0 in L^2(Ω;m) one can bound u_0^h_L^2(Ω;m) by a constant, independent of h ∈ (0,h_0].

Having shown that ⟨ν^h,ϵ,ã⟩_L^α(Q_T) is bounded by a constant,
independent of h and ϵ, we now compute


    ν^h,ϵ_L^α(Q_T;E_p/α')   =sup_γ_L^α'(Q_T;E_p/α)=1 ∫_Q_T⟨ν^h,ϵ_t,x,γ(t,x,·)⟩xt
       =sup_γ_L^α'(Q_T;E_p/α)=1 ∫_Q_T⟨ν^h,ϵ_t,x,γ(t,x,·)(1+|·|^p/α)/1+|·|^p/α⟩xt
       =sup_γ_L^α'(Q_T;E_p/α)=1 ∫_Q_T∫_A ∈ℝ^m × nγ(t,x,A)(1+|A|^p/α)/1+|A|^p/αν^h,ϵ_t,x(A) xt
       ≤sup_γ_L^α'(Q_T;E_p/α)=1∫_Q_Tγ(t,x,·)_E_p/α
    |⟨ν^h,ϵ_t,x,ã⟩| xt
       ≤⟨ν^h,ϵ,ã⟩_L^α(Q_T),


which we can further bound from above by a positive constant C, independent of h ∈ (0,h_0] and ϵ∈ (0,1], using (<ref>). Hence, for each ϵ∈ (0,1] there is a subsequence of the sequence of measures {ν^h,ϵ} indexed by h_j and a measure ν^ϵ∈ L^α(Q_T;E_p/α') such that for every γ∈ E_p/α we have that


    ⟨ν^h_j,ϵ,γ⟩⇀⟨ν^ϵ,γ⟩


in L^α(Q_T;m× n) as j→∞. In particular this holds for γ=a ∈ E_p/α, 1<α<p/q-1, that is:


    ⟨ν^h_j,ϵ,a⟩⇀⟨ν^ϵ,a⟩


in L^α(Q_T;m× n) for any such α>1. Similar calculations hold for the family of measures μ^h,ϵ.

Taking the final limit in ϵ 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,ϵ) and pass to the limit along this subsequence, obtaining the existence of a two measures μ and ν satisfying, for all ψ∈ W^1,1_0(0,T;W^1,α'_0(Ω;m)) (contained in W^1,1_0(0,T;W^1,min{q̂,α'}_0(Ω;m))⊂ L^min{q̂,α'}(0,T;W^1,min{q̂,α'}_0(Ω;m))):


    ∫_0^T∫_Ω⟨μ_t,x,ξ⟩·(∂_tψ(t,x)-B^ Tψ(t,x))-⟨ν_t,x,a⟩ : Dψ(t,x)+F(t,x)·ψ(t,x)xt=0.





Step 3: Interpretation as a Young measure solution. Next, we discuss how the pair of measures (μ,ν) satisfying (<ref>) can be can be interpreted as a Young measure solution. Recall that μ^h,ϵ is the law of the random variable u^h=u^h,ϵ(ω;t,x), and that, for each ω∈𝔛 and each t ∈ [0,T], the function u^h,ϵ(ω;t,·) ∈ V^h_m ⊂ W^1,∞_0(Ω;m) is a linear combination of finite element basis functions over the domain Ω. Therefore,


    x ∈Ω↦⟨μ^h,ϵ_t,x,ξ⟩=∫_𝔛 u^h,ϵ(ω;t,x)P(ω),      t ∈ [0,T],


is weakly differentiable with respect to x ∈Ω, and thanks to the definition of weak derivative and Fubini's theorem,


    D⟨μ^h,ϵ_t,x,ξ⟩=∫_𝔛Du^h,ϵ(ω;t,x)P(ω)=⟨ν^h,ϵ_t,x,ξ⟩,     (t,x) ∈ (0,T] ×Ω.


Using the weak-star convergence of μ^h,ϵ and ν^h,ϵ to, respectively, μ and ν established in Step 2, in conjunction with the regularity of u^h,ϵ coming from (<ref>)–(<ref>), and the fact that 𝔛 is a space of finite P-measure, we see that by passing h,ϵ→0_+ along subsequences, we get the existence of a function


    u∈ L^∞(𝔛;_i=1^mL^p_i(0,T;W_0^1,p_i(Ω;m)))


such that, for (t,x) ∈ (0,T] ×Ω,


    ⟨μ_t,x,ξ⟩=∫_𝔛 u(ω;t,x)P(ω)


and


    ⟨ν_t,x,ξ⟩=∫_𝔛 Du(ω;t,x)P(ω).


This comes from the fact that, for example,


    ⟨μ^h,ϵ,ξ⟩⇀⟨μ,ξ⟩,


and we have that


    ∫_0^T∫_Ω∫_𝔛u^h,ϵ(ω;t,x)·ϕ(ω;t,x)P(ω)xt
    →∫_0^T∫_Ω∫_𝔛u(ω;t,x)·ϕ(ω;t,x)P(ω)xt


for all ϕ∈ L^1(𝔛; L^α(0,T; W_0^1,α(Ω;m))). However, as 𝔛 has finite P-measure, we can choose test functions which are constant in the parameter ω to see that 


    ∫_0^T∫_Ω(∫_𝔛u^h,ϵ(ω;t,x)P(ω))·ϕ(t,x)xt
    →∫_0^T∫_Ω(∫_𝔛u(ω;t,x)P(ω))·ϕ(t,x)xt


for all ϕ∈ L^α(0,T;W_0^1,α(Ω;m)), which gives the desired representation for ⟨μ_t,x,ξ⟩ (the corresponding term for ⟨ν_t,x,ξ⟩ follows similarly).

Now we define the function


    U(t,x):=∫_𝔛 u(ω;t,x)P(ω)=⟨μ_t,x,ξ⟩.


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)=⟨ν_t,x,ξ⟩ in a suitable sense, then we will have shown that the pairing (U,ν) is a Young measure solution of the system (<ref>)–(<ref>) 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 ∈ L^p(0,T;W^1,p_0(Ω;m)) satisfies u|_∂Ω=0, and so we can extend the function u (and hence the function U) by zero from [0,T]×Ω to the whole of [0,T]×n so that the extended
function (still denoted by u) belongs to L^p(0,T;W^1,p(n;m)). We consider the difference quotient of U, defined by


    D^δ_i U=U(t,x+δ e_i)-U(t,x)/δ,


where e_i is the unit vector in the i-th co-ordinate direction. We then have that


    D^δ_i U(t,·)_L^p_i(Ω;m)   =U(t,·+δ e_i)-U(t,·)/δ_L^p_i(Ω;m)
       =∫_𝔛u(ω;t,·+δ e_i)-u(ω;t,·)/δP(ω)_L^p_i(Ω;m)
       ≤∫_𝔛u(ω;t,·+δ e_i)-u(ω;t,·)/δ_L^p_i(Ω;m)P(ω)
       ≤∫_𝔛D_iu(ω;t,·)_L^p_i(Ω;m)P(ω),     ∀ δ∈ (0,1),     i=1,…,n.


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


    u(ω;t,·+δ e_i)-u(ω;t,·)/δ_L^p_i(Ω;m)^p_i= ∫_0^1 D_i u(ω;t,·+sδ)  s^p_i_L^p_i(Ω;m)≤∫_0^1  D_i u(ω;t,·+sδ)^p_i_L^p_i(Ω;m) s
           = ∫_0^1 ∫_Ω |D_i u(ω;t,x+sδ)|^p_i x  s ≤∫_0^1 ∫_Ω |D_i u(ω;t,x)|^p_i x  s = D_i u(ω;t,·)^p_i_L^p_i(Ω;m).


Therefore, by a standard characterization of Sobolev functions in terms of difference quotients, D_iU(t,·)∈ L^p_i(Ω;m × n), i=1,…,n, for almost every t∈[0,T]. By integrating over t ∈ [0,T] in the above inequality it follows that D_iU∈ L^p_i(Q_T;m × n), i=1,…,n, and in particular U ∈ L^p(0,T;W^1,p_0(Ω;m)), as required. By an analogous argument,


    U ∈ L^∞(0,T;L^2(Ω;m)),    ∂_t U ∈ L^q̂'(0,T;W^-1,q̂'(Ω;m)),


thanks to the regularity results


    u(ω;·,·) ∈ L^∞(0,T;L^2(Ω;m)),    ∂_t u(ω;·,·) ∈ L^q̂'(0,T;W^-1,q̂'(Ω;m)),


which follow from the regularity properties of u^h,ϵ stated in  (<ref>)–(<ref>) using the weak lower-semi-continuity of the norm function.

Next, we write


    V(t,x):=⟨ν_t,x,ξ⟩=∫_𝔛 Du(ω;t,x)P(ω),


and let φ∈ C^∞_0(Q_T;m × n) be a test function. We then compute, using Fubini's Theorem, with D = D_x and = _x,


    ∫_0^T∫_Ω DU(t,x) : φ(t,x)xt   =-∫_0^T∫_Ω U(t,x) · φ(t,x)xt
       =-∫_0^T∫_Ω∫_𝔛 u(ω;t,x) · φ(t,x)P(ω)xt
       =-∫_𝔛∫_0^T∫_Ω u(ω;t,x) · φ(t,x)xtP(ω)
       =∫_𝔛∫_0^T∫_Ω Du(ω;t,x) : φ(t,x)xtP(ω)
       =∫_0^T∫_Ω∫_𝔛 Du(ω;t,x) : φ(t,x)P(ω)xt
       =∫_0^T∫_Ω V(t,x) : φ(t,x)xt     ∀ φ∈ C^∞_0(Q_T;m × n).


Hence, DU(t,x) = V, and because V:= ⟨ν, ξ⟩, it follows that DU = ⟨ν, ξ⟩, as desired. Thus we have shown that the pair (U,ν) is a Young measure solution of the system (<ref>)–(<ref>) under consideration.



Step 4: Attainment of the initial condition. As


    ∂_tU^h,ϵ∈ L^q̂'(0,T;W^-1,q̂'(Ω;m)),


and


    U^h,ϵ∈_i=1^m L^p_i(0,T;W_0^1,p_i(Ω;m))∩ L^∞(0,T;L^2(Ω;m)),


it follows from Lemma <ref> that there is a subsequence (h_j,ϵ_j) of (h,ϵ) along which


    ∫_Ω U^h_j,ϵ_j(t,x)·φ(x)x→∫_Ω U(t,x)·φ(x)x


uniformly in C([0,T]) for all φ∈ L^2(Ω;m). Therefore to show that U(0,·) = u_0 we let φ∈ L^2(Ω;m) and compute:


    ∫_Ω (U(0,x)-u_0(x))·φ(x)x   =∫_Ω (U(0,x)-U^h_j,ϵ_j(0,x))·φ(x)x
            + ∫_Ω (U^h_j,ϵ_j(0,x)-U^h_j,ϵ_j(t,x))·φ(x)x
            + ∫_Ω(U^h_j,ϵ_j(t,x)-∫_𝔛u_0^h_j,ϵ_j(ω,x)P(ω))·φ(x)x
            + ∫_Ω(∫_𝔛u_0^h_j,ϵ_j(ω,x)P(ω)-u_0(x))·φ(x)x
       =:I_j+II_j(t)+III_j(t)+IV_j.


The term I_j converges to zero as j→∞ by the uniform convergence in C([0,T]) stated in (<ref>). The term II_j(t) converges to zero as t→0_+ as it holds that u^h_j,ϵ_j(ω;t,x)→ u^h_j,ϵ_j(ω,0,x) by continuity of t ∈ [0,T] ↦ u^h_j,ϵ_j(ω;t,x), ω∈𝔛, x ∈Ω, and we can apply the Dominated Convergence Theorem to interchange this limit with the integral over 𝔛 that appears in the definition of U^h_j,ϵ_j. The term III_j(t) also converges to zero as t→0_+ as we have that the numerical functions u^h_j,ϵ_j satisfy the initial condition, and we may once again apply the Dominated Convergence Theorem to interchange the limit with the integral over 𝔛. Finally for the term IV_j we compute:


    |∫_𝔛 u_0^h_j,ϵ_j(ω,x)P(ω)-u_0(x)|
       =|∫_𝔛u_0^h_j+ϵ_jυ^h_j(ω;x)-u_0(x)P(ω)|
       ≤ |u_0^h(x)-u_0(x)|+ϵ_j∫_𝔛|υ^h_j(ω,x)|P(ω).


Integrating this over Ω we get


    ∫_Ω|∫_𝔛 u_0^h_j,ϵ_j(ω,x)P(ω)-u_0(x)|x   ≤u_0^h_j-u_0_L^1(Ω;m)+ϵ_j ∫_𝔛υ^h_j(ω,·)_L^1(Ω;m)P(ω)
       ≤u_0^h_j-u_0_L^2(Ω;m)+ϵ_j |Ω|^1/2∫_𝔛υ^h_j(ω,·)_L^2(Ω;m)P(ω)
       ≤u_0^h_j-u_0_L^2(Ω;m)+ϵ_j |Ω|^1/2,


which converges to zero as j→∞ by the assumed strong convergence of our discretized initial condition u_0^h to u_0 in L^2(Ω;m). In summary then,  lim_j→∞ I_j =0, lim_t→ 0_+ II_j(t)=0 for all j ≥ 1,  lim_t→ 0_+ III_j(t)=0 for all j ≥ 1, and  lim_j→∞ IV_j =0. Passing to the limit t → 0_+ in (<ref>) with j ≥ 1 kept fixed and then with j →∞ we deduce that


    ∫_Ω (U(0,x)-u_0(x))·φ(x)x = 0     ∀ φ∈ L^2(Ω;m).


Thus we have that U(0,x)=u_0(x) for almost every x∈Ω, showing that the initial condition is satisfied.

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



Algorithms C and D converge to a Young measure solution of the problem (<ref>)–(<ref>) in the event that the initial Young measure is a Dirac mass, δ_u_0(x), concentrated on an initial datum u_0∈ L^2(Ω;m).





If we have a (nonatomic) measure-valued initial datum σ, 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 σ 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 σ. The existence of such a random variable is guaranteed by Proposition <ref>.




 In <cit.> and <cit.>, 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 ϵ tending to zero and M→∞ are taken. For us here the passages to the limits in h→ 0_+, ϵ→ 0_+ and Δ t → 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 <ref> 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,ϵ. This required that the solution operator 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,ϵ for each fixed h, and were unable to transfer this property to the limiting function as we pass h→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 S_t was measurable when considered as a mapping from the nondiscretized initial condition u_0∈ L^2(Ω;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→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].


plain


