Homogenization of time-fractional diffusion equations with periodic coefficients

TL;DR

This paper studies the homogenization of time-fractional diffusion equations with periodic coefficients, deriving convergence rates for the approximation and validating results through numerical tests.

Contribution

It provides the first order approximation and convergence rates for homogenized solutions of time-fractional diffusion equations with periodic coefficients.

Findings

Convergence rate of O(ε^{1/2}) for dimensions d≤2

Convergence rate of O(ε^{1/6}) for dimension d=3

Numerical tests confirm the effectiveness of the approximation

Abstract

We consider the initial boundary value problem for the time-fractional diffusion equation with a homogeneous Dirichlet boundary condition and an inhomogeneous initial data $a(x)\in L^{2}(D)$ in a bounded domain $D\subset \mathbb{R}^d$ with a sufficiently smooth boundary. We analyze the homogenized solution under the assumption that the diffusion coefficient $\kappa^{\epsilon}(x)$ is smooth and periodic with the period $\epsilon>0$ being sufficiently small. We derive that its first order approximation has a convergence rate of $\mathcal{O}(\epsilon^{1/2})$ when the dimension $d\leq 2$ and $\mathcal{O}(\epsilon^{1/6})$ when $d=3$. Several numerical tests are presented to show the performance of the first order approximation.

Figures

Tables

Table 1. Table 1: The convergence history of the first order approximation to Problem ( 1 ) with κ 𝜅 \kappa in ( 25 ) and ϵ := 1 8 assign italic-ϵ 1 8 \epsilon:=\frac{1}{8} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	2.1235e-9	1.3488e-5	1.4463e-7	2.0673e-4
0.5	4.5471e-10	1.3637e-5	3.0813e-8	2.0794e-4
1	2.411e-10	1.3656e-5	1.6325e-8	2.0809e-4

Table 2. Table 2: The convergence history of the first order approximation to Problem ( 1 ) with κ 𝜅 \kappa in ( 26 ) and ϵ := 1 8 assign italic-ϵ 1 8 \epsilon:=\frac{1}{8} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	1.6061e-8	1.16596e-4	1.0832e-6	1.7659e-3
0.5	3.4581e-9	1.1753e-4	2.3239e-7	1.7737e-3
1	1.8342e-9	1.1765e-4	1.2321e-7	1.7747 e-3

Table 3. Table 3: The convergence history of the first order approximation to Problem ( 1 ) with κ 𝜅 \kappa in ( 25 ) and ϵ := 1 16 assign italic-ϵ 1 16 \epsilon:=\frac{1}{16} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	5.2208e-10	3.3162e-6	7.1850e-8	1.0270e-4
0.5	1.1117e-10	3.3341e-6	1.5262e-8	1.0300e-4
1	5.8895e-11	3.3363e-6	8.0828e-9	1.0303e-4

Table 4. Table 4: The convergence history of the first order approximation to Problem ( 1 ) with κ 𝜅 \kappa in ( 25 ) and ϵ := 1 32 assign italic-ϵ 1 32 \epsilon:=\frac{1}{32} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	1.4377e-10	9.1322e-7	3.8067e-8	5.4412e-5
0.5	3.0547e-11	9.1614e-7	8.0778e-9	5.4515e-5
1	1.6179e-11	9.1649e-7	4.2776e-9	5.4527e-5

Table 5. Table 5: The convergence history of the first order approximation to Problem ( 1 ) with κ := κ 1 assign 𝜅 subscript 𝜅 1 \kappa:=\kappa_{1} in ( 27 ) and ϵ := 1 8 assign italic-ϵ 1 8 \epsilon:=\frac{1}{8} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	8.3877e-10	2.5129e-2	5.6728e-9	3.8002e-2
0.5	1.8973e-10	2.5130e-2	1.2832e-9	3.8003e-2
1	1.0121e-10	2.5130e-2	6.8450e-10	3.8003e-2

Table 6. Table 6: The convergence history of the first order approximation to Problem ( 1 ) with κ := κ 2 assign 𝜅 subscript 𝜅 2 \kappa:=\kappa_{2} in ( 28 ) and ϵ := 1 8 assign italic-ϵ 1 8 \epsilon:=\frac{1}{8} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	5.5466e-9	1.6618e-1	3.6977e-8	2.4771e-1
0.5	1.2546e-9	1.6618e-1	8.3640e-9	2.4772e-1
1	6.6929e-10	1.6618e-1	4.4618e-9	2.4772e-1

Table 7. Table 7: The convergence history of the first order approximation to Problem ( 1 ) with a ​ ( x ) 𝑎 𝑥 a(x) defined in ( 29 ), κ 𝜅 \kappa in ( 25 ) and ϵ := 1 8 assign italic-ϵ 1 8 \epsilon:=\frac{1}{8} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	2.0111e-7	1.8092e-4	5.1832e-6	8.3829e-4
0.5	4.4356e-8	1.8618e-4	1.1545e-6	8.6128e-4
1	2.3592e-8	1.8678e-4	6.1475e-7	8.6391e-4

