Determination of a Nonlinear Coefficient in a Time-Fractional Diffusion Equation

Abstract

In this paper, we study direct and inverse problems for a nonlinear time fractional diffusion equation. We prove that the direct problem has a unique weak solution and the solution depends continuously on the coefficient. Then we show that the inverse problem has a quasi-solution. The direct problem is solved by the method of lines using an operator approach. A quasi-Newton optimization method is used for the numerical solution to the inverse problem. The Tikhonov regularization is used to overcome the ill-posedness of the inverse problem. Numerical examples with noise-free and noisy data illustrate the applicability and accuracy of the proposed method to some extent.

1. Introduction

In recent years, different approaches have been provided for the modeling of diffusion in fractal geometry [1,2,3,4,5,6,7,8,9,10]. Those approaches are mainly based on fractional differential equations. In this article we also adopted a similar approach. For the sake of completeness, we give a brief background. One-dimensional mass transport due to diffusion is formulated by

$$\frac{\partial u}{\partial t}=-\frac{\partial }{\partial x}\left(J_u\right),$$

(1)

where $J_u$ is the diffusive mass flux. For Fickian diffusion, $J_u$ is defined by the following equation:

$$J_u=-D^{\left(m\right)}\frac{\partial u}{\partial x},$$

(2)

where $D^{\left(m\right)}$ is the diffusion coefficient. Defining the new non-dimensional variables $X=\frac{x}{x_0}$, $\tau =\frac{t}{t_0}$, $U=\frac{u}{u_0}$, where $x_0$, $t_0$ and $u_0$ are the characteristic scales, the Equations (1) and (2) can be given in the non-dimensional form as follows:

$$\frac{1}{t_0}\frac{\partial U}{\partial \tau }=\frac{1}{x_0^2}\frac{\partial }{\partial X}\left(D^{\left(m\right)}\frac{\partial U}{\partial X}\right).$$

(3)

If $x_0=t_0$, then the Equations (1) and (2) preserve their original form. However, many experiments with fractal objects show that this correlation does not hold. In this case, it is proved that the-mean square displacement of a random walker $<x^2>\sim t^{\frac{2}{2+\theta }}$, where $\theta$ is the index of the anomalous diffusion. Thus, in our case we have the correlation $<x_0^2>\sim t_0^{\frac{2}{2+\theta }}$ and it is clear that there are many expressions for the mass flux that correspond to that correlation. For example, the diffusion coefficient may be defined as $D^{\left(m\right)}\left(x\right)=D_f{x}^{-\theta }$, where the constant $D_f$ is the effective diffusion coefficient. Therefore, we obtain the following diffusion equation:

$$\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(D_f{x}^{-\theta }\frac{\partial u}{\partial x}\right).$$

(4)

As another example, the mass flux may be taken proportional to the fractional derivative of concentration with respect to spatial coordinate of the order $\theta +1$. In this case, the obtained mass flux is doubtful since the order of the corresponding diffusion equation is greater than 2. For this reason, the following expression for the mass flux is considered:

$$J_u=D_f{\partial }_t^{1-\beta }\left({\partial }_x^{\gamma }u\right),\beta >0,\gamma <1,$$

(5)

where $\beta$ and $\gamma$ are the order of the temporal and spatial fractional derivatives, respectively. In the expression (5), ${\partial }_t^{\beta }$ and ${\partial }_x^{\gamma }$ are spatial and temporal Caputo fractional derivatives and are defined as the following:

$${\partial }_x^{\gamma }u(x,t)=\frac{{\partial }^{\gamma }u(x,t)}{\partial x^{\gamma }}=\frac{1}{\Gamma (1-\gamma )}{\int }_0^x{(x-\xi )}^{-\gamma }\frac{\partial u(\xi ,t)}{\partial \xi }d\xi ,$$

$${\partial }_t^{\beta }u(x,t)=\frac{{\partial }^{\beta }u(x,t)}{\partial t^{\beta }}=\frac{1}{\Gamma (1-\beta )}{\int }_0^t{(t-\xi )}^{-\beta }\frac{\partial u(x,\xi )}{\partial \xi }d\xi ,$$

where $\Gamma (·)$ is the Gamma function. We also note that there is another frequently used fractional derivative, so called the Riemann-Liouville fractional derivative, which is defined by

$${}^R{\partial }_t^{\beta }u(x,t)=\frac{1}{\Gamma (1-\beta )}\frac{\partial }{\partial t}{\int }_0^t{(t-\xi )}^{-\beta }u(x,\xi )d\xi .$$

These two fractional derivatives agree when the initial condition is zero. Kilbas et al. [11] and Podlubny [12] can be referred for further properties of the Caputo and Riemann-Liouville fractional derivatives. Since the initial condition is zero in our problem, any result found in the literature for one of these holds for the other one. Substituting the temporal and spatial fractional Caputo derivatives in the expression (5), we obtain:

$$\frac{\partial u}{\partial t}=\frac{\partial }{\partial x}\left(D_f{\partial }_t^{1-\beta }\frac{{\partial }^{\gamma }u}{\partial x^{\gamma }}\right),$$

(6)

where $\gamma$ and $\beta$ are coupled in such a way that $<x_0^2>\sim t_0^{\frac{2}{2+\theta }}$ is satisfied. It can also be shown that the correlation $t_0=x_0^{\frac{1+\gamma }{\beta }}$ holds. We can then conclude that $\theta +2=\frac{1+\gamma }{\beta }$. By applying the fractional integral operator to both sides of (6), we obtain the following useful form:

$$\frac{{\partial }^{\beta }u}{\partial t^{\beta }}=\frac{\partial }{\partial x}\left(D_f\frac{{\partial }^{\gamma }u}{\partial x^{\gamma }}\right).$$

(7)

The Equation (7) can be found in the literature on the diffusion phenomena in the chaotic migration of the particles and anomalous contaminant diffusion from a fracture into a porous rock matrix with an alteration zone bordering the fracture, see [13,14,15,16,17] and some of the references cited therein for further reading. In addition, in [7,18], the case where $\gamma =1$ is considered under the assumption the porous medium has a comb-like structure of fractal geometry. We refer the readers to [15] for more details about the Equation (7).

When a porous medium equation is considered, the pressure is taken to be a monotone function of the concentration u. In this case, by using the corresponding Darcy’s Law and the continuity equation, the following nonlinear equation is obtained:

$$u_t+\nabla ·\left(d\left(u\right)\nabla u\right)=0.$$

(8)

The Equation (8) is also known as the Richards’ equation in hydrology [19].

For the convenience of the reader, we present the main steps of the derivation of the governing equation by following [20]. If the fluid particles are trapped in some region for several periods of time $s_1,\cdots ,s_n$, the continuity equation becomes

$$u_t=-\sum _{i=1}^n{w}_i\nabla ·q(x,t-s_i),$$

(9)

where $w_i$ are some weights. If we consider the case with $w_i=w\left(s_i\right)\Delta s_i$ and $\Delta s_i=s_i-s_{i-1}$ for some weight density $w=w\left(s\right)$, and calculate the limit of $u_t$ as n→∞, the Equality (9) becomes

$$u_t=-{\int }_0^{t}w(t-s)\nabla ·q(x,s)ds.$$

(10)

The Equation (10) accounts for the fluid particles that can be trapped for any period of time. The function $w\left(s\right)$ represents the amount of flux of the particles that waits for the time s. Using the choice for w used in [20], we conclude the following:

$$u_t=-\frac{1}{\Gamma \left(\beta \right)}\frac{\partial }{\partial t}{\int }_0^t{(t-s)}^{\beta -1}\nabla ·q(x,s)ds={-}^R{\partial }_t^{1-\beta }\nabla ·q.$$

(11)

If we apply the Riemann-Liouville fractional operator $I_t^{1-\beta }$ to both sides of (11) and take into account the composition formula for the functions with vanishing initial conditions, we obtain

$$\frac{{\partial }^{\beta }}{\partial t^{\beta }}u=\nabla ·\left(d\left(u\right)\nabla u\right).$$

(12)

The Equation (12) is also called the generalized Richards equation [21] and can be found in many papers. In [21], the Equation (12) is solved numerically and the solution to the data on horizontal water transport is fitted. The numerical solution is also studied by some authors, see [22,23,24,25]. In [26], magnetic resonance imaging was employed to study water ingress in fine zeolite powders compacted by high pressure. The measured moisture profiles indicate sub-diffusive behavior with a spatio-temporal scaling variable $\eta =\frac{x}{t^{\gamma /2}}$. In [26], the Equation (12) is used to analyze the data, and an expression that yields the moisture dependence of the generalized diffusivity is derived and applied to their measured profiles. In [27], the authors use a time-fractional diffusion equation for the modeling of the probability density function of displacements. In [28], a method of approximating equations with the Erdelyi-Kober fractional operator which arise in mathematical descriptions of anomalous diffusion has been introduced. A theorem is proven on the exact form of the approximating series. The authors also provide an illustration obtained by the fractional porous-medium equation that is used to model moisture diffusion in building materials.

