Social contacts, epidemic spreading and health system. Mathematical modeling and applications to COVID-19 infection

1.   Introduction

Fractional diffusion involves phenomena that have spatial and temporal correlations ^[1,2]. Anomalous diffusion through fractional equations is associated with super-statistics and can be linked to a generalized random walk ^[3]. The phenomenon of anomalous diffusion has received widespread attention in the fields of natural sciences, engineering, technology, and mathematics ^[4,5,6]. The fractional diffusion equations which serve as models for describing this phenomenon are of utmost importance ^{[7,8,9,10,11]}. Numerous publications have been dedicated to this field so far (e.g., Sakamoto and Yamamoto ^[12]). In contrast to classical parabolic equations, the time fractional diffusion equations replace the traditional local partial derivative $\partial_{t}$ with the nonlocal fractional derivative $\partial_{t}^{\alpha}$. The fractional equations are highly regarded in mathematical physics and present distinct properties that challenge conventional differential equations. Nevertheless, some properties, such as the maximum principle, remain valuable in our research. This paper plans to describe the behavior of time fractional diffusion equations.

In this paper we consider the following initial boundary value problem (IBVP) of a time fractional diffusion equation with a period structure

where $0 < \alpha < 1$, $T > 0$, $\Omega\subset \mathbf{R}^{d}$ is bounded domain with $C^{2}-$class boundary $\partial\Omega$, $\varepsilon > 0$ is a scale parameter, $B^{\varepsilon}(x) = B(\frac{x}{\varepsilon})$ is a diffusion matrix which satisfies some appropriate conditions, $B(y)$ is periodic, and $f^{\varepsilon}(x, t)$ and $u^{\varepsilon}_{0}(x)$ are the source function and the initial function, respectively.

The existence and uniqueness of solutions to the initial boundary value problem (1.1) have been investigated widely. Sakamoto and Yamamoto ^[12] derived a kind of solution in terms of the Fourier series; Kubicam, Ryszewska and Yamamoto ^[13] gave the variational formulation; Hu and Li ^[14] gave the formally homogenized equation by the multiple scale expansion as $\varepsilon\rightarrow 0^{+}$ and Kawamoto, Machida and Yamamoto ^[15] gave the homogenized equation

where $B^{0}$ is the homogenized coefficient matrix, and then proved the precise homogenization theorem; they also discussed the inverse problem between different structures in the one dimensional case and in the layered material case where $B^{\varepsilon}(x)$ is a diagonal matrix with an unknown element when $f = 0$. The aim of this paper is to generalize this result from the case with only one unknown element to the case with multiple unknown elements.

The rest of this paper is organized as follows. In Section $2$, we introduce some necessary tools, including the well-posedness and homogenization of fractional diffusion equations with oscillating diffusion matrix, the eigenvalue problem, and the Mittag-Leffler function. In Section 3 and Section 4 we state main results and prove them. In Section 5, we draw concluding remarks.

2.   Preliminaries

In this section, we state some basic tools to investigate the inverse problems of the initial boundary value problem (1.1) and its homogenized equation (1.2), including the well-posedness, homogenization theory, the eigenvalue problem of the corresponding elliptic operator, and the Mittag-Leffler function, see ^[13,16].

2.1.   Weak solution for fractional diffusion equation

We recall the Riemann-Liouville fractional integral operator

then the domain $\mathcal{D}(J^{\alpha}) = L^{2}(0, T)$ and the range $\mathcal{R}(J^{\alpha}) = H_{\alpha}(0, T)$ with

where $H^{\alpha}(0, T)$ is the Sobolev space. Moreover, $J^{\alpha}: L^{2}(0, T)\rightarrow H_{\alpha}(0, T)$ is a homeomorphism with

Therefore, the general fractional derivative of the Caputo type is defined by

Obviously, $\partial^{\alpha}_{t}$ is also a homeomorphism.