Table 8. Table 8: The convergence history of the first order approximation to Problem ( 1 ) with a ​ ( x ) 𝑎 𝑥 a(x) defined in ( 29 ), κ 𝜅 \kappa in ( 25 ) and ϵ := 1 16 assign italic-ϵ 1 16 \epsilon:=\frac{1}{16} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	5.3061e-8	4.7739e-5	2.7039e-6	4.3744e-4
0.5	1.1742e-8	4.9292e-5	6.0373e-7	4.5053e-4
1	6.24783e-9	4.9471e-5	3.2156e-7	4.5204e-4

Table 9. Table 9: The convergence history of the first order approximation to Problem ( 1 ) with a ​ ( x ) 𝑎 𝑥 a(x) defined in ( 29 ), κ 𝜅 \kappa in ( 25 ) and ϵ := 1 32 assign italic-ϵ 1 32 \epsilon:=\frac{1}{32} .

$t$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{L^2\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{L^2\left(D\right)}}$	${\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}$	$\frac{{\Vert u_h^{ϵ}-U_{1,h}^{ϵ}\Vert }_{H^1\left(D\right)}}{{\Vert u_h^{ϵ}\Vert }_{H^1\left(D\right)}}$
0.1	1.3154e-8	1.1835e-5	1.3839e-6	2.2390e-4
0.5	2.9144e-9	1.2235e-5	3.0918e-7	2.3074e-4
1	1.5510e-9	1.2281e-5	1.6469e-7	2.3153e-4

Equations