The authors in [29] look for self-similar solutions for the Equation (12). The resulting similarity equations are of nonlinear integro-differential type. They approximate these equations by an expansion of the integral operator and look for solutions in a power function form. Several applications of the Equation (12) are presented in [30]. In addition to what is already mentioned above, there has been a growing interest in inverse problems with fractional derivatives. These problems are physically and practically very important. Interested readers can be referred to [31,32,33,34,35,36,37,38,39,40,41,42,43] for a deeper understanding of the topic.

In [30], the authors proposed a generalization of Richards’ equation to estimate super-diffusion and sub-diffusion by means of the power-law ruler. Resulting solutions to the fractional Richards’ equation display anomalous non-Boltzmann scaling as a result of fractal medium of heterogeneous media. In [28], a fractional porous-medium equation to model moisture diffusion in building materials is considered. To achieve this goal, the author improved a method of approximating equations by means of Erdélyi-Kober fractional operator which is observed in equations describing anomalous diffusion. In [26], a time-fractional diffusion equation model of anomalous diffusion is adopted to analyze the data of moisture profiles. The authors consider a nonlinear diffusion equation with a time-fractional derivative which is a generalization of the porous medium equation and they find self-similar solutions to that problem in [29].

In this paper, we study an inverse coefficient problem for the nonlinear time-fractional diffusion Equation (12). The difference of our study from the studies in the references [31,32,33,34,35,36,37,38,39,40,41,42,43] is that the unknown of the inverse problem is non-linear, i.e., it depends on the solution u. This is a relatively new topic and there are only a few works, see [44,45,46]. In [44], the unknown coefficient depends on the gradient of the solution and belongs to a set of admissible coefficients. The authors prove that the direct problem has a unique solution and set up the continuous dependence of the solution of the corresponding direct problem on the unknown coefficient. Then, the existence of a quasi-solution of the inverse problem is shown in the appropriate class of admissible coefficients.

In [44], the numerical solutions of the direct and the inverse problems have been introduced and in [45] an application of the governing equation in the materials sciences has been mentioned. An inverse problem for the nonlinear time-fractional diffusion (12) is studied in [46]. In this paper, the authors prove that the direct problem has a unique solution. The existence of a quasi-solution is also proved. However, neither the uniqueness of the solution is proved nor the direct and the inverse problems are solved numerically. Therefore, our study can be considered as the continuation of the series of work in [46] and the works mentioned above on fractional inverse problems.

It is also worth mentioning that some inverse problems are studied for the case $\beta =1$ in (12). For example, in [47], the determination of the unknown coefficient $d\left(u\right)$ from over-specified data measured at the boundary has been studied. The inverse problem is reformulated as an auxiliary inverse problem and it is shown that this auxiliary problem has at least one solution in a specified admissible class. Finally, the auxiliary problem is approximated by an associated identification problem and some numerical results are presented. In [48], an operator approach is improved by an input-output mapping and it is shown that the mapping is isotonic. This result is used to derive a uniqueness result for the inverse problem.

This paper is organized as follows: In Section 2 we formulate the direct and the inverse problems. Existence and uniqueness for the direct and the inverse problems are discussed in Section 3. The numerical solutions of the direct and the inverse problems are studied in Section 4 and Section 5, respectively. The conclusion and the potential directions of improvement on the problem are presented in Section 6.

2. Formulation of the Direct and the Inverse Problems

In this section we formulate the direct and the inverse problems. First, we consider the following problem:

$$\left\{\begin{array}{cc}\frac{{\partial }^{\beta }u}{\partial t^{\beta }}=\nabla ·\left(d\left(u\right)\nabla u\right)+f(x,y,t),\hfill & (x,y,t)\in {\Gamma }_T,\\ u(x,y,t)=0,\hfill & (x,y,t)\in {\Gamma }_{3T}\cup {\Gamma }_{4T},\\ -d\left(u\right)u_y(x,y,t)=g_1(x,t),\hfill & (x,y,t)\in {\Gamma }_{1T},\\ -d\left(u\right)u_x(x,y,t)=g_2(y,t),\hfill & (x,y,t)\in {\Gamma }_{2T},\\ u(x,y,0)=0,\hfill & (x,y)\in \overline{\Omega },\end{array}\right.$$

(13)

where $\beta$ is the order of the Caputo fractional time derivative $,\Omega :=(0,1)\times (0,1)$, ${\Gamma }_T:=\Omega \times (0,T)$, ${\Gamma }_{iT}:={\Gamma }_i\times (0,T),$ for $i=1,2,3,4$ and $T>0$ is the end time. We assume that $\Omega$ is a bounded simply connected domain with a piece-wise smooth boundary $\partial \Omega ={\Gamma }_1\cup {\Gamma }_2\cup {\Gamma }_3\cup {\Gamma }_4$ with ${\Gamma }_i\cap {\Gamma }_j=\varnothing$, $i\ne j$. We define ${\Gamma }_i,i=1,2,3,4$ as follows:

$$\begin{array}{c}\hfill {\Gamma }_1:=(0,1)\times \left\{1\right\},{\Gamma }_2:=\left\{1\right\}\times (0,1),\\ \hfill {\Gamma }_3:=(0,1)\times \left\{0\right\},{\Gamma }_4:=\left\{0\right\}\times (0,1).\end{array}$$

The direct problem consists of determining the function $u(x,y,t)$ from the problem (13) for the given inputs $\beta$, $d\left(u\right)$, $f(x,y,t)$, $g_1(x,t)$ and $g_2(y,t)$. We solve the direct problem (13) in Section 4 by converting the equation in (13) into a fractional system using the method of lines. Next, we define the class of admissible coefficients and the weak solution of the problem (13). We note that, throughout the paper, $\parallel ·\parallel$ and $〈,〉$ denote the usual $L^2(\Omega )$ norm and inner product, respectively, ${\parallel ·\parallel }_X$ denotes the norm in a Hilbert space X. $C\left(I\right)$ denotes the set of continuous functions defined on I.

Definition 1.

Let $I:=[0,\infty )$. A set $\mathbb{D}$ satisfying the following conditions is called the class of admissible coefficients for the problem (13):

$$d\in C\left(I\right),c_0\le d\left(s\right)\le c_1,\forall s\in I,$$

(14)

$$\left(d\left(u_1\right)\nabla u_1-d\left(u_2\right)\nabla u_2\right)·\nabla (u_1-u_2)\ge c_2\left|\right|\nabla (u_1-u_2){\left|\right|}^2,\forall u_1,u_2\in H_0^1(\Omega ),$$

(15)

where $c_0,c_1,c_2$ are positive constants.

Definition 2.

A weak solution for the problem (13) is a function $u\in S^{\beta }\left({\Gamma }_T\right):=L^2\left(0,T;H_0^1(\Omega )\right)\cap W_2^{\beta }\left(0,T;L^2(\Omega )\right)$ such that the integral identity

$$\begin{array}{cc}\hfill \underset{\Omega }{\int }\frac{{\partial }^{\beta }u}{\partial t^{\beta }}vdxdy& +\underset{\Omega }{\int }d\left(u\right)\nabla u·\nabla vdxdy\hfill \\ \hfill & =\underset{\Omega }{\int }fvdxdy+{\int }_{{\Gamma }_3}{g}_{1}vdx+{\int }_{{\Gamma }_4}{g}_{2}vdy,\hfill \end{array}$$

(16)

holds for almost every $t\in [0,T]$ and for each $v\in S^{\beta }\left({\Gamma }_T\right)$, where

$$W_2^{\beta }(0,T):=\left\{u\in L^2[0,T]:\frac{{\partial }^{\beta }u}{\partial t^{\beta }}\in L^2[0,T]\mathrm{and}u\left(0\right)=0\right\}$$

is the fractional Sobolev space of order β. We note that $S^{\beta }\left({\Gamma }_T\right)$ is a Banach space with the norm:

$${\parallel u\parallel }_{S^{\beta }\left({\Gamma }_T\right)}={\left({\parallel u\parallel }_{W_2^{\beta }\left(0,T;L^2(\Omega )\right)}^2+{\parallel u\parallel }_{L^2\left(0,T;H_0^1(\Omega )\right)}^2\right)}^{\frac{1}{2}}.$$

We consider $u(x,y,t)$ as a mapping from $t\in (0,T)$ to $L^2(\Omega )$ and write as $u\left(t\right)=u(·,·,t)$.

A weak solution of the problem (13) can also be defined as a solution of the following abstract operator equation

$$Lu+Au=F,$$

(17)

where $Lu:=〈\widehat{L}u,v〉,\widehat{L}u:=\frac{{\partial }^{\beta }u}{\partial t^{\beta }}$, $\widehat{L}:D\left(\widehat{L}\right)\subset V\to V^{*}$ with the domain $D\left(\widehat{L}\right)=\left\{u\in V:\frac{{\partial }^{\beta }u}{\partial t^{\beta }}\in V^{*}\right\}$, $V:=L^2\left(0,T;H_0^1(\Omega )\right)$, the nonlinear operator $A:V\to V^{*}$ is defined by

$$〈Au,v〉:=\underset{\Omega }{\int }\frac{{\partial }^{\beta }u}{\partial t^{\beta }}vdxdy+\underset{\Omega }{\int }d\left(u\right)\nabla u·\nabla vdxdy,$$

(18)

and F on V is defined by

$$〈F,v〉:=\underset{\Omega }{\int }fvdxdy+{\int }_{{\Gamma }_3}{g}_{1}vdx+{\int }_{{\Gamma }_4}{g}_{2}vdy.$$

(19)

The inverse problem here consists of determining the pair of functions $\left\{u(x,y,t),d\left(u\right)\right\}$ from the problem (13) by means of the additional data $u(x,y,t)={\widehat{g}}_1(x,t),(x,y,t)\in {\Gamma }_{1T}$ and $u(x,y,t)={\widehat{g}}_2(y,t),(x,y,t)\in {\Gamma }_{2T}$. We solve the inverse problem by minimizing the cost functional defined by (22). For the consistency of the additional data with the data of (13) on ${\Gamma }_{1T}$, it is assumed that $-d\left({\widehat{g}}_1(x,t)\right){\left(\widehat{g_1}\right)}_y(x,y,t)=g_1(x,t)$ on ${\Gamma }_{1T}$ and $-d\left({\widehat{g}}_2(x,t)\right){\left(\widehat{g_2}\right)}_y(x,y,t)=g_1(x,t)$ on ${\Gamma }_{2T}$.

We denote the solution of the direct problem (13) for a given function $d\in \mathbb{D}$ by $u(x,y,t;d)$. If the function $u(x,y,t;d)$ also satisfies the additional data above, it is called a strict solution of the inverse problem. Now, we reformulate the inverse problem. For this purpose, we define the input-output map $\Phi :L^2\left({\Gamma }_T\right)\to L^2\left({\Gamma }_{3T}\right)\times L^2\left({\Gamma }_{4T}\right)$ as

$${\Phi \left(d\right):=u(x,y,t;d)|}_{{\Gamma }_{3T}\times {\Gamma }_{4T}}:=\left({u(x,y,t;d)|}_{{\Gamma }_{3T}}{,u(x,y,t;d)|}_{{\Gamma }_{4T}}\right).$$

(20)

Then, the inverse problem is defined as a solution of the following operator equation:

$$\Phi \left(d\right)=g,g=({\widehat{g}}_1(x,t),{\widehat{g}}_2(y,t))\in L^2\left({\Gamma }_{3T}\right)\times L^2\left({\Gamma }_{4T}\right).$$

(21)

However, due to the measurement errors in practice, the exact equality in (21) is usually not achieved. Hence, one needs to introduce the auxiliary (error) functional:

$$I\left(d\right):={\int }_0^T\underset{{\Gamma }_3}{\int }|u(x,y,t;d)-{\widehat{g}}_1(x,t){|}^{2}dxdt+{\int }_0^T\underset{{\Gamma }_4}{\int }|u(x,y,t;d)-{\widehat{g}}_2(y,t){|}^{2}dydt,$$

(22)

and consider the following minimization problem:

$$I\left(\widehat{d}\right)=\underset{d\in \mathbb{D}}{min}I\left(d\right).$$

(23)

A solution of the minimization problem (23) is called a quasi-solution (or approximate solution) of the inverse problem. Evidently, if $I(\widehat{d})=0$, then the quasi-solution $\widehat{d}$ is also a strict solution of the inverse problem.

3. Analysis of the Direct and the Inverse Problems

In this section, we analyse both the direct and the inverse problems. The theoretical aspect of the direct problem (13) is studied in [46]. In this study, the authors prove that the direct problem (13) is well-posed in the sense of Hadamard. For the sake of the reader, we provide some relevant results from [46]. The following theorems state that the direct problem (13) has a unique weak solution that depends continuously on the coefficient $d\left(u\right)$. We refer the readers to [46] for detailed proofs.

Theorem 1.

Let $d\in \mathbb{D}$. Then, the direct problem (13) has a unique weak solution $u\in S^{\beta }\left({\Gamma }_T\right)$. Moreover, for almost every $t\in [0,T]$ there exist some constants $c,C>0$ such that

$$\frac{{\partial }^{\beta }{\parallel u\parallel }^2}{\partial t^{\beta }}+c{\parallel u\parallel }_{H_0^1(\Omega )}^2\le C\left[{\parallel f\parallel }^2+\parallel g_1{\parallel }_{L^2\left({\Gamma }_1\right)}^2+{\parallel g_2\parallel }_{L^2\left({\Gamma }_2\right)}^2\right].$$

Theorem 2.

Suppose that a sequence of coefficients $\left\{d_m\right\}\subset \mathbb{D}$ converges pointwise in $[0,\infty )$ to a function $d\in \mathbb{D}$. Then, the sequence of solutions $u_m:=u(x,y,t;d_m)$ converges to the solution $u:=u(x,y,t;d)\in S^{\beta }\left({\Gamma }_T\right)$, where $u:=u(x,y,t;d)$ denotes the solution of the direct problem (13) for a given coefficient $d\in \mathbb{D}$.

Next, we prove a theorem for the existence of the solution to the inverse problem. There are two methods in the literature to prove the existence of the solution of inverse problems. The first method is called the monotonicity method, which is based on the continuity and the monotonicity of the input-output mapping. The second method is called quasi-solution method, which is based on minimizing an error functional between the output data and the additional data. We adopt the quasi-solution approach to the inverse problem under consideration. For this purpose, we show that the cost functional defined by (22) is continuous and we construct a compact subset of the class of the admissible coefficients.

Theorem 3.

Assume that a sequence of coefficients $\left\{d_m\right\}\subset \mathbb{D}$ converges pointwise in $[0,\infty )$ to a function $d\in \mathbb{D}$. Then, $|I\left(d_m\right)-I\left(d\right)|\to 0$ as $n\to \infty$.

Proof. 

Let $\left\{d_m\right\}\subset \mathbb{D}$ be a sequence of coefficients that converges pointwise in $[0,\infty )$ to a function $d\in \mathbb{D}$, $u_m:=u(x,y,t;d_m)$ and $u:=u(x,y,t;d)$. Then, we have:

$$\begin{array}{cc}\hfill |I\left(d_m\right)-I\left(d\right)|& =|{\int }_0^T\underset{{\Gamma }_3}{\int }|u_m-{\widehat{g}}_1(x,t){|}^{2}dxdt-{\int }_0^T\underset{{\Gamma }_3}{\int }|u-{\widehat{g}}_1(x,t){|}^{2}dxdt\hfill \\ \hfill & +{\int }_0^T\underset{{\Gamma }_4}{\int }|u_m-{\widehat{g}}_2(y,t){|}^{2}dydt-{\int }_0^T\underset{{\Gamma }_4}{\int }|u-{\widehat{g}}_2(y,t){|}^{2}dydt|.\hfill \end{array}$$

(24)

For the first two terms in (24), we have:

$$\begin{array}{cc}\hfill & {\int }_0^T\underset{{\Gamma }_3}{\int }|u_m-{\widehat{g}}_1(x,t){|}^{2}dxdt-{\int }_0^T\underset{{\Gamma }_3}{\int }|u-{\widehat{g}}_1(x,t){|}^{2}dxdt\hfill \\ \hfill & =|\parallel u_m-{\widehat{g}}_1(x,t){\parallel }_{L^2\left({\Gamma }_{3T}\right)}^2-\parallel u-{\widehat{g}}_1(x,t){\parallel }_{L^2\left({\Gamma }_{3T}\right)}^2|\hfill \\ \hfill & \le \parallel u_m-u{\parallel }_{L^2\left({\Gamma }_{3T}\right)}\times (\parallel u_m-{\widehat{g}}_1(x,t){\parallel }_{L^2\left({\Gamma }_{3T}\right)}\hfill \\ \hfill & +\parallel u-{\widehat{g}}_1(x,t){\parallel }_{L^2\left({\Gamma }_{3T}\right)})\hfill \\ \hfill & \le \parallel u_m-u{\parallel }_{L^2\left({\Gamma }_T\right)}\times (\parallel u-{\widehat{g}}_1(x,t){\parallel }_{L^2\left({\Gamma }_{3T}\right)}\hfill \\ \hfill & +\parallel u-{\widehat{g}}_1(x,t){\parallel }_{L^2\left({\Gamma }_{3T}\right)}).\hfill \end{array}$$

(25)

By using the fact that $\parallel u{\parallel }_{L^2\left(0,T;H_0^1(\Omega )\right)}\le K\parallel u{\parallel }_{S^{\beta }\left({\Gamma }_T\right)}$ for some $K>0$ (see [49] for details), we conclude that the following holds for the first term on the right hand side of (25):

$$\parallel u_m-u{\parallel }_{L^2\left({\Gamma }_T\right)}\le \parallel u_m-u{\parallel }_{L^2\left(0,T;H_0^1(\Omega )\right)}\le K\parallel u_m-u{\parallel }_{S^{\beta }\left({\Gamma }_T\right)}.$$

(26)

From (26) and Theorem 2, we conclude that the first two terms in (24) goes to zero as $n\to \infty$. Similarly, we can prove that the last two terms in (24) tend to zero as $n\to \infty$. This completes the proof. □

The conditions (14) and (15) arise in the solvability of the direct problem (13) and can be found in some papers, for example see the condition H3 in [50]. In virtue of Theorem 3, it is natural to construct a compact set of admissible coefficients in $C\left(I\right)$. For this reason, in addition to the assumptions (14) and (15), we assume that there is a subset ${\mathbb{D}}_c$ of $\mathbb{D}$ which is equicontinuous, i.e., ${\mathbb{D}}_c\subset \mathbb{D}$ and for every $ϵ>0$ there exists a $\delta >0$ such that if $d\in {\mathbb{D}}_c,s_1,s_2\in I$ and $|s_1-s_2|<\delta$, then $|d\left(s_1\right)-d\left(s_2\right)|<ϵ$. By following Theorem 3 in [51] it can be proved that ${\mathbb{D}}_c$ is compact. Then, we obtain the following existence theorem.

Theorem 4.

The inverse problem has at least one quasi-solution in the set of admissible coefficients ${\mathbb{D}}_c$.

Proof. 

The Weierstrass theorem and Theorem 3 guarantees the existence of the solution to the inverse problem is proved. □

The following theorem shows that the input-output operator defined by (20) is a compact operator.

Theorem 5.

Let the conditions (14) and (15) hold. Then, the input-output operator defined by (20) is a compact operator [46].

Since nonlinear equations with compact operators are ill-posed [52], the inverse problem under consideration is also ill-posed.

The following example shows that there exists a sequence $\left\{d_n\right\}$ such that $u(x,y,t;d_n),(x,y,t)$ $\in {\Gamma }_{1T}$ and $u(x,y,t;d_n),(x,y,t)\in {\Gamma }_{2T}$ converges to zero as $n\to \infty$, but $d_n\to \infty$ as $n\to \infty$.

Example 1.

For $d_n\left(u\right)=n^2$, $\beta =1/2$ and $f_n(x,t)=\frac{8x^3{y}^2{t}^{3/2}}{3\sqrt{\pi }{n}^5}-\frac{2x{t}^2}{n^3}(3y^2+x^2)$, the inverse problem (13) becomes

$$\left\{\begin{array}{c}\frac{{\partial }^{\beta }u}{\partial t^{\beta }}={\left(n^2{u}_x\right)}_x+{\left(n^2{u}_y\right)}_y+\frac{8x^3{y}^2{t}^{3/2}}{3\sqrt{\pi }{n}^5}-\frac{2x{t}^2}{n^3}(3y^2+x^2),(x,y,t)\in {\Gamma }_T,\hfill \\ u(x,y,t)=0,(x,y,t)\in {\Gamma }_{3T}\cup \in {\Gamma }_{4T},\hfill \\ -n^2{u}_y(x,y,t)=-\frac{2}{n^3}{x}^3{t}^2,(x,y,t)\in {\Gamma }_{1T},\hfill \\ -n^2{u}_x(x,y,t)=-\frac{3}{n^3}{y}^2{t}^2,(x,y,t)\in {\Gamma }_{2T},\hfill \\ u(x,y,0)=0,(x,y)\in \overline{\Omega }.\hfill \end{array}\right.$$

(27)

It can easily be verified that the function $u_n(x,y,t;d_n)=\frac{x^3{y}^2}{n^5}{t}^2$ is the solution of the corresponding direct problem. Obviously, $u_n(x,1,t;d_n)$ and $u_n(1,y,t;d_n)$ converge to zero as $n\to \infty$, but $d_n\left(u\right)\to \infty$ as $n\to \infty$.

4. Numerical Solution of the Direct Problem

In this section we introduce the methodology used for solving the direct problem numerically. The main idea is to convert the fractional partial differential equation into a system of fractional ordinary differential equations using the method of lines and vectorization and then solve the resulting system of fractional ordinary differential equations. In addition to the classical method of lines, we adopt the operator approach to approximate derivatives, which reduces computations and memory demand of the algorithm. We first illustrate the methodology over a one-dimensional heat equation. For this purpose, we consider the following one-dimensional problem:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{\partial u}{\partial t}=k\frac{{\partial }^{2}u}{\partial x^2},x\in (0,1),t\in (0,1),\hfill \\ u(x,0)=h\left(x\right),x\in (0,1),\hfill \\ u(0,t)=f\left(t\right),u(1,t)=g\left(t\right),t\in (0,1).\hfill \end{array}\right.\end{array}$$

For a given positive integer M, let $x_i=i\Delta x$ for $i=0,1,2,\cdots ,M$ with $\Delta x=1/M$. Additionally, let $u_i\left(t\right)$ denote the approximation of the solution at the node $(x_i,t)$ for $i=0,1,2,\cdots ,M$, where $u_0\left(t\right)=f\left(t\right)$ and $u_M\left(t\right)=g\left(t\right)$. We approximate the given equation by the following system of ordinary differential equations for $u_i$’s:

$$\left\{\begin{array}{c}\frac{d{u}_i}{dt}=k\frac{u_{i+1}-2u_i+u_{i-1}}{{\Delta x}^2},t\in (0,1),i=1,\cdots ,M-1,\hfill \\ u_i\left(0\right)=h\left(x_i\right),i=1,\cdots ,M-1.\hfill \end{array}\right.$$

(28)

Given the vector of approximate solutions at each node without boundaries, $\left[u_i\right]=[u_1,u_2,\cdots ,u_{M-1}]$, we define the left and right shift operators as

$$\begin{array}{cc}\hfill & LS\left(\left[u_i\right]\right)=\left[u_{i+1}\right]\mathrm{and}RS\left(\left[u_i\right]\right)=\left[u_{i-1}\right],i=0,1,2,\cdots ,M-1.\hfill \end{array}$$

(29)

Then, the system of ordinary differential equations given in (28) can be written as follows:

$$\begin{array}{c}\hfill \left[\frac{d{u}_i}{dt}\right]=\frac{LS\left(\left[u_i\right]\right)-2\left[u_i\right]+RS\left(\left[u_i\right]\right)}{{\Delta x}^2},\end{array}$$

(30)

with $\left[\frac{d{u}_i}{dt}\right]=\left[\frac{d{u}_1}{dt},\frac{d{u}_2}{dt},\cdots ,\frac{d{u}_{M-1}}{dt}\right]$ being the vector of time derivatives at the discretized nodes. To illustrate how we generalize the shift-operator approach to higher dimensional partial differential equations, we consider the following two-dimensional problem:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}\frac{\partial u}{\partial t}=\frac{{\partial }^{2}u}{\partial x^2}+\frac{{\partial }^{2}u}{\partial y^2}+f(x,y,t),(x,y,t)\in {\Gamma }_T,\hfill \\ u(x,y,0)=h(x,y),(x,y)\in \Omega ,\hfill \\ u(x,y,t)=g(x,y,t),(x,y)\in \partial \Omega ,t\in (0,1).\hfill \end{array}\right.\end{array}$$

(31)

Let $x_i=i\Delta x,i=0,\cdots ,M$ and $y_j=j\Delta y,j=0,1,2,\cdots ,N$ with $\Delta x=1/M$ and $\Delta y=1/N$. Additionally, let $u_{ij}\left(t\right)$ denote the solution at point $(x_i,y_j)$ at a time $t\in (0,1)$, where $u_{0j}=g(0,y_j,t)$, $u_{Mj}=g(1,y_j,t)$, $u_{i0}=g(x_i,0,t)$ and $u_{iN}=g(x_i,1,t)$. Then, the centered difference approximation of the time derivative at point $(x_i,y_j)$ will be

$$\begin{array}{cc}\hfill & \frac{d{u}_{ij}}{dt}=\frac{u_{i+1j}-2u_{ij}+u_{i-1j}}{\Delta x^2}+\frac{u_{ij+1}-2u_{ij}+u_{ij-1}}{\Delta y^2}+f(x_i,y_j,t),\hfill \end{array}$$

with $i=1,\cdots ,M-1\mathrm{and}j=1,\cdots ,N-1.$ This approximation can be vectorized by defining the solution matrix at the interior points $\left[u_{ij}\right]=u(x_i,y_j)$ with $i=1,\cdots ,M-1$ and $j=1,\cdots ,N-1$. We define the left and right shift operators on the matrix $\left[u_{ij}\right]$ of solution approximations at the interior points as follows:

$$\begin{array}{cc}\hfill & LS\left(\left[u_{ij}\right]\right)=\left[u_{i+1j}\right]\mathrm{and}RS\left(\left[u_{ij}\right]\right)=\left[u_{i-1j}\right],\hfill \end{array}$$

with $i=1,\cdots ,M-1\mathrm{and}j=1,\cdots ,N-1.$ Then the matrices of the centered difference approximations to the first order derivatives can be expressed as follows:

$$\begin{array}{c}\hfill \left[u_{xij}\right]=\frac{LS\left(\left[u_{ij}\right]\right)-RS\left(\left[u_{ij}\right)\right]}{2\Delta x},\end{array}$$

(32)

$$\begin{array}{c}\hfill \left[u_{yij}\right]=\frac{LS\left({\left[u_{ij}\right]}^{\prime }\right)-RS\left({\left[u_{ij}\right]}^{\prime }\right)}{2\Delta y},\end{array}$$

(33)

where ${\left[a_{ij}\right]}^{\prime }$ denotes the transpose of the matrix $\left[a_{ij}\right]$. The matrices of the centered difference approximations to the second order derivatives, $\left[u_{xxij}\right],\left[u_{yyij}\right]$, and $\left[u_{xyij}\right]$ can be obtained by applying $LS$ and $RS$ operators to the matrices of the first order derivative approximations given in (32) and (33). Then, (31) can be expressed in matrix form as follows:

$$\begin{array}{c}\hfill \begin{array}{c}\left[\frac{d{u}_{ij}}{dt}\right]=\left[u_{xxij}\right]+\left[u_{yyij}\right]+\left[f(x_i,y_j,t)\right],\hfill \end{array}\end{array}$$

(34)

with initial condition matrix $\left[u_{ij}\left(0\right)\right]=\left[h(x_i,y_j)\right]$ of size $M-1\times N-1$.

Next, we consider the following time-fractional problem:

$$\left\{\begin{array}{c}\frac{{\partial }^{\beta }u}{\partial t^{\beta }}=f(x,y,t,u,u_x,u_y,u_{xx},u_{yy},u_{xy}),(x,y,t)\in {\Gamma }_T,\hfill \\ u(x,y,0)=h(x,y)(x,y)\in \Omega ,\hfill \\ u(x,y,t)=g(x,y,t),(x,y)\in \partial \Omega ,t\in (0,1).\hfill \end{array}\right.$$

(35)

The vectorized method of lines approach described in examples above results in the following difference approximation for the equation in (35):

$$\left[\frac{{\partial }^{\beta }{u}_{ij}}{\partial t^{\beta }}\right]=f(\left[x_{ij}\right],\left[y_{ij}\right],t,\left[u_{ij}\right],\left[u_{xij}\right],\left[u_{yij}\right],\left[u_{xxij}\right],\left[u_{yyij}\right],\left[u_{xyij}\right])],$$

(36)

with $\left[u_{ij}\left(0\right)\right]=\left[h(x_i,y_j)\right]$$\left[x_{ij}\right]={\left[x_i\right]}_{1\times M-1}\otimes {\left[1\right]}_{N-1\times 1}$ and $\left[y_{ij}\right]={\left[y_j\right]}_{N-1\times 1}\otimes {\left[1\right]}_{M-1\times 1}$, where ⊗ denotes the Kronecker matrix product. We solve this system of nonlinear fractional ordinary differential equations using a Matlab implementation of the Adam-Bashfort-Moulton(ABM) type predictor-corrector method given in [53]. Detailed convergence and stability analysis are considered both in [54,55]. In [54], it is concluded that if the solution under consideration is sufficiently smooth, the method has uniform convergence of order $h^2$ for $\beta >1$, and of order $h^{1+\beta }$ for $\beta <1$, respectively. It is further shown by numerical examples that, these bounds are strict and cannot be improved.

The ABM is a Predict-Evaluate-Correct-Evaluate (PECE)-type method. That is, for the approximation at the time nodes $t_k$, and corresponding approximations, $y_j\cong y\left(t_k\right)$ at each $k_{th}$ step, there are two approximations computed for the next node, namely, predictor, $yp\left(t_{k+1}\right)$, and using the predictor, the corrector approximation $yc\left(t_{k+1}\right)$ is obtained and used in the calculation. The approximation error is obtained by finding the difference of predictor and corrector approximations. There are two main advantages of using a PECE-type method compared to the classical equivalent-order explicit methods. The first benefit of using an algorithm of type PECE is the high accuracy and stability, see [55,56,57,58]. For fractional ordinary differential equations, it has been shown that the stability and accuracy remains high compared to equivalent-order numerical methods, see [53,54,55]. The second advantage of using a PECE-type numerical algorithm is that these methods use variable time steps that reduce the computational cost of the approximation. The method can control the time steps by using the difference between the corrector and the predictor approximations. When the difference is smaller than the desired level of accuracy with the current time step, this can be used as an indication that the solver is in a non-stiff area, and consequently time steps are increased in an adaptive manner.

The idea of combining the method of lines approach to reduce the given fractional partial differential equation to a system of fractional ordinary differential equations and using shift operators in the evaluation of the right-hand side of the PDE can prove to be useful in terms of memory and computation compared to similar operator approaches, such as that of Podlubny’s intuitive matrix operator approach, see [59,60], depending on the problem under consideration. In terms of memory, the method of lines approach uses matrices of size $M\times N$. Whereas in Podlubny’s Matrix Operator approach, the matrices under consideration are of size $M\times N\times K$ where K is the number of mesh points in t direction used in the calculation. So, the method of lines approach is able to improve the memory demand of the calculations. However, this improvement in memory demand comes at the cost of calculation of the right-hand side function at each time step. Similarly, Podlubny’s matrix operator approach may also face the computational challenges depending on the complexity of the fractional partial differential equation. This is because the method requires solving nonlinear algebraic matrix equation of very high dimensions. Hence, the very nature of the question under consideration is the sole factor in choosing the numerical method to apply.

The first series of the numerical simulations is related to numerical solution of the direct problem. For this purpose, the function $u(x,y,t)=t{x}^2{y}^2$ is taken to be the analytic solution of the equation $\frac{{\partial }^{\beta }u}{\partial t^{\beta }}={\left(d\left(u\right)u_x\right)}_x+{\left(d\left(u\right)u_y\right)}_y+f(x,y,t),$ with the function $d\left(u\right)=1+u$ and appropriately chosen source function $f(x,y,t)$. The boundary conditions are found from the trace of the function $u(x,y,t)=t{x}^2{y}^2$ on $\Omega$. First, we check the difference between the numerical solution $u_{\mathrm{num}}$ and the exact solution $u_{exact}=u(x,y,t)=t{x}^2{y}^2$. The absolute error between $u_{\mathrm{num}}$ and $u_{\mathrm{exact}}$ is defined by $\parallel u_{\mathrm{num}}-u_{\mathrm{exact}}{\parallel }_{\infty }$, where ${\parallel ·\parallel }_{\infty }$ denotes the supremum norm, which is taken over all $x,y,t$ where $x\in [0,1],y\in [0,1]$ and $t\in [0,0.3]$. The software is simulated for different values of the system parameters to evaluate their effect on the absolute error and the simulation time. The time step is taken $\Delta t={10}^{-5}$. We note that a decay in the time step may lead to and increase in the simulation time and has almost no effect on the computations. However, there is a subtle relation between the number of time nodes and x-y nodes. Increasing $\Delta t={10}^{-5}$ when $t\in [0,0.3]$ affects the stability of the simulation when higher number of x-y nodes is used. We simulated the solution for $M=N=5$, $M=N=10$ and $M=N=20$ for $\beta =0.6$ and $\beta =0.7$. We show the results with computation times in

Table 1. Absolute error between $u_{\mathrm{num}}$ and $u_{\mathrm{exact}}$ for $\beta =0.7$.

M (=N)	Absolute Error	CPU Time (seconds)
20	$0.2347$	17.22
10	$0.0847$	8.49
5	$0.021$	6.14

and

Table 2. Absolute error between $u_{\mathrm{num}}$ and $u_{\mathrm{exact}}$ for $\beta =0.6$.

M (=N)	Absolute Error	CPU Time (seconds)
20	$0.3135$	23.139
10	$0.0949$	8.397
5	$0.0232$	6.12

. We observe that $M=N=5$ serves the best time and the least absolute error. In addition, it appears that as $\beta$ increases towards 1, the absolute error decreases. Hence, the time step and number of discretization nodes for the remaining simulations are chosen to be seemingly the optimal values of ${10}^{-5}$ and $M=N=5$.

Next, we consider the ill-posedness of the inverse problem. A problem is called to be well-posed in the Hadamard sense if the solution exists, unique and depends continuously on the input data. Failure to comply any of the mentioned properties makes the given problem ill-posed. To show that the inverse problem under consideration is ill-posed, we simulate the direct problem for $a_1\left(u\right)=1+u$ and $a_2\left(u\right)=3+u$. The functions $u_1(x,y,1)$ and $u_2(x,y,1)$ are found through the numerical solution. Figure 1 shows that the graph of the difference $\mu u_1(x,y,1)-u_2(x,y,1)$ where we see the solutions on both ${\Gamma }_{3T}$ and ${\Gamma }_{4T}$ are very close to each other with an absolute difference less than $0.02$.

5. Numerical Solution of the Inverse Problem

We assume that for a fixed $\beta$, the following problem has the solution $\tilde{u}(x,y,t)$ for the specific functions $\tilde{d},f(x,y,t),g_1(x,t)$ and $g_2(y,t)$. In the inverse problem, we aim to try to find the function $\tilde{d}$ when only $\tilde{u}(x,1,t)$ and $\tilde{u}(1,y,t)$ are known about the solution $\tilde{u}$ while the boundary conditions and $f(x,y,t)$ are known and fixed. The boundaries are the same as in the Equation (13).

$$\left\{\begin{array}{cc}\frac{{\partial }^{\beta }u}{\partial t^{\beta }}=\nabla ·\left(d\left(u\right)\nabla u\right)+f(x,y,t),\hfill & (x,y,t)\in {\Gamma }_T,\\ u(x,y,t)=0,\hfill & (x,y,t)\in {\Gamma }_{3T}\cup \in {\Gamma }_{4T},\\ -d\left(u\right)u_y(x,y,t)=g_1(x,t),\hfill & (x,y,t)\in {\Gamma }_{1T},\\ -d\left(u\right)u_x(x,y,t)=g_2(y,t),\hfill & (x,y,t)\in {\Gamma }_{2T},\\ u(x,y,0)=0,\hfill & (x,y)\in \overline{\Omega }.\end{array}\right.$$

(37)

In our experiments, first we fix our solution $\tilde{u}(x,y,t)$ as $t{x}^2{y}^2$ and the specific function $d\left(u\right)$ beforehand. By using them we find the exact values of $g_1(x,t)$, $g_2(x,t)$, whereas $f(x,y,t)$ is found numerically. Now $g_1(x,t)$, $g_2(x,t)$ and $f(x,y,t)$ at hand, we take $\tilde{u}(1,y,t)=t{y}^2$ and $\tilde{u}(x,1,t)=t{x}^2$ as additional conditions for the inverse problem and we try to find $d\left(u\right)$. To do this, we set up the following error functional for each function d as in (23):

$$I\left(d\right):={\int }_0^T{\int }_0^1|u(x,y,t;d)-\tilde{u}(x,1,t){|}^{2}dxdt+{\int }_0^T{\int }_0^1|u(x,y,t;d)-\tilde{u}(1,y,t){|}^{2}dydt$$

where $u(x,y,t;d)$ denotes the solution of the direct problem for the function d. We can easily deduce that $I\left(\tilde{d}\right)=0$. Thus, we expect to find $\tilde{d}\left(u\right)$ as the minimizer of the error functional.

We carry out the experiment for the noise-free data and the noisy data. For the noisy data, we add the noise $0.1{\varphi }_1(x,t)$ and $0.1{\varphi }_2(x,t)$ to $g_1(x,t)$ and $g_2(x,t)$, respectively. Here, ${\varphi }_1(x,t)$ and ${\varphi }_2(x,t)$ assume the values in the form of $0.1\times m(x,t)$ where m takes random values from $[-1,1]$.

We will use polynomials to approximate the minimum $\tilde{d}$ in the light of Theorem 2. For $n>0$, we assume that $d\left(u\right)=d_n{u}^n+\cdots +d_0$. Thus, $I\left(d\right)$ is a function of the vectors $(d,d_0)\in {\mathbb{R}}^{n+1}$ where $d=(d_n,\cdots d_1).$

We apply the BGFS method, a Quasi-Newton method, to minimize the function. The method requires the evaluation of the gradient of the function $I\left(d\right)$. The BGFS method approximates the Hessian of the error function with a cubic line search procedure in each step. For further reading, we refer [24]. We choose the stopping criteria for the algorithm as $\left|\right|\nabla I\left(d\right)\left|\right|<{10}^{-6}$. The integration in $I\left(d\right)$ is approximated by the trapezoid rule.

In the experiments, in accordance with the observations given in Section 4 we assume $\beta =0.7$ for the order of fractional time derivative, $M=N=5$ are used to make a mesh-grid for x and y on $[0,1]\times [0,1]$ and time step is taken as $▵t={10}^{-5}$ with upper bound $T=0.3$.

We consider the variable $d\left(u\right)$ as a second degree polynomial. Hence, $I\left(d\right)$ is treated as a function of three variables. The initial points are chosen to be close to the coefficients of the second degree Taylor polynomial of $\tilde{d}\left(u\right)$, a random number between 0 and 1 is added to or subtracted from each component. In the

Table 3. Results for noise-free data for $\tilde{d}\left(u\right)=2sin\left(u\right)+cos\left(u\right)$.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error
−0.27408	2.0855	1.0987	−0.30691	1.9416	1.0172	$6.5695\times {10}^{-8}$	0.0053729
−0.32929	2.2625	0.73813	−0.41756	1.9666	1.0105	$2.3175\times {10}^{-8}$	0.0032792
−1.4234	2.7303	0.27877	−0.56634	2.0006	1.0014	$3.0948\times {10}^{-10}$	0.00055506
−0.93021	1.5114	1.1068	−0.79997	2.0539	0.98705	$4.4387\times {10}^{-8}$	0.0039033
−0.68482	1.4215	0.34624	−0.60656	2.0101	0.99864	$9.2876\times {10}^{-10}$	0.00053536
0.40488	2.2373	1.4942	−0.58484	2.0048	1.0003	$2.2985\times {10}^{-10}$	0.00029774
0.47975	2.4588	0.22095	−0.58776	2.0056	0.99998	$2.7193\times {10}^{-10}$	0.00034198
−0.06113	1.0369	1.715	−0.58706	2.0054	1.0001	$2.6604\times {10}^{-10}$	0.00034775
−0.38888	1.4532	0.096279	−0.35227	1.9509	1.0152	$4.5877\times {10}^{-8}$	0.0044349
−0.75806	2.5211	1.8909	−0.79146	2.0521	0.9875	$4.1074\times {10}^{-8}$	0.0037083
−0.76221	1.3759	0.80219	−0.63609	2.0172	0.99653	$3.2795\times {10}^{-9}$	0.0008257
−0.27825	1.6326	1.5	−0.19536	1.9171	1.0234	$1.3041\times {10}^{-7}$	0.0076122
−0.79668	2.0377	1.9047	−0.74601	2.0417	0.99027	$2.5509\times {10}^{-8}$	0.0028553
−0.92417	1.0867	1.6177	−0.68659	2.0283	0.99378	$1.0782\times {10}^{-8}$	0.0017539

,

Table 4. Results for noisy data for $\tilde{d}\left(u\right)=2sin\left(u\right)+cos\left(u\right)$.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error
−0.27408	2.0855	1.0987	−0.3141	1.9232	1.0344	0.0020018	0.0052846
−0.32929	2.2625	0.73813	−0.42473	1.9484	1.0275	0.0020017	0.0032272
−1.4234	2.7303	0.27877	−0.67733	2.0063	1.0119	0.0020016	0.0021874
−0.93021	1.5114	1.1068	−0.80617	2.0373	1.0027	0.0020016	0.0037943
−0.68482	1.4215	0.34624	−0.61309	1.9932	1.0148	0.0020016	0.0021119
0.40488	2.2373	1.4942	−0.69746	2.0106	1.0109	0.0020016	0.0022959
0.47975	2.4588	0.22095	−0.70218	2.012	1.0104	0.0020016	0.0023236
−0.06113	1.0369	1.715	−0.70344	2.0122	1.0104	0.0020016	0.0023527
−0.38888	1.4532	0.096279	−0.35922	1.9326	1.0323	0.0020017	0.0043445
−0.75806	2.5211	1.8909	−0.79917	2.0341	1.0044	0.0020016	0.0037813
−0.76221	1.3759	0.80219	−0.64273	1.9991	1.0135	0.0020016	0.002096
−0.27825	1.6326	1.5	−0.20201	1.899	1.0403	0.0020018	0.0075475
−0.79668	2.0377	1.9047	−0.75291	2.0236	1.0072	0.0020016	0.0029034
−0.92417	1.0867	1.6177	−0.69291	2.0101	1.0107	0.0020016	0.0022322

,

Table 5. Results for noisy data for $\tilde{d}\left(u\right)=2sin\left(u\right)+cos\left(u\right)$ with regularization parameter.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error	$\mathit{\lambda }$
−0.32929	2.2625	0.73813	−0.41297	1.936	1.0355	0.0020067	0.0037234	${10}^{-6}$
0.40488	2.2373	1.4942	0.229	1.7102	1.1339	0.0020458	0.021791	${10}^{-5}$
−0.68482	1.4215	0.34624	0.22395	1.4118	1.3433	0.0024042	0.083568	${10}^{-4}$
−0.92417	1.0867	1.6177	0.16702	1.1677	1.3914	0.0055819	0.31648	${10}^{-3}$

,

Table 6. Results for noise-free data for $\tilde{d}\left(u\right)=u^2+u+1$.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error
1.8143	1.5688	1.1067	0.98685	1.0027	0.99947	$1.8662\times {10}^{-10}$	0.00021487
1.2435	1.4694	0.038102	1.0944	0.981	1.0043	$8.5119\times {10}^{-9}$	0.0014074
1.9293	0.9881	1.0046	0.98892	1.0026	0.9992	$1.3021\times {10}^{-10}$	0.00018254
0.8034	1.1622	1.8173	0.81833	1.0366	0.99164	$3.1636\times {10}^{-8}$	0.0027013
1.3517	1.1656	1.2599	1.3088	0.9365	1.0151	$9.1374\times {10}^{-8}$	0.0046914
1.8308	1.602	0.19993	0.98007	1.0051	0.99838	$5.4332\times {10}^{-10}$	0.00039664
1.5853	0.73703	1.4314	0.9837	1.0037	0.99902	$3.2042\times {10}^{-10}$	0.00029842
0.082806	1.6892	1.1818	1.0041	0.99907	1.0003	$1.7019\times {10}^{-11}$	6.6238e-05
0.2428	0.54946	1.1455	1.0075	0.99842	1.0004	$5.5279\times {10}^{-11}$	0.00011151
0.24627	0.91618	0.86393	1.0279	0.99108	1.0038	$1.5989\times {10}^{-9}$	0.00061966
1.0759	0.84762	1.5499	1.1222	0.97718	1.0042	$1.438\times {10}^{-8}$	0.0017248
1.054	0.17418	0.85505	1.1964	0.96141	1.0082	$3.6805\times {10}^{-8}$	0.0028801
0.22083	1.9961	0.37794	1.0108	0.99797	1.0004	$1.1477\times {10}^{-10}$	0.00014608
0.065989	0.92182	1.351	1.0119	0.9973	1.0008	$1.4633\times {10}^{-10}$	0.00018962

,

Table 7. Results for noisy data for $\tilde{d}\left(u\right)=u^2+u+1$.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error
1.8143	1.5688	1.1067	0.87821	1.0042	1.0132	0.0020016	0.0015718
1.2435	1.4694	0.038102	1.0871	0.96223	1.0227	0.0020017	0.0023355
1.9293	0.9881	1.0046	0.86725	1.0068	1.0123	0.0020016	0.00169
0.8034	1.1622	1.8173	0.81146	1.0196	1.0086	0.0020016	0.002464
1.3517	1.1656	1.2599	1.3018	0.9178	1.0335	0.0020018	0.0050523
1.8308	1.602	0.19993	0.86061	1.0086	1.0118	0.0020016	0.0018171
1.5853	0.73703	1.4314	0.87637	1.0045	1.0131	0.0020016	0.001589
0.082806	1.6892	1.1818	0.88511	1.0027	1.0136	0.0020016	0.001566
0.2428	0.54946	1.1455	0.37128	1.1072	0.98935	0.0020019	0.0089087
0.24627	0.91618	0.86393	0.89739	0.99902	1.0151	0.0020016	0.0016656
1.0759	0.84762	1.5499	1.1151	0.95842	1.0225	0.0020017	0.0023623
1.054	0.17418	0.85505	1.1896	0.94251	1.0266	0.0020017	0.0032631
0.22083	1.9961	0.37794	0.89089	1.0018	1.0136	0.0020016	0.0015688
0.065989	0.92182	1.351	0.91539	0.99126	1.0192	0.0020016	0.0019827

,

Table 8. Results for noisy data for $\tilde{d}\left(u\right)=u^2+u+1$ with regularization parameter.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error	$\mathit{\lambda }$
0.8034	1.1622	1.8173	0.79386	1.0152	1.0136	0.0020043	0.0023183	${10}^{-6}$
1.3517	1.1656	1.2599	0.23178	1.0905	1.0144	0.0020249	0.0095046	${10}^{-5}$
0.22083	1.9961	0.37794	0.15297	0.96968	1.0883	0.0022225	0.025866	${10}^{-4}$
1.3517	1.1656	1.2599	0.11759	0.81796	1.0163	0.0039331	0.1176	${10}^{-3}$

,

Table 9. Results for noise-free data for $\tilde{d}\left(u\right)=e^u$.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error
0.90181	1.5752	1.4868	0.79218	0.94905	1.0105	$5.6575\times {10}^{-8}$	0.0046883
0.57597	1.0598	0.56414	0.54444	0.99668	1.0006	$1.4739\times {10}^{-10}$	0.000262
0.73992	0.76522	1.4468	0.78835	0.95044	1.0099	$5.4814\times {10}^{-8}$	0.0047249
0.62332	0.64684	0.69365	0.67086	0.97351	1.0048	$1.5113\times {10}^{-8}$	0.0024087
0.26005	1.0154	0.48923	0.25222	1.0548	0.98759	$7.5481\times {10}^{-8}$	0.0059219
0.45035	0.83101	0.20517	0.45166	1.0131	0.99813	$7.8171\times {10}^{-9}$	0.0022299
1.4027	1.6491	1.6443	0.53726	0.99922	0.99945	$6.6822\times {10}^{-11}$	0.00026575
1.4448	1.7317	0.62139	0.51571	1.0057	0.99697	$1.064\times {10}^{-9}$	0.00053638
−0.40005	1.2963	0.060998	0.5487	0.99689	1	$9.9924\times {10}^{-11}$	0.00019071
0.13075	0.25531	1.8759	0.32007	1.0428	0.98936	$4.4431\times {10}^{-8}$	0.0045091
0.88974	1.1835	0.41296	0.81012	0.94438	1.0121	$6.5201\times {10}^{-8}$	0.0048764
0.40355	1.7802	1.4709	0.28431	1.0486	0.98887	$5.966\times {10}^{-8}$	0.0052538
−0.44205	0.070614	1.8443	0.55639	0.99609	0.99981	$2.7998\times {10}^{-10}$	0.00020745
−0.45613	0.22429	0.80524	0.55508	0.99551	1.0004	$2.3236\times {10}^{-10}$	0.00026031

,

Table 10. Results for noisy data for $\tilde{d}\left(u\right)=e^u$.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error
0.90181	1.5752	1.4868	0.78517	0.93079	1.0288	0.0020018	0.0050131
0.57597	1.0598	0.56414	0.53751	0.97834	1.019	0.0020016	0.0022175
0.73992	0.76522	1.4468	0.78169	0.93207	1.0283	0.0020018	0.0051061
0.62332	0.64684	0.69365	0.66409	0.9552	1.0231	0.0020017	0.0028501
0.26005	1.0154	0.48923	0.24532	1.0365	1.0059	0.0020017	0.0055843
0.45035	0.83101	0.20517	0.44478	0.99464	1.0166	0.0020016	0.0023482
1.4027	1.6491	1.6443	0.40493	1.0054	1.0124	0.0020016	0.0024947
1.4448	1.7317	0.62139	0.41451	1.0035	1.013	0.0020016	0.0022539
−0.40005	1.2963	0.060998	0.42752	1.001	1.0135	0.0020016	0.0021903
0.13075	0.25531	1.8759	0.3134	1.0244	1.0078	0.0020016	0.0041734
0.88974	1.1835	0.41296	0.80315	0.92604	1.0305	0.0020018	0.0051862
0.40355	1.7802	1.4709	0.27712	1.0303	1.0072	0.0020017	0.0049314
−0.44205	0.070614	1.8443	0.43707	0.99816	1.0147	0.0020016	0.0022583
−0.45613	0.22429	0.80524	0.43069	1.0003	1.0138	0.0020016	0.002201

and

Table 11. Results for noisy data for $\tilde{d}\left(u\right)=e^u$ with regularization parameters.

	${\mathit{d}}_{\mathbf{initial}}$			${\mathit{d}}_{\mathbf{final}}$		$\mathit{I}\left({\mathit{d}}_{\mathbf{final}}\right)$	Relative Error	$\mathit{\lambda }$
0.57597	1.0598	0.56414	0.52797	0.97607	1.0214	0.0020039	0.0023258	${10}^{-6}$
0.26005	1.0154	0.48923	0.22214	1.0051	1.0289	0.002023	0.012937	${10}^{-5}$
0.98925	0.54908	0.46717	0.13735	0.92151	1.0695	0.002208	0.051983	${10}^{-4}$
0.36803	1.0811	0.76951	0.10825	0.7698	0.96251	0.003772	0.27405	${10}^{-3}$

, we provide the value of $I\left(d\right)$ for the final d, i.e., $d_{\mathrm{final}}$ calculated by the algorithm. We also provide the relative error between the solution of the direct problem for the calculated value of $d_{\mathrm{final}}$ and $u=t{x}^2{y}^2$. The relative error is given by

$$\frac{\left|\right|u_{\mathrm{final}}(x,y,t)-t{x}^2{y}^2{\left|\right|}_{\infty }}{\left|\right|t{x}^2{y}^2{\left|\right|}_{\infty }},$$

where $\left|\right|·{\left|\right|}_{\infty }$ is estimated by the maximum value of the function on the meshgrid on $[0,1]\times [0,1]\times [0,0.3]$.

The inverse problem is an ill-posed problem and it is sensitive to the noise. To deal with the noisy data, we use Tikhonov regularization in the error functional and define it as:

$$I\left(d\right):={\int }_0^T{\int }_0^1|u(x,y,t;d)-\tilde{u}(x,1,t){|}^{2}dxdt+{\int }_0^T{\int }_0^1|u(x,y,t;d)-\tilde{u}(1,y,t){|}^{2}dydt+{\lambda \left|\right|d\left|\right|}^2,$$

(38)

where $\left|\right|·\left|\right|$ is the Euclidean norm. We run the algorithm for different regularization parameter $\lambda$ with the same initial values and present the best results according to the relative error.

Experiment 1.

The correct $d(·)$ is $\tilde{d}\left(u\right)=2sin\left(u\right)+cos\left(u\right)$, whose second degree Taylor polynomial is $-0.5u^2+2u+1$. Table 3 and Table 4 show the results for the noise-free and noisy data, respectively. Table 5 shows the best results for different values of the regularization parameter $\lambda$.

Experiment 2.

The correct $d(·)$ is $\tilde{d}\left(u\right)=u^2+u+1$. Table 6 and Table 7 show the results for the noise-free and noisy data, respectively. Table 8 shows the best results for different values of the regularization parameter $\lambda$.

Experiment 3.

The correct $d(·)$ is $\tilde{d}\left(u\right)=e^u$ whose second degree Taylor polynomial is $0.5u^2+u+1$. Table 9 and Table 10 show the results for the noise-free and noisy data, respectively. Table 11 shows the best results for different values of the regularization parameter $\lambda$.

Since the algorithm is based on approaching $d\left(u\right)$ with polynomials and the error function is a function of three variables, the second degree Taylor polynomials of the correct $d\left(u\right)$’s are expected to be attained by the algorithm for each initial point.

In all experiments with noise-free data, we notice that the algorithm finds the linear coefficients of the target Taylor polynomial for almost all initial points, while it mostly fails to reach the nonlinear coefficient, i.e., the coefficient of $u^2$ of the target Taylor polynomial for most of the initial points. However, the coefficients of the final $d\left(u\right)$’s with the least relative errors almost match the coefficient of $u^2$ of the target Taylor polynomial. In all of the experiments, the relative error for the solution u corresponding to each final polynomial $d\left(u\right)$ is observed to fluctuate between 0.0001 and 0.007. Another observation for noise-free experiments is that the relative error for the corresponding u increases with $I\left(d\right)$.

In the experiments that are carried out with noisy data, there is a remarkable behavior in the error function: $I\left(d_{\mathrm{final}}\right)$ stays around 0.002 for almost all initial points. The relative error for all initial points and experiments fluctuate between 0.001 and 0.009. So the relative error increases up to 10 times. One observation about the experiments with noisy data is that even though the linear coefficients of $d_{\mathrm{final}}$ are more distanced from the linear coefficients of the expected Taylor polynomials, they still seem to accumulate around the expected linear coefficients. The nonlinear coefficients look far from the expected.

The experiments show that the linear part of the $d\left(u\right)$ is stable with respect to the noise. This is most probably due to the constraints imposed by the method used in the solution of the direct problem. Note that we have a stable solution for t in $[0,0.3]$ with $▵={10}^{-5}$.

For the noisy data, we apply the Tikhonov regularization. The regularization parameters are very close to zero, i.e., the values that are less than ${10}^{-7}$ results in reliable results. In addition, it has been observed that the linearization of $d\left(u\right)$ is robust against the noise. It appears that the robustness of the linear part increases as $\beta$ approaches to 1. In the experiments with $\beta =0.6$, (not shown in this article) the results have shown that the error is higher and that the final points are found to be beyond the expected. So, we can conclude that the problem becomes more ill-posed as $\beta$ decreases.

We also note that one can also use trust region methods to minimize the error functional. Indeed, some preliminary tests have been carried out with trust region method and the BGFS method. A significant difference in the computation time was observed between two optimization methods. For this reason, the BGFS method was chosen as the primary algorithm for this study.

6. Concluding Remarks

We study an inverse problem for the nonlinear time fractional equation $u_t^{\beta }=\nabla ·(d\left(u\right)\nabla u)+f(x,y,t)$. We solve the direct problem by converting it into a system of fractional ODEs using the method of lines and vectorization. In doing so, we adopted the operator approach to approximate the derivatives that reduces computational and memory demand of the algorithm. After we solve the direct problem, we reformulate the inverse problem as a minimization problem. Then we solve the minimization problem by using the BGFS method and a Quasi-Newton method. The inverse problem is solved for both noise-free and noisy data. When solving the direct problem, we observe that as the order of fractional time derivative approaches to the value 1, the numerical solution approximates the exact solution in much faster and reliable manner. Additionally, it appears that there is a subtle relation between the size of the interval and time steps used in the numerical scheme. We observe that $M=N=5$ provides the best performance in terms of the minimum simulation time and the absolute error. The numerical results of the experiments show that solving the inverse problem that we examined in this article using the error functional with the help of polynomials and solving the direct problem with the method of lines work well together. However, the computational constraints resulting from the discussion of the solution of the direct problem makes the problem a bit restrictive on the meshgrid and time interval. The experiments with the noise-free data verifies the theoretical results, i.e., the existence of a unique quasi-solution.The experiments with noisy data show that the linear part of the $d\left(u\right)$ shows a robust behavior against the noise. Finally, the values of $\beta$ away from 1 results in unstable results, that is, the problem becomes more ill-posed.

The authors of this paper plan to study a set of three partial differential equations involving time-fractional derivatives and nonlinear diffusion operators. After we prove that the direct problem is a well-posed problem, we define the inverse problem. The inverse problem may be either an inverse parameter problem, or an inverse coefficient problem, or an inverse source problem. Then we study the existence and uniqueness of the solution and solve the inverse problem by the method presented in this paper.