We now consider the initial boundary value problem

where $u_{0}(x)\in L^{2}(\Omega)$, $f(x, t)\in L^{2}(0, T; H^{-1}(\Omega))$, and the matrix $B(x) = (b_{ij}(x))_{d\times d}$ satisfies that

for $0 < \nu < \mu$. From ^[13,15], we know that there exists a weak solution $u\in L^{2}(0, T; H_{0}^{1}(\Omega))$ satisfying $u-u_{0}\in H_{\alpha}(0, T; H^{-1}(\Omega))$, and

for a.e. $t\in(0, T)$ and $\forall\phi\in H_{0}^{1}(\Omega).$

2.2.   Homogenization for fractional diffusion equation

For $Y = (0, l_{1})\times\dots\times(0, l_{d})$, we say that a function $f(x)$ is Y-periodic if

Theorem 2.1. ^[15] For $B^{\varepsilon}(x) = B(\frac{x}{\varepsilon})$, assume that $B(y) = (b_{ij}(y))_{d\times d}$ is Y-periodic and satisfies Eq (2.6), $u_{0}^{\varepsilon}\in L^{2}(\Omega)$, and $f^{\varepsilon}\in L^{2}(0, T; H^{-1}(\Omega))$. If

and $u^{\varepsilon}$ is the weak solution of IBVP (1.1), then

where $u^{0}$ is the weak solution of the homogenized problem (1.2), and $B^{0}$ is the homogenized coefficient matrix. Furthermore, for the layered material, that is, $B(y)$ is a diagonal matrix

then

where $\mathcal{M}_{\Omega}(b) = \frac{1}{|\Omega|}\int_{\Omega}b(x)dx.$

2.3.   Eigenvalue problem

For the diagonal matrix $B_{p} =$diag$\{p_1, \cdots, p_d\}$, $p_i$ are constants and $\nu\le p_i \le \mu$, denote a vector $p = (p_1, p_2, \cdots, p_d)\in \mathbb{R}^{d}$ and an operator $\mathcal{B}_{p}(\cdot) = -$div$(B_p\nabla\cdot)$, we consider an eigenvalue problem of the operator $\mathcal{B}_{p}$ on $\Omega = \prod_{i = 1}^{d}(0, \delta_i)$.

According to the domain $\Omega$, we consider the sub-eigenvalue problems

Then, we can verify that $\varphi(x) = \prod_{i = 1}^d\varphi_i(x_i)$ is a solution of eigenvalue problem (2.11) with $\lambda = \sum_{i = 1}^d p_i\sigma_i$, i.e., $\mathcal{B}_{p}\varphi = (\sum_{i = 1}^d p_i\sigma_i)\varphi.$ Denote by $\varphi_i^{k_i}$ the $k_i-$th simple eigenvalue of the $i-$th sub-eigenvalue problem, that is, $-\frac{d^{2}}{dx^{2}_{i}}\varphi_i^{k_i}(x_i) = \sigma_i^{k_i}\varphi_i^{k_i}(x_i).$ It is known that $\varphi_i(x_i)$ are the sine functions. Since the eigenfunctions $\{\varphi_i^k\}_{k = 1}^{\infty}$ are an orthonormal basis of $L^2(0, \delta_i)$, so $\{\prod_{i = 1}^d\varphi_i^{k_i}(x_i)\}_{k_1, \dots, k_d \in \mathbb{N}}$ is an orthonormal basis of $L^2(\Omega)$. Then, we can prove that the all eigenvalues $\lambda$ of $\mathcal{B}_{p}$ have following form:

In fact, for $\mathcal{B}_{p}\varphi = \lambda\varphi$, $\varphi \neq 0$, there exists $\prod_{i = 1}^d\varphi_i^{k_i}$ such that $(\varphi, \prod_{i = 1}^d\varphi_i^{k_i})_{L^2(\Omega)}\neq 0.$ Taking inner product of $L^2(\Omega)$ with respect to $\prod_{i = 1}^d\varphi_i^{k_i}$ on both sides of equation (2.11) and integration by parts, we can complete the proof.