$$\{\begin{array}{rlrl}{\partial }_t^{\alpha }{u}^{ϵ}(x,t)& =\nabla \cdot \left({\kappa }^{\epsilon }(x)\nabla u^{ϵ}(x,t)\right)& & \text{in }D,t\in \left(0,T\right]\\ u^{ϵ}& =0& & \text{on }\partial D,t\in \left(0,T\right]\\ u^{ϵ}(0)& =a(x)& & \text{in }D.\end{array}$$

$${\kappa }_{ij}(x){\xi }_i{\xi }_j\ge \mu |\xi {|}^2\text{ for all }\xi \in {\mathbb{R}}^d\text{ and }x\in D.$$

$${\partial }_t^{\alpha }w(t)=\frac{1}{\Gamma (1-\alpha )}{\int }_0^t\frac{1}{(t-s{)}^{\alpha }}{w}^{\prime }(s)\mathrm{ds}.$$

$$u^{ϵ}(x,t)=u_0(x,y,t)+ϵu_1(x,y,t)+{ϵ}^2{u}_2(x,y,t)+\cdots$$

$$A^{ϵ}=-\frac{\partial }{\partial x_i}\left({\kappa }_{ij}\left(x/ϵ\right)\frac{\partial }{\partial x_j}\right).$$

$$A^{ϵ}={ϵ}^{-2}{A}_1+{ϵ}^{-1}{A}_2+{ϵ}^0{A}_3,$$

$$A_1$$

$$A_2$$

$$A_3$$

$$A_1{u}_0$$

$$A_1{u}_1$$

$$A_1{u}_2$$

$$-\frac{\partial }{\partial y_i}\left({\kappa }_{ij}(y)\frac{\partial }{\partial y_j}\right)u_0(x,y,t)=0.$$

$$u_0(x,y,t)=u_0(x,t).$$

$$u_0(0)$$

$$u_0$$

$$-\frac{\partial }{\partial y_i}\left({\kappa }_{ij}(y)\frac{\partial }{\partial y_j}\right)u_1=\left(\frac{\partial }{\partial y_i}{\kappa }_{ij}(y)\right)\frac{\partial u_0}{\partial x_j}(x,t).$$

$$-\frac{\partial }{\partial y_i}\left({\kappa }_{ij}(y)\frac{\partial }{\partial y_j}\right){\chi }_j=\frac{\partial }{\partial y_i}{\kappa }_{ij}(y)\text{ in }Y.$$

$$u_1(x,y,t)={\chi }_j(y)\frac{\partial u_0}{\partial x_j}(x,t)+{\tilde{u}}_1(x,t).$$

$$\frac{\partial }{\partial y_i}\left({\kappa }_{ij}(y)\frac{\partial }{\partial y_j}\right)u_2=A_2{u}_1+A_3{u}_0+{\partial }_t^{\alpha }{u}_0.$$

$${\int }_Y(A_2{u}_1+A_3{u}_0+{\partial }_t^{\alpha }{u}_0)\mathrm{d}y=0.$$

$${\int }_Y\frac{\partial }{\partial y_i}F(x,y,t)\mathrm{d}y=0$$

$${\partial }_t^{\alpha }{u}_0-\frac{\partial }{\partial x_i}\left(⟨{\kappa }_{ij}(y)⟩\frac{\partial }{\partial x_j}{u}_0\right)-\frac{\partial }{\partial x_i}\left(⟨{\kappa }_{ij}(y)\frac{\partial }{\partial y_j}{u}_1⟩\right)=0.$$

$${\partial }_t^{\alpha }{u}_0-\frac{\partial }{\partial x_i}\left({\kappa }_{ij}^{\ast }\frac{\partial }{\partial x_j}\right)u_0=0,$$

$${\kappa }_{ij}^{\ast }=\frac{1}{|Y|}\left({\int }_Y{\kappa }_{ij}+{\kappa }_{ik}\frac{\partial {\chi }_j}{\partial y_k}\mathrm{d}y\right).$$

$$\{\begin{array}{rlrl}{\partial }_t^{\alpha }{u}_0& =\frac{\partial }{\partial x_i}\left({\kappa }_{ij}^{\ast }\frac{\partial }{\partial x_j}\right)u_0& & \text{ in }D,t\in \left(0,T\right]\\ u_0& =0& & \text{ on }\partial D,t\in \left(0,T\right]\\ u_0(0)& =a(x)& & \text{ in }D.\end{array}$$

$$t^{\alpha }{\Vert u_0(\cdot ,t)\Vert }_{H^2(D)}+{\Vert u_0(\cdot ,t)\Vert }_{L^2(D)}+t^{\alpha }{\Vert {\partial }_t^{\alpha }{u}_0(\cdot ,t)\Vert }_{L^2\left(D\right)}$$

$${\Vert \nabla {\chi }_j\Vert }_{L^2(Y)}$$

$$\Vert {\chi }_j{\Vert }_{H^2(Y)}+\Vert {\chi }_j{\Vert }_{L^{\infty }(Y)}$$

$${\int }_Y\kappa (y)\nabla {\chi }_j\cdot \nabla v\mathrm{d}y=-{\int }_Y\kappa (y)e_j\cdot \nabla v\mathrm{d}y,\text{ for all }v\in V_{\#}(Y).$$

$${\int }_Y\kappa (y)\nabla {\chi }_j\cdot \nabla {\chi }_j\mathrm{d}y=-{\int }_Y\kappa (y)e_j\cdot \nabla {\chi }_j\mathrm{d}y.$$

$$\Vert {\chi }_j{\Vert }_{H^2(Y)}\lesssim \Vert \nabla \cdot (\kappa \nabla {\chi }_j){\Vert }_{L^2(Y)}+\Vert \Delta {\chi }_j{\Vert }_{L^2(Y)}\lesssim 1.$$

$$||D^{\ell }{\zeta }^{ϵ}|{|}_{L^{\infty }(D)}\lesssim {\epsilon }^{-\ell }.$$

Full text

Keywords: time-fractional diffusion; homogenization; 2-scale asymptotic expansion; first order approximation; error estimate

1 Introduction

Let $D$ be a bounded domain in $\mathbb{R}^{d}\,(d=1,2,3)$ with a sufficiently smooth boundary $\partial D$. We consider a partial differential equation with the fractional derivative in time $t$, satisfying:

[TABLE]

Here, $0<\alpha<1$ is a given fixed parameter and the matrix $\kappa^{\varepsilon}({x}):=\kappa(x/\epsilon)$ is periodic with period $\varepsilon$. Let $(\kappa_{ij}(x))_{i,j=1}^{d}$ be symmetric with $\kappa_{ij}(x)\in C^{\infty}(D)$, $i,j=1,\ldots,d$. We assume for some constant $\mu>0$, there holds

[TABLE]

The initial data $a(x)\in L^{2}(D)$ is a given macro-scale function and $T>0$ is a fixed value.

In the model problem (1), ${\partial^{\alpha}_{t}}w$ refers to the left-sided Caputo fractional derivative of order $\alpha$ of the function $w(t)$, defined by (see, e.g. [6, p. 91, (2.4.1)] or [13, p. 78])

[TABLE]

Fractional diffusion equations were introduced in physics with the aim of describing diffusions in media with fractal geometry [11]. They have been applied to many fields, e.g., in engineering, physics, biology and finance. Their practical applications include electron transport in Xerox photocopier, visco-elastic materials, and protein transport in cell membranes [15, 1, 7]. In this paper, we are concerned with the fractional diffusion problem (1) in a heterogeneous periodic media $\kappa^{\epsilon}(x)$, which is utilized in many important applications, for instance, porous media and composite material modeling. Most recently, periodic structures are utilized in metamaterials [16] to design novel materials. The main challenge in the classical numerical treatment of these applications is that it becomes prohibitively expensive and even intractable at the microscale as $\epsilon\to 0$.

The goal of this paper is to construct an efficient numerical solver for (1) based on homogenization theory [3]. The main idea of homogenization is to obtain the effective or homogenized problem by solving $d$ cell problems, and then the corresponding homogenized solution $u_{0}$ serves as a good approximation to original unknown $u^{\epsilon}$ as the period $\epsilon\to 0$. To the best of our knowledge, there has been no such result for time-fractional diffusion problems so far. However, there are quite a few results on the parabolic equations, i.e., $\alpha=1$ in (1), for the same setting, cf. [12, 17]. Nevertheless, it is nontrivial to generalize the results for the parabolic equations to the time-fractional equations with the same technique, mainly due to the lack of the product rule for the fractional derivative.

We prove in Theorem 3.4 the error between the exact solution $u^{\epsilon}$ and its first order approximation $u^{\epsilon}_{1}$ is of $\mathcal{O}(\epsilon^{\min\{1/2,2/d-1/2\}})$ in the Bochner space $L^{p}([s,T];H^{1}(D))$ for any $s\in(0,T]$ and $p\in(1,\frac{1}{\alpha})$ as $\epsilon\to 0$. Thus, the first order approximation $u^{\epsilon}_{1}$ achieves optimal convergence rate of $\mathcal{O}(\epsilon^{1/2})$ when $d\leq 2$ and $\mathcal{O}(\epsilon^{1/6})$ when $d=3$. For the latter case, it is unclear whether or not this rate is optimal. The proof of the result relies on the regularity of time-fractional diffusion problems [14] and the introduction of cut-off functions to handle boundary layers and initial data [10, 17]. When the initial data has a better regularity, e.g., $a(x)\in H^{1}_{0}(D)$, there is no need to introduce a proper cut-off function for $a(x)$.

Furthermore, we present in Section 4 a number of numerical tests for $d=2$ to verify Theorem 3.4. Specifically, we test the cases when $\kappa(x)$ is smooth and discontinuous, and when $\kappa(x)$ admits large deviations. Our numerical tests demonstrate a higher convergence rate, namely, $\mathcal{O}(\epsilon)$, compared to the theoretical result. We also observe that a larger deviation or discontinuity in the coefficient $\kappa(x)$ results in a larger error.

The remainder of this paper is organized as follows. In Section 2 we present the homogenized equation using two-scale asymptotic expansion for problem (1), as well as regularity results on time-fractional diffusion problems. We then introduce auxiliary functions and define the first order approximation in Section 3. The main error estimate is derived therein. To verify our theoretical findings, we present in Section 4 extensive numerical tests with smooth (or nonsmooth) diffusion coefficient and with smooth (or nonsmooth) initial data $a(x)$. Finally, we conclude in Section 5 with a summary and discussion of our results.

2 Two-scale asymptotic expansion

This section is concerned with the two-scale asymptotic expansion of the solution $u^{\epsilon}$ to (1). The approach is standard and can be found, e.g., in [3]. We recall the general procedure for the sake of completeness.

First, analogous to the standard homogenization theory, we denote $y:=x/\epsilon$ as the fast variable, and $x$ is referred as the slow variable. Let L^{2}_{\#}(Y):=\{u\in L^{2}(Y):u\text{ is y-periodic}\} with $Y$ being a unit cell in $\mathbb{R}^{d}$. Similarly, we can define $H^{1}_{\#}(Y)$. Denote $V_{\#}(Y):=\{v\in H^{1}_{\#}(Y):\langle v\rangle=0\}$. The notation $A\lesssim B$ denotes $A\leq CB$ for some constant $C$ independent of the microscale $\epsilon$. Throughout the paper, we follow the Einstein summation convention.

Under the assumption that the fast variable $y$ and the slow variable $x$ are independent when $\epsilon\to 0$, we seek for an asymptotic expansion of the solution $u^{\epsilon}(x,t)$ as follows:

[TABLE]

with the functions $u_{j}(x,y,t)\in H^{1}_{\#}(Y)$ being periodic in the fast variable $y$ with period 1. The leading order term $u_{0}(x,y,t)$ is referred as the homogenized solution and the following terms $\epsilon^{k}u_{k}(x,y,t)$ are the $k^{\rm th}$ order corrector for $u^{\epsilon}(x,t)$ for $k=1,2,\cdots$.

Denote by $A^{\epsilon}$ the second order elliptic operator

[TABLE]

With this notation, we can expand $A^{\epsilon}$ as follows

[TABLE]

where

[TABLE]

Substituting the expansions for $u^{\epsilon}$ and $A^{\epsilon}$ into the differential equation (1), and equating the terms with the same power of $\epsilon$, we get

[TABLE]

Next we examine these equations one by one. Equation (3) is equivalent to seeking $u_{0}(x,y,t)\in H^{1}_{\#}(Y)$, satisfying

[TABLE]

The theory of second order elliptic PDEs implies that $u_{0}(x,y,t)$ is independent of $y$. Consequently, we obtain

[TABLE]

Since there is no micro-scale $\epsilon$ in either the boundary or the initial conditions, those conditions are imposed on the leading order term $u_{0}$ directly, i.e.,

[TABLE]

Thanks to (6), the second equation (4) can be formulated as seeking for $u_{1}\in H^{1}_{\#}(Y)$, satisfying

[TABLE]

Note that $u_{1}$ can be defined alternatively independent of the slow variable $x$. To this end, let $\chi_{j}\in V_{\#}(Y)$ be the solution to the following cell problem:

[TABLE]

The general solution of equation (8) for $u_{1}$ then admits the expression

[TABLE]

For simplicity, we take $\tilde{u}_{1}(x,t)=0$.

Finally, we deal with the third term $u_{2}$. Because of (5), we can seek for $u_{2}\in H^{1}_{\#}(Y)$, such that,

[TABLE]

The solvability condition implies that the right hand side of (10) must have mean zero in $y$ over the unit cell $Y=[0,1]^{d}$, i.e.

[TABLE]

We note that

[TABLE]

for any $F(x,y,t)$ which is periodic with respect to $y$. This can be easily verified using the Divergence Theorem. After integrating (10) over $Y$, the average over the terms starting with ${\partial\over\partial y_{i}}$ disappears and we arrive at

[TABLE]

Substituting the expression for $u_{1}$ into this equation, we obtain the homogenized equation:

[TABLE]

where

[TABLE]

Last, we derive the a priori estimate for the homogenized solution $u_{0}$, which will be utilized below. We obtain from (7) and (11) that

[TABLE]

Then by application of [14, Theorem 2.1], we obtain that $u_{0}\in C([0,T];L^{2}(D))\cap C((0,T];H^{2}(D)\cap H^{1}_{0}(D))$. Furthermore, the following estimate holds:

[TABLE]

Estimate (14) implies that the homogenized solution $u_{0}(x,t)$ is singular at $t=0$. The strength of this singularity depends on the regularity of the initial data $a(x)$. This singularity can disappear when the initial data $a(x)$ has certain regularity. Furthermore, the regularity of the solution $u_{0}$ from time-fractional diffusion problem has only limited regularity over the space domain $D$, c.f. [5].

Remark 2.1**.**

Note that [14, Theorem 2.1] still holds when the permeability coefficient $\kappa^{\varepsilon}$ admits high oscillation. Furthermore, the a priori estimate (14) is stable with respect to the parameter $\varepsilon$.

3 First order approximation estimate

We present in this section the first order approximation to $u^{\varepsilon}$, and then derive its error estimate.

To this end, we will first provide the a priori estimate to the solutions of the cell problem (8). This result can be found, e.g., in [12, Section 8]. For the completeness, we also present the proof:

Lemma 3.1**.**

Let $\chi_{j}$ be the solution to the cell problem (8) for all $j=1,\cdots,d$. Then there holds

[TABLE]

Proof.

Let $e_{j}$ be the $j^{\rm th}$ canonical unit vector in $\mathbb{R}^{d}$. The weak formulation associated to Problem (8) is to seek $\chi_{j}\in V_{\#}(Y)$ such that

[TABLE]

Testing with $v:=\chi_{j}$, we obtain

[TABLE]

Then the boundedness of $\kappa(y)$ and an application of the Hölder’s inequality reveal the assertion (15).

Furthermore, by application of [2, 8], we obtain

[TABLE]

This, together with the weak maximum principle, yields the second assertion (16). ∎

Note that both the exact solution $u^{\epsilon}$ and the homogenized solution $u_{0}$ satisfy a homogeneous Dirichlet boundary condition. However, the first order corrector $\epsilon\chi_{j}({x}/{\epsilon})\frac{\partial u_{0}}{\partial x_{j}}(x,t)$ admits $\epsilon$-oscillation over the global boundary $\partial D$. To account for this fact, we must adapt the corrector near the boundary. We employ an approach developed in [10].

To this end, we introduce the cut-off function $\zeta^{\varepsilon}$ corresponding to $\partial\Omega$. Here $\zeta^{\varepsilon}=1$ on $\partial D$ and $\mbox{supp}(\zeta^{\varepsilon})\subset\{x\in\bar{D}:\;\mbox{dist}(x,\partial D)\leq\varepsilon\}$. With regularity condition $\zeta^{\varepsilon}\in C^{2}\bar{(D)}$, for $0\leq\ell\leq 2$, we have

[TABLE]

In addition, using (18) and the fact that $|\mbox{supp}(\zeta^{\varepsilon})|=O(\varepsilon)$, and also by applying Hölder’s inequality, we have the following estimate for derivatives of the cut-off function $\zeta^{\varepsilon}$ in $L^{q}({D})$ given by

[TABLE]

To avoid higher regularity condition on the initial data $a(x)$, we employ the trick in [17] and introduce another cutoff function $\eta(t;\theta)$ in the time domain $[0,T]$ for any parameter $\theta\in(0,T]$. It is defined as follows:

[TABLE]

Let the first order approximation to $u^{\epsilon}(x,t)$, and its modification be defined by

[TABLE]

Then by definition,

[TABLE]

Denote $u^{\epsilon}=u_{1}^{\epsilon}+R(x,t;\theta)$, where $R(x,t;\theta)$ is the residual between the exact solution $u^{\epsilon}(x,t)$ and its modified first order approximation. Then plugging this expression of $u^{\epsilon}$ into (1) yields

[TABLE]

Collecting the terms and using the homogenized equation (11), the boundary and initial conditions for $u_{0}(x,t)$, we arrive at the following error equation

[TABLE]

Here,

[TABLE]

Remark 3.2** (Regularity of the modified first order approximation (20)).**

Thanks to the condition $a\in L^{2}(D)$, we can obtain $u_{1}^{\epsilon}(x,t;\theta)\in C((0,T];H^{1}_{0}(D))$ for all $\theta\in(0,T]$.

A fundamental ingredient of the homogenization theory will be a proper a priori estimate for the residual $R(x,t;\theta)$. To this end, we first estimate the term $f_{1}$:

Lemma 3.3**.**

For all $s\in(0,T]$, $f_{1}\in C([s,T];H^{-1}(D))$. Moreover, there holds

[TABLE]

Proof.

The following identity holds

[TABLE]

To estimate $\left\|f_{1}\right\|_{H^{-1}(D)}$, we only need to estimate $\left\|T_{1}\right\|_{H^{-1}(D)}$ and $\left\|T_{2}\right\|_{H^{-1}(D)}$. To this end, we can further split the first term into the summation of the two terms

[TABLE]

with

[TABLE]

To estimate $T_{1,1}$, we apply the argument in [3, Section 1.4], together with the definition of the cutoff function $\eta(t;\theta)$, and obtain

[TABLE]

Then we estimate $T_{1,2}$. A direct calculation results in

[TABLE]

Then by the triangle inequality and the chain rule, we deduce

[TABLE]

This, together with the generalized Hölder’s inequality and Sobolev embedding, leads to

[TABLE]

where

[TABLE]

Here, the parameter $\delta\in(0,1/2)$ is an arbitrary constant, and $d$ is the dimension of the domain $D$.

Together with the estimates (15), (19) and (14), by letting $\delta\to 0$, we obtain

[TABLE]

Next we estimate the second term $T_{2}$. The estimate (13) implies

[TABLE]

By exchanging the fractional derivative with respect to time $t$ and the derivative with respect to the space variable $x$, we obtain

[TABLE]

Then by application of (14), we arrive at

[TABLE]

This, together with (22), (23) and an application of the triangle inequality, proves the desired assertion after letting $\theta\to s$. ∎

Finally, we are ready to present the error estimate for the residual $R$ defined in the error equation (21):

Theorem 3.4** (Modified first order approximation).**

For any $p\in[1,\frac{1}{\alpha})$ and $s\in(0,T]$, there exists a unique weak solution $R\in{L}^{p}([s,\,T];{H}^{1}_{0}(D))$ to Problem (21) such that ${\partial^{\alpha}_{t}}R\in{L}^{p}([s,\,T];{H}^{-1}(D))$, satisfying

[TABLE]

Proof.

First, Lemma 3.3 gives

[TABLE]

Therefore, $f_{1}\in L^{p}([s,T];H^{-1}(D))$ for any $p\in[1,\frac{1}{\alpha})$, and there holds

[TABLE]

Meanwhile, following the proof of [14, Theorems 2.1 and 2.2] and [4, Theorem 2.3], we can arrive at

[TABLE]

Consequently, a combination of the previous two estimates proves the desired assertion. ∎

We present in the next result the error estimate with the boundary layer effect:

Corollary 3.5** (First order approximation).**

Let $a(x)\in L^{2}(D)$, $p\in[1,\frac{1}{\alpha})$, $s\in(0,T]$ and $\epsilon$ sufficiently small. Let $u^{\epsilon}$ be the solution to Problem (1), then the following estimate holds

[TABLE]

Proof.

The proof follows from the argument in Theorem 3.4 and the fact that the estimate (23) dominates the convergence rate. ∎

Remark 3.6** (Comparision with homogenization results to parabolic equations with periodic coefficients).**

Due to the boundary layer effect in the bounded domain, we can obtain from [12, 17] that the optimal convergence rate is $\mathcal{O}(\epsilon^{1/2})$. Corollary 3.5 shows that the convergence is $\mathcal{O}(\epsilon^{\min\{1/2,2/d-1/2\}})$. One main restriction to apply similar technique employed in [17] to our current problem of a time-fractional diffusion problem (1) is the fact that there is no product rule for fractional derivative, but there is for the first derivate. For the same reason, one can not derive the pointwise estimate over the time domain for time-fractional diffusion problems.

4 Numerical experiments

In this section, we conduct a series of numerical experiments to demonstrate the performance of the first order corrector introduced in Section 3. Furthermore, we will validate the convergence result presented in Corollary 3.5 corresponding to different permeability fields and fractional power $\alpha$.

Consider the time-fractional diffusion equation (1) in the unit square $D=[0,1]^{2}$ with total time $T=1$ and $\alpha:=0.9$. We will use scalar coefficient $\kappa^{\epsilon}(x_{1},x_{2})$ in the following numerical tests. The smooth initial data tested in Sections 4.1 and 4.2 is

[TABLE]

We refer to Figure 1 for an illustration.

In Sections 4.1 and 4.2, we are concerned with the convergence rate of the first order approximation $U_{1}^{\epsilon}$ for diffusion coefficients $\kappa^{\epsilon}(x_{1},x_{2})$ of different regularities. To this end, we test two kinds of permeability fields, namely, the smooth and nonsmooth permeability fields in these two sections, respectively. Since Corollary 3.5 is also valid for rough initial data $a(x)\in L^{2}(D)$, we present the convergence history of the first order approximation $U_{1}^{\epsilon}$ with a rough initial data in Section 4.3.

4.1 Numerical tests with smooth permeability fields

To define the smooth permeability field, we take

[TABLE]

as the periodic smooth function defined over the unit square $Y$. Here, $\{\cdot\}$ means taking the fractional part.

Recall that $\kappa^{\epsilon}(x_{1},x_{2})$ is a periodic function with period $\epsilon$. Its one cell after stretching over one unit cell is $\kappa(y_{1},y_{2})$, i.e., $\kappa^{\epsilon}(x_{1},x_{2}):=\kappa(x_{1}/\epsilon,x_{2}/\epsilon)$. See Figure 2 for an illustration. Note that the contrast for this permeability is $\frac{11}{9}$. The main aim of this section is to investigate the convergence rate of the first order corrector $U_{1}^{\epsilon}$. We will present the absolute and relative errors between $U_{1}^{\epsilon}$ and $u^{\epsilon}$ in $L^{2}$-norm and $H^{1}$-norm.

Let $\mathcal{T}_{h}$ be a decomposition of the domain $D$ into non-overlapping shape-regular rectangular elements with maximal mesh size $h:=2^{-9}$. Let $V_{h}$ be the conforming piecewise affine finite element space associated with the partition $\mathcal{T}_{h}$:

[TABLE]

where $\mathcal{Q}_{1}(T)$ denotes the space of affine polynomials on each element $T\in\mathcal{T}_{h}$. Let $V_{h}^{0}:=V_{h}\cap H^{1}_{0}(D)$. We discretize the time interval $[0,1]$ with a time step $\Delta t:=1/100$. Let $t_{n}=(n-1)\Delta t$ for $n=1,2,\cdots,N$ with $N:=\Delta t^{-1}+1$.

We adopt one popular scheme [9] to discretize the time variable $t$ in (1) and (13), and apply the conforming Galerkin method to approximate the exact solution $u^{\epsilon}$. We denote its approximation at $t=t_{k+1}$ by $u^{\epsilon,k+1}_{h}$ for $k=0,1,\cdots,N-1$.

To this end, we seek for $u^{\epsilon,k+1}_{h}\in V_{h}^{0}$ for $k=0,1,\cdots,N-1$, satisfying

[TABLE]

Here, the parameters $b_{j}$ are given by

[TABLE]

The numerical solutions $u^{\epsilon,k}_{h}$ for $k=11,51\text{ and }101$ are depicted in Figure 3.

We denote $U^{\epsilon,k}_{1,h}$ (or $u_{0,h}^{k}$) for $k=1,2\cdots,N$ as the numerical approximation to $U^{\epsilon}_{1}(t)$ (or $u_{0}(t)$) for $t=0,\Delta t,\cdots,1$. In order to obtain the first order approximation $U^{\epsilon,k}_{1,h}$, we first solve the cell problem (8). To this end, we first divide the computational domain $Y$ into non-overlapping shape-regular rectangular elements with a maximal mesh size $h:=2^{-6}$. Then we solve the cell problem (8) with continuous piece-wise bilinear Lagrange Finite Element Method using the conforming Galerkin formulation. We plot the two cell solutions $\chi_{1}$ and $\chi_{2}$ in Figure 4 with the smooth permeability field in Figure 2.

Utilizing the solutions to the cell problem (8), i.e., $\chi_{1}$ and $\chi_{2}$, we can obtain the effective coefficient $\kappa^{*}$ from (12), and then solve for $u_{0}$ by (13). Note that there is no microscale oscillation in the effective coefficient $\kappa^{*}$, thus the solution $u_{0}$ can be solved in a much coarser mesh compared to the mesh $\mathcal{T}_{h}$ associated to the original problem (1). To simplify our notations, we adopt the same mesh $\mathcal{T}_{h}$ and finite element space $V_{h}^{0}$ as before, and utilize the conforming Galerkin formulation to solve for $u_{0,h}^{k+1}$ for $k=0,1,\cdots,N-1$.

We present the fine scale approximate solutions $u^{k}_{0,h}$ for $k=11,51\text{ and }101$ in Figure 5.

Finally, the first order approximation $U^{\epsilon}_{1}$ can be estimated using formula (20). We present the graphs of the approximate solutions $U^{\epsilon,k}_{1,h}$ for $k=11,51\text{ and }101$ in Figure 6.

We present the absolute error and relative error in $L^{2}$ norm and $H^{1}$ norm in Table 1. The first column displays the discrete time steps at which we calculate the error. The next two columns display the absolute error and relative error under the $L^{2}$ norm between the fine-scale solution $u_{h}^{\epsilon}$ and the first order approximation $U^{\epsilon}_{1,h}$. The last two columns display the absolute error and relative error under the $H^{1}$ norm between the fine-scale solution $u_{h}^{\epsilon}$ and the first order approximation $U^{\epsilon}_{1,h}$.

Furthermore, one numerical experiment is conducted with larger variation in the coefficient $\kappa$ compared to (25) to see the influence of the variation on the accuracy of homogenization. In this experiment, we take the same initial data $a(x_{1},x_{2})$ and the period $\epsilon:=\frac{1}{8}$. We set

[TABLE]

Note that the variation in this coefficient $\kappa^{\varepsilon}$ is much larger than that defined in (25). The convergence history of the first order approximation with $\kappa(y_{1},y_{2})$ in (26) is presented in Table 2. Compared with the results in Table 1, a larger variation in the diffusion coefficient $\kappa(y_{1},y_{2})$ results in a larger error in the first order approximation.

In order to verify the convergence rate of the first approximation $U_{1}^{\epsilon}$, we test two different values of the parameter $\epsilon$ being $\frac{1}{16}$ and $\frac{1}{32}$, respectively. Their corresponding convergence histories are shown in Tables 3 and 4.

One can calculate directly from Tables 1, 3 and 4 that the first order approximation maintains a convergence rate of $\mathcal{O}(\epsilon^{0.9623})$, $\mathcal{O}(\epsilon^{0.9657})$ and $\mathcal{O}(\epsilon^{0.9661})$ for $t=0.1$, $t=0.5$ and $t=1$, respectively.

We also test different values of the parameter $\alpha\in(0,1)$ and they all exhibit a similar convergence rate, as proved in Corollary 3.5. For the brevity of presentation, we will not present these results.

4.2 Numerical tests with non-smooth permeability fields

Even though the theoretical result presented in Corollary 3.5 is proved under the assumption that the diffusion coefficient $\kappa(y_{1},y_{2})$ is sufficiently smooth, we investigate in this section how well the first order approximation $U_{1}^{\epsilon}(x,t)$ performs when the coefficient $\kappa(y_{1},y_{2})$ is rough and can admit large variation.

Firstly, we define two rough permeability fields $\kappa(x_{1},x_{2})$ with different variations. Let $U:=[\frac{1}{5},\frac{4}{5}]^{2}$ be a rectangle. Recall that $Y=[0,1]^{2}$ is a unit square. We set the variation in the first permeability field to be 1.1, which is defined by

[TABLE]

See Figure 7 for an illustration.

We set the variation of the second rough permeability field to 2 and define

[TABLE]

We present the graphs of such $\kappa(y_{1},y_{2})$ and $\kappa^{\epsilon}(x_{1},x_{2})$ in Figure 8.

Let the parameter $\epsilon:=1/8$. We present the convergence histories of the first order approximation $U_{1}^{\epsilon}(x,t)$ with those two nonsmooth permeability fields in (27) and (28) in Tables 5 and 6. One can observe that the former outperforms the latter.

4.3 Numerical tests with rough initial data

Now we study the convergence rate of the first order approximation $U^{\epsilon}_{1}(x,t)$ when the initial data $a(x)$ is nonsmooth. To this aim, we take a rough initial data $a(x_{1},x_{2})$ defined by

[TABLE]

which is depicted in Figure 9.

To emphasize the effect of a rough initial data on the convergence rate of the first order approximation, we take the smooth permeability field $\kappa(x_{1},x_{2})$ as defined in (25). We adopt the same numerical scheme to calculate the numerical solutions as in Section 4.1. Like before, we test $\epsilon=\frac{1}{8},\frac{1}{16}$ and $\frac{1}{32}$ to validate the convergence rate proved in Corollary 3.5. Their corresponding convergence histories of the first order approximation to Problem (1) are presented in Tables 7, 8 and 9.

Comparing the results presented in Table 7 with that in Table 1, one can observe that the latter admits a first order approximation with a much smaller error as expected, since a rough initial data produces singularity in the solution when the time $t$ is small.

One can calculate from Tables 7, 8 and 9 that the first order approximation has a convergence rate of $\mathcal{O}(\epsilon^{0.9523})$, $\mathcal{O}(\epsilon^{0.9502})$ and $\mathcal{O}(\epsilon^{0.9499})$ for $t=0.1$, $t=0.5$ and $t=1$, respectively.

5 Conclusion

This paper is concerned with constructing a homogenization theorem for the time-fractional diffusion equations with a homogeneous Dirichlet boundary condition and an inhomogeneous initial data $a(x)\in L^{2}(D)$ in a bounded convex polyhedral domain $D$. Under the assumption that the diffusion coefficient $\kappa^{\epsilon}(x)$ is smooth and periodic with a period $\epsilon>0$ being sufficiently small, we proved that its first order approximation has a convergence rate of $\mathcal{O}(\epsilon^{1/2})$ when the dimension $d\leq 2$ and $\mathcal{O}(\epsilon^{1/6})$ when $d=3$. Several numerical tests are presented to show the performance of the first order approximation. Future studies include deriving the homogenization results for the case of existence of extra time scale in the diffusion coefficient $\kappa^{\epsilon,\eta}(x,t)$ with $\epsilon$ and $\eta$ being small parameters.

Acknowledgement

GL acknowledges the support from the Royal Society through a Newton International Fellowship.