For example, for $d = 3, p_i = 1, i = 1, \cdots, d$, we have $\mathcal{B}_{p} = -\Delta$. Taking $\delta_1 = \delta_2 = \delta_3,$ we have $0 < \sigma_i^1 < \sigma_i^2 < \cdots\to+\infty,$ $i = 1, 2, 3,$ $\sigma_1^k = \sigma_2^k = \sigma_3^k,$ and $\varphi_1^k = \varphi_2^k = \varphi_3^k,$ $k\in\mathbb{N}$. Then we can write the eigenvalue of Eq (2.11) as

and the corresponding eigenfunctions are

Thus, we know that some eigenvalues of problem (2.11) have more than one geometric multiplicity, which is different from the eigenvalues of problem (2.12) such that all eigenvalues are simple.

Returning to the general operator $\mathcal{B}_{p}$ on a bounded domain $\Omega\subset \mathbb{R}^{d}$, we rearrange the eigenvalues of $\mathcal{B}_{p}$ without multiplicity, $0 < \lambda_1 < \lambda_2 < \cdots\to+\infty,$ and rearrange the eigenvalues $\{\sigma_i^{k}\}_{k = 1}^{\infty}(i = 1, \cdots, d)$ such that

where $m_n$ is the multiplicity of the eigenvalue $\lambda_n = \lambda_n^{1} = \lambda_n^{2} = \cdots = \lambda_n^{m_n}$. Note that $\lambda_1$ is simple, i.e. $m_1 = 1$ and $\sigma_i^{1_1} < \sigma_i^{n_j}, \; n > 1, \; j = 1, \cdots, m_n.$ Set

where $\varphi_i^{n_j}(x_{i})$ is the eigenfunction corresponding to eigenvalue $\sigma_i^{n_j}$, and $\{\varphi^{n_j}\}_{j = 1}^{d_n}$ is the orthonormal basis of $ker(\mathcal{B}_{p}-\lambda_n I)$.

We introduce a projection operator $P_n: L^2(\Omega)\to ker(\mathcal{B}_{p}-\lambda_n I)$ such that

is an eigenprojection. We note that the eigenfunctions of $-\Delta$ and $\mathcal{B}_p$ are identical indeed, but their eigenvalues are not identical.

2.4.   Mittag-Lefler function

From ^[17], the Mittag-Leffler function is defined by

which is an entire function in the complex plane.

When $\alpha = \beta = 1$, $E_{\alpha, \beta}$ is precisely an exponential function. What is more, we have the asymptotic expansion and estimate

where $K\in\mathbb{N}$, $0 < \alpha < 2$, and $\frac{\pi\alpha}{2} < \theta <$min$\{\pi, \alpha\pi\}$.

3.   Statement of main results

We now consider problem (1.1) on the domain $\Omega = \prod_{i = 1}^d(0, \delta_i), \; \delta_i > 0, d < 4$. In order to state the main results, we first give a definition.

Definition 13.1. For the matrices $A = (a_{ij})_{N\times M}$ and $B = (b_{ij})_{N\times M}$, we say $A\ge B$ if

We say $A > B$ if $A\ge B$ and there exists an index $i$ or $j$ such that $a_{ij} > b_{ij}.$

We consider the IBVP with a periodic structure

with unknown initial function $u^{\varepsilon}_{0, p}$ and unknown diffusion matrix $B_p^{\varepsilon}$ with layer structure

satisfying

Due to Theorem 2.1, we get $u_p^{\varepsilon}\rightarrow u^0_p$ in $L^2(0, T; L^2(\Omega))$, where $u^0_p(x)$ is a weak solution of the homogenized equation

where $u_{0, p}(x)$ is the $L^2(\Omega)$ limit of $u^{\varepsilon}_{0, p}(x)$, $B^{0}_p =$diag$\{p_1^0, p_2^0, \dots, p_d^0\}$ and

satisfying $\nu\le p_i^0\le\mu, \; i = 1, \dots, d.$ By Eq (3.5) and simple calculation, we have

with $C = C(\nu, \mu, l_1) > 0$. Moreover, if

we have

with $C = C(\nu, \mu, l_1)\geq1$. These can be seen in the proof of [15, Lemma 3.11].

For simplicity of notation, we set $p = (p_1, \cdots, p_d)\in\mathbb{R}^d_+$. We first consider several inverse problems of determining the diffusion matrix of the following IBVP:

where $A_p = \text{diag}\{p_1, \cdots, p_d\}$ and $\nu \le p_1, \cdots, p_d\le\mu$.

Inverse problem Ⅰ: Let $x_0\in\Omega\text{, }t_0\in(0, T)$. We will determine the diffusion matrix $A_p$ by the single data point $u_p(x_0, t_0)$ of problem (3.9).

Theorem 3.1. Let $u_{0, p} = u_{0, q} = u_0\in H^2(\Omega)\cap H^1_0(\Omega)$ and $u_{p}(x, t)$ be a solution of problem (3.9). Assume $p\ge q$ and

Then there exists a constant $C(\nu, \mu) > 0$ such that

Inverse problem Ⅱ: Let $\omega\subset\Omega\text{, }I\subset(0, T)$. We will determine the diffusion matrix $A_p$ by the data $\int_{I}\int_{\omega}u_p^{}(x, t)dxdt$ of problem (3.9). Note that the measurement data is an integral expression. Thus, it is more useful for applications.

Theorem 3.2. Let $u_{0, p} = u_{0, q} = u_0\in H^2(\Omega)\cap H^1_0(\Omega)$ and $u_{p}(x, t)$ be a solution of problem (3.9). Assume $p\ge q$ and

Then there exists a constant $C(\nu, \mu) > 0$ such that

Inverse problem Ⅲ: Let $x_0\in\Omega\text{, }t_1\in(0, T)$. We determine the diffusion matrix $A_p$ by the time trice data $u_p(x_0, t), \; 0 < t < t_1$ of problem (3.9). We can only get uniqueness here.

Theorem 3.3. Let $u_{p}(x, t)$ be a solution of problem (3.9) and $u_{0, p}, u_{0, q}\in H^2(\Omega)\cap H^1_0(\Omega)$ such that

If $u_p(x_0, t) = u_q(x_0, t), \; 0 < t < t_1$, then $A_p = A_q.$

Remark 3.1. Let $\Omega = (0, \pi)\times(0, \pi)$. The eigenfunction corresponding to $\lambda_1 = 2$ is $\varphi_1 = \frac{2}{\pi}sinxsiny$ and the eigenfunction corresponding to $\lambda_2 = 5$ is $\varphi_2 = \frac{2}{\pi}sin2xsiny$ and $\varphi_3 = \frac{2}{\pi}sinxsin2y$. We take $u_{0, p} = \varphi_1$ and $u_{0, q} = \varphi_2$. We have

and

If $p = (4, 1), q = (1, 1)$, and $(x_0, y_0) = (\frac{\pi}{3}, \frac{\pi}{2})$, we have $u_p(x_0, y_0, t) = u_q(x_0, y_0, t), \; t > 0$, but $A_p\neq A_q$. Thus, condition (3.14) is necessary.

We will present the proof of these theorems later in Section 4.

Following Theorems 3.1–3.3 and the estimates (3.8), we can immediately obtain the following corollaries of the inverse problem determining the diffusion matrix between different structures to problem (3.1) and problem (3.4). First, we present the inverse problem of determining the period coefficient matrix by the homogenized data.

Corollary 3.1. Let $u_{p}^{0}(x, t)$ be a solution of problem (3.4) and $u_{0, p} = u_{0, q} = u_0\in H^2(\Omega)\cap H^1_0(\Omega)$. Under the condition (3.7) and

Then there exists a constant $C(\nu, \mu, l_1) > 0$ such that

We see that the condition for $u_0$ is from Condition 3.10 and $p(y), q(y)$ are vector-valued functions over $(0, l_1)$ portraying the period structure. Further, we must guarantee $u^{\epsilon}_{0, p}$ and $u^{\epsilon}_{0, p}$ have the same limit $u_0$. Similarly, the following result also follows.

Corollary 3.2. Let $u_{0, p} = u_{0, q} = u_0\in H^2(\Omega)\cap H^1_0(\Omega)$, $u_{p}^{0}(x, t)$ be a solution of problem (3.4) and $u_{p}(x, t)$ be a solution of problem (3.9). Under condition (3.7) and

Then there exists a constant $C(\nu, \mu, l_1) > 0$ such that

Corollary 3.3. Let $u_{p}^{0}(x, t)$ be a solution of problem (3.4) and $u_{0, p}, u_{0, q}\in H^2(\Omega)\cap H^1_0(\Omega)$ such that

Then if $u^0_p(x_0, t) = u^0_q(x_0, t), \; 0 < t < t_1$, we have $B_p(y) = B_q(y)\; a.e.\; y\in Y.$

Limited by our approach, as in Theorem 3.3, we can only obtain the uniqueness of this inverse problem. We can also use the periodic structure data to determine the homogenized coefficient matrix as the following result of asymptotic stability.

Corollary 3.4. Let $u_{0, p} = u_{0, q} = u_0\in H^2(\Omega)\cap H^1_0(\Omega)$, $u_{p}^{\varepsilon}(x, t)$ be a solution of problem (3.1) and $u_{p}^{0}(x, t)$ be a solution of problem (3.4). Under the condition (3.7) and

Then there exists a constant $C(\nu, \mu, l_1) > 0$ such that

where $\theta(\varepsilon)\rightarrow 0$ as $\varepsilon\rightarrow 0$.

4.   Proof of Theorems 3.1–3.3

In this section, we give the proof of Theorems 3.1–3.3.

Proof of Theorem 3.1. We split the proof into the following three steps, and prove separately for each step.

Step 1. We first prove that $p\ge q$ implies that $u_p(x, t)\ge u_q(x, t)$. Set $y = u_p-u_q$, then

Denote $v = \sum_{i = 1}^d(p_i-q_i)\partial_i^2u_q =$div$(A_{p-q}\nabla u_q): = -\mathcal{A}_{p-q}u_q$. Since

we get

Thus, by [12, Theorem 2.1], we know that $v$ is a weak solution of the following problem

By [18, Theorem 2.1] and $\sum_{i = 1}^d(p_i-q_i)\partial_i^2u_0\ge 0$, we have $v\ge 0$. Applying [18, Theorem 2.1] to Eq (4.1), we get $y\ge 0$, i.e., $u_p\ge u_q$.

Step 2. For $p = (p_1, \cdots, p_d)$, we prove the analyticity of $u_p(x_0, t_0)$ with respect to every $p_i > 0$. Observe that

hence

By the Sobolev embedding $H^2(\Omega)\hookrightarrow C(\overline{\Omega})$, we have

Thus, we get the convergent series

and since $E_{\alpha, 1}(z)$ is holomorphic in the complex plane, we see that $h(p) = u_p(x_0, t_0)$ is analytic with respect to $p_i, \; i = 1, \dots, d$.

Step 3. We prove that $p > q$ means $h(p) > h(q)$. First, $\partial_i^2u_0\ge 0, \; i = 1, \cdots, d$, so $\Delta u_0\ge0$. Since $u_0\in H^2(\Omega)\subset C(\overline{\Omega})$ and $u_0\in H^1_0(\Omega)$, by the strong maximum principle ^[20], we have $u_0\le 0$. On the basis of [19, Theorem 9], we know that $u_p(x_0, t_0) < 0$ for all $p\in\mathbb{R}^d_+$.

If there exist $p^0, q^0$ such that $p^0 > q^0$ and $u_{p^0}(x_0, t_0) = u_{q^0}(x_0, t_0)$, then there is $i\in\{1, \cdots, d\}$ such that $p_i^0 > q_i^0$. Therefore, when $q_i^0 < s < p_i^0$,

Since $h(p)$ is analytic with respect to $p_i$, we have

Moreover,

where we use the fact that

and the series $\sum_{n = 1}^{\infty}\sum_{j = 1}^{m_n}\sigma_1^{-n_j}$ converges. Passing to the limit as $s\rightarrow \infty$ in Eq (4.4), we have $u_{p^0}(x_0, t_0) = 0$, a contradiction. Therefore, $p > q$ means that $h(p) > h(q)$.

Step 4. By the last step, we have $\partial_i h(p) > 0$ for all $p\in\mathbb{R}^d_+, i = 1, \cdots, d.$ Set $g(t) = h(q+t(p-q))$, then $g'(t) = \nabla h(q+t(p-q))\cdot(p-q)$. By the mean value theorem we get

where

This ends the proof of Theorem 3.1. □

Proof of Theorem 3.2. Referring to Theorem 3.8(ⅱ) in ^[15], we can similarly verify that the function $H(p) = \int_{I}\int_{\omega}u_p^{}(x, t)dxdt$ satisfies $\partial_i H(p) > 0$ for all $p\in\mathbb{R}^d_+, \; i = 1, \cdots, d.$ Thus, we can complete the proof similarly to that of Theorem 3.1. □

Proof of Theorem 3.3. From $u_p(x_0, t) = u_q(x_0, t)$ for $0 < t < t_1$ and the analyticity of $u_p(x, t)$ with respect to $t$, we have $u_p(x_0, t) = u_q(x_0, t)$ for $t > 0.$ Since

and by the asymptotic expansion (2.13), we have that

holds for $K\in\mathbb{Z}^+$. We equate the coefficients of $t^{-\alpha k}$, which yields

If $q > p$, we have

We observe that

Taking $k\rightarrow \infty$ in Eq (4.5), we get

this yields a contradiction. Thus, $q\leq p$. Similarly, we can prove $p\leq q$. This ends the proof of Theorem 3.3. □

5.   Conclusions

We should mention that our work is done under some constraints. We require the diffusion coefficient matrix to be diagonal and the area considered to be a rectangular area. This is because our method relies on the expansion of the eigenfunction, and only under these constraints can we clarify our eigenfunction, which is very advantageous for our proof. In addition, we only consider the case of layered matter, that is, the diffusion matrix only depends on one variable. Only in this way can we obtain formula $(3.8)$. In order to break through these limitations, we believe that we can only seek other more difficult methods. Besides, we can also consider other more general problems. For example, we can consider the inverse problem of the variable-order time-fractional equation ^[21] in the current frame, but we do not discuss these problem in this paper. Moreover we also consider the inverse problem of determining the variable order and diffusion matrix simultaneously, such as in article ^[22], which investigates this problem in one space dimension. We also hope to expand their results to the situation of high-dimensional situations in the future.

Acknowledgments

We would like to thank the Editor and Referees for valuable comments and contributions. The research is supported by the National Natural Science Foundation of China (No. 12171442). Feiyang Peng carried out the inverse problem and the homogenization theory of fractional diffusion equations, and Yanbin Tang carried out the reaction diffusion equations and the perturbation theory of partial differential equations. All authors carried out the proofs and conceived the study. All authors read and approved the final manuscript.

Conflict of interest

The authors declare that there is no conflict of interest.