Anomalous Diffusion with Caputo-Fabrizio Time Derivative: an Inverse Problem

ABSTRACT

In this work we identify the source in a 1D anomalous diffusion equation, from measurements of the concentration at a finite number of points. We use Caputo-Fabrizio time fractional derivative to model the phenomenon. Separating variables, we arrive to a linear system which provides approximate values for the Fourier coefficients of the unknown source. Numerical examples show the efficiency of the method, as well as some of its practical limitations.

Keywords:
inverse problems; fractional calculus; anomalous diffusion

1 INTRODUCTION

The motion of microscopic particles in a fluid has been investigated for a long time. The “Brownian motion” (named after Robert Brown, 1827) consists of random displacements of the suspended particles. The molecules of the fluid, much smaller than the solid particles, knock them incessantly, both driving and damping their movement. Macroscopically, those knockings give rise to the “fluid viscosity”.

The traditional mathematical model for this diffusion phenomenon is based on Einstein’s theory by which the mean square displacement of a diffusion particle is proportional to time. This idea leads to the classical diffusion equation,

$u{\text{'}}_t(\mathbf{\text{x}},t)-k{\nabla }^{2}u(\mathbf{\text{x}},t)=0,$

where u(x ,t) is the quantity of particles per unit volume at position x and time t, and k is a diffusion coefficient. If there is also a source of additional incoming particles, the equation becomes non homogeneous:

$u{\text{'}}_t(\mathbf{\text{x}},t)-k{\nabla }^{2}u(\mathbf{\text{x}},t)=f(\mathbf{\text{x}},t).$

This equation may be derived by a statistical reasoning, or else, by macroscopic considerations of conservation of matter and Fick’s experimental law (“the net flow of diffusion particles occurs from higher to lower concentrations, and is proportional to the concentration gradient”), which is similar to Fourier’s heat law.

However, there are experimental results - like anomalous diffusion of particles in porous or fractal media, in biological media, in turbulent plasma, in polymers, etc. - that show that, in some cases, the mean square displacement of the particles must be considered to be proportional not to the time but to a fractional power of time, to fit the empirical data. This fractional order may be less than unity (subdiffusion) or greater than one (superdiffusion).

Fractional differential equations have been proposed to model this anomalous diffusion phenomena:

${\mathcal{D}}_t^{\alpha }u(\mathbf{\text{x}},t)-k{\nabla }^{2}u(\mathbf{\text{x}},t)=f(\mathbf{\text{x}},t),$ (1.1)

where α is a not integer order of derivation with respect to time; 0 < α < 1 for subdiffusion, α > 1 for superdiffusion.

There are several definitions of the fractional derivative 𝒟_t ^α: Riemann-Liouville’s, Caputo’s, Atangana-Banleanu’s, Caputo-Fabrizio’s, among others. These derivatives are defined as integral operators which take account of the past “history” of the function to which they are applied (see ³⁴ for an exhaustive list). In the last decades fractional calculus has been successfully used to model phenomena in different areas: diffusion problems, hydraulics, potential theory, control theory, electrochemistry, electromagnetism, viscoelasticity and nanotechnology (see for example ^{1), (2), (6), (7), (8), (9), (10), (16), (22), (28), (30), (38} among others). Numerical schemes to calculate approximate solutions to fractional differential equations have also been introduced (see ^{5), (24), (35), (39), (40}).

In this work we have chosen the Caputo-Fabrizio fractional derivative ¹² to describe a 1-dimensional model of subdiffusion of the form

${\mathcal{D}}_t^{\alpha }u(x,t)-k\frac{{\partial }^2}{\partial x^2}u(x,t)=s\left(x\right)h\left(t\right),$

where $x\in (0,1)$ and $t\in (0,T)$.

Note that the source, f(x,t), is supposed to be a separable function of variables x and $t:f(x,t)=s\left(x\right)h\left(t\right)$.

In this context we pose the inverse problem that consists in finding s(x) and u(x,t), for a known h(t), with additional data: measurements at the final time T, u(x _i , T), with $x_i\in (0,1)$ for i = 1,· · · , N.

The Caputo-Fabrizio fractional derivative (CFFD) has been proposed to describe physical phenomena in different fields as thermodynamics, electromagnetism and continuum mechanics (see, for example, ^{13),(14),(17),(33}). The reason for choosing this fractional derivative is that it has a non- singular kernel (with the benefits that this property implies) while maintaining the appropriate mathematical and memory properties to describe these kind of dynamics. However, it is necessary to point out that the definition of CFFD imposes limitations to the boundary conditions; we will describe some of them briefly in the next section.

This work is organized as follows: in the next section we present some mathematical definitions and properties of the CFFD; the proposed solution to the inverse problem appears in Section 3. In Section 4 we present some numerical examples. Finally, in Section 5, we state some conclusions.

2 MATHEMATICAL PRELIMINARIES

2.1 The Caputo-Fabrizio Derivative

We denote by H ¹(a, b) the Sobolev space W ^1,2(a, b) of functions $f:(a,b)\to \mathrm{ℝ}$, with (weak) derivative $f\text{'}\in L^2(a,b)$.

For $f\in H^1(a,b)$, the Caputo-Fabrizio fractional derivative (CFFD) of order α, with 0 < α < 1, is defined by

${}_a{}^{CF}{\mathcal{D}}_t^{\alpha }\left[f\right]\left(t\right):=\frac{M\left(\alpha \right)}{1-\alpha }{\int }_a^{t}f\text{'}\left(s\right)e^{-\frac{\alpha (t-s)}{1-\alpha }}ds$ (2.1)

where $-\infty \le a<b$ and M(α) is a normalizing positive factor verifying M(0) = M(1) = 1.

Remark: Observe ${}_a{}^{CF}{\mathcal{D}}_t^{\alpha }\left[f\right]\left(a\right)=0$. ◇

2.2 Solutions of two elementary inverse problems with CFFD

In this subsection we solve two simple fractional inverse problems that will be useful to approximate the solution to the inverse source problem in which we are interested.

• Let us consider the following basic inverse problem involving CFFD:

$\left\{\begin{array}{c}{}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[f\right]\left(t\right)=h\left(t\right)0<t<b\\ f\left(0\right)=f_0,\end{array}\right\}$

for $f\in H^1(0,b)$ unknown, $f_0\in \mathrm{ℝ}$ and h(t) a known function, continuous for $t\in (0,b)$, with $\underset{t\to 0^{+}}{\lim }h\left(t\right)=0$.

From

$\frac{M\left(\alpha \right)}{1-\alpha }{\int }_0^{t}f\text{'}\left(s\right)e^{-\frac{\alpha (t-s)}{1-\alpha }}ds=h\left(t\right)$

we have

$\begin{array}{c}{\int }_0^{t}f\text{'}\left(s\right)e^{\frac{\alpha s}{1-\alpha }}ds=\frac{1-\alpha }{M\left(\alpha \right)}{e}^{\frac{\alpha t}{1-\alpha }}h\left(t\right).\end{array}$

As $f\text{'}\left(s\right)e^{\frac{\alpha s}{1-\alpha }}$ is integrable on (0, b), $\frac{1-\alpha }{M\left(\alpha \right)}{e}^{\frac{\alpha t}{1-\alpha }}h\left(t\right)$ is absolutely continuous in (0, b) and, consequently, h(t) is also absolutely continuous in (0, b).

Then, almost for every $t\in (0,b)$ where derivatives exist,

$\begin{array}{c}f\text{'}\left(t\right)e^{\frac{\alpha t}{1-\alpha }}=\frac{1-\alpha }{M\left(\alpha \right)}\left[\frac{\alpha }{1-\alpha }{e}^{\frac{\alpha t}{1-\alpha }}h\left(t\right)+e^{\frac{\alpha t}{1-\alpha }}h\text{'}\left(t\right)\right],\end{array}$

or

$\begin{array}{c}f\text{'}\left(t\right)=\frac{1-\alpha }{M\left(\alpha \right)}\left[\frac{\alpha }{1-\alpha }h\left(t\right)+h\text{'}\left(t\right)\right].\end{array}$

Then

$\begin{array}{c}f\left(t\right)-f\left(0\right)=\frac{\alpha }{M\left(\alpha \right)}{\int }_0^{t}h\left(s\right)ds+\frac{1-\alpha }{M\left(\alpha \right)}h\left(t\right),\end{array}$

and

$\begin{array}{c}f\left(t\right)=f_0+\frac{\alpha }{M\left(\alpha \right)}{\int }_0^{t}h\left(s\right)ds+\frac{1-\alpha }{M\left(\alpha \right)}h\left(t\right).\end{array}$ (2.2)

• Now let us consider this another simple inverse problem involving CFFD:

$\left\{\begin{array}{c}{}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[f\right]\left(t\right)+\lambda f\left(t\right)=g\left(t\right)0<t<b\\ f\left(0\right)=f_0,\end{array}\right.$ (2.3)

for $f\in H^1(0,b)$ unknown, f ₀ and λ in ℝ and g(t) a known function, continuous for $t\in (0,b)$.

Remark: Note that, since ${}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[f\right]\left(0\right)=0$, it must be λ f₀ = g(0) for compatibility reasons. ◇

To solve (2.3), if $h\left(t\right)=g\left(t\right)-\lambda f\left(t\right)$ in (2.2), we have

$\begin{array}{c}f\left(t\right)=\underset{\psi \left(t\right)}{\underbrace{f_0+\frac{\alpha }{M\left(\alpha \right)}{\int }_0^{t}g\left(s\right)ds+\frac{1-\alpha }{M\left(\alpha \right)}g\left(t\right)}}-\lambda \left[\frac{\alpha }{M\left(\alpha \right)}{\int }_0^{t}f\left(s\right)ds+\frac{(1-\alpha )}{M\left(\alpha \right)}f\left(t\right)\right].\end{array}$

Then, for every $t\in (0,b)$ where derivatives exist,

$\begin{array}{c}f\text{'}\left(t\right)=\psi \text{'}\left(t\right)-\lambda \left[\frac{\alpha }{M\left(\alpha \right)}f\left(t\right)+\frac{(1-\alpha )}{M\left(\alpha \right)}f\text{'}\left(t\right)\right].\end{array}$

• - For $\lambda =\frac{M\left(\alpha \right)}{\alpha -1}$ ,

$\begin{array}{r}f\left(t\right)=\frac{\alpha -1}{\alpha }\psi \text{'}\left(t\right)=-\frac{1-\alpha }{M\left(\alpha \right)}\left[g\left(t\right)+\frac{1-\alpha }{\alpha }g\text{'}\left(t\right)\right].\end{array}$

• - If $\lambda \ne \frac{M\left(\alpha \right)}{\alpha -1}$ ,

$f\text{'}\left(t\right)+\frac{\lambda \alpha }{M\left(\alpha \right)+\lambda (1-\alpha )}f\left(t\right)=\frac{M\left(\alpha \right)}{M\left(\alpha \right)+\lambda (1-\alpha )}\psi \text{'}\left(t\right).$

Denoting ${\delta }_{\lambda }=M\left(\alpha \right)+\lambda (1-\alpha )$, (2.3) is equivalent to

$\left\{\begin{array}{rcl}f\text{'}\left(t\right)+\frac{\lambda \alpha }{{\delta }_{\lambda }}f\left(t\right)& =& \frac{1}{{\delta }_{\lambda }}\left[\alpha g\right(t)+(1-\alpha )g\text{'}(t\left)\right]\\ f\left(0\right)=f_0.& & \end{array}\right.$ (2.4)

Its solution and, consequently, the solution to (2.3) is

$\begin{array}{r}f\left(t\right)=f_0{e}^{-\frac{\lambda \alpha t}{{\delta }_{\lambda }}}+\frac{\alpha }{({\delta }_{\lambda }{)}^2}{\int }_0^{t}g\left(s\right)e^{-\frac{\lambda \alpha (t-s)}{{\delta }_{\lambda }}}ds\\ +\frac{1-\alpha }{{\delta }_{\lambda }}\left[g\left(t\right)-g\left(0\right)e^{-\frac{\lambda \alpha t}{{\delta }_{\lambda }}}\right]\end{array}$ (2.5)

(in ³² the author arrives to the same expression, when f ₀ = 0, using another procedure).

Remark: Note that, if g ≡ 0, we are looking for the eigenfunctions of CFFD, i.e., $f\in H^1(0,b)$ such that ${}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[f\right]\left(t\right)=-\lambda f\left(t\right)$, but this operator has no eigenfunctions. Actually, the equivalent problem (2.4) becomes

$\left\{\begin{array}{r}f\text{'}\left(t\right)+\frac{\lambda \alpha }{{\delta }_{\lambda }}f\left(t\right)=0\\ f\left(0\right)=0\end{array}\right.$

which has only the trivial solution. ◇

3 APPROXIMATE SOLUTION TO AN INVERSE PROBLEM FOR A FRACTIONAL DIFFUSION EQUATION

Identification of sources in diffusion equations appear in mathematical modelling of different processes in diverse areas of science and engineering. The objective is to identify unknown sources from measurable boundary and / or final output data. Uniqueness of the source can only be guaranteed for certain types of sources, depending on the equations that govern the process ⁴. Analytical methods, as well as accurate numerical schemes, have been proposed to efficiently solve these inverse problems.

Inverse source problems for different type of sources appeared in ^{11),(18),(20),(31),(37}. In ²⁷ there is a complete list of references where identification of various type of sources are listed. The case of fractional diffusion equation was studied in ^{3), (15), (25), (26), (29), (36}, among others.

Separable source terms of the form h(t)s(x) arise in various physical model; for example, the

heat process of radioactive decay ¹⁹, and have been extensively studied. In ²¹ and ²³, identification of this type or sources, for the heat conduction equation, is studied.

In this work we will consider the following boundary value problem of anomalous diffusion in one dimension:

$\left\{\begin{array}{r}{}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[u\right](x,t)-u\text{'}{\text{'}}_{xx}(x,t)=s\left(x\right)h\left(t\right)0<t<T,0<x<1\\ u(0,t)=u(1,t)=00\le t\le T\end{array}\right.$

where ${}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[u\right]$ stands for the temporal CFFD of order α, with 0 < α < 1.

This equation attempts to model the diffusion of particles in a linear path of length 1, constituted by an heterogeneous medium which behaves as having a sort of “memory”, due to the fluctuations introduced by elements at different dimension scales ¹². The product s(x)h(t), with $s\in C[0,1]$ and $h\in C[0,T]$, represents an external source of particles; only the temporal factor, h(t), is known.

We look for smooth solutions $u\in C^2\left(\right[0,1]\times [0,T\left]\right)$ and s(x), the spatial component of the source. We have some additional data: measurements u(x _i , T) in several positions $x_i\in [0,1]$, for i = 1, · · · , N, at final time T.

Similarly to what happens in the standard case, this equation can be solved by separating variables: we propose

$u(x,t)=\sum _{k\in \mathrm{\mathbb{N}}}{u}_k\left(t\right)\sin \left(k\pi x\right),$ (3.1)

with u _k (t) the Fourier coefficients of u(x,t) for each $t\in [0,T]$.

Remark: If $u\in C^2\left(\right[0,1]\times [0,T\left]\right)$, the derivatives $u\text{'}{\text{'}}_{xx}(x,t)$ and ${}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[u\right](x,t)=\frac{{\displaystyle M\left(\alpha \right)}}{{\displaystyle 1-\alpha }}{\int }_0^t{u}_t(x,s)e^{-\frac{{\displaystyle \alpha (t-s)}}{{\displaystyle 1-\alpha }}}ds$ are in $C\left(\right[0,1]\times [0,T\left]\right)$. If $u_k^{**}\left(t\right)$ and $u_k^{*}\left(t\right)$ are, respectively, the Fourier coefficients of $u\text{'}{\text{'}}_{xx}(x,t)$ and ${}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[u\right](x,t)$ for each $t\in [0,T]$, it is easy to prove that $u_k^{**}\left(t\right)=-k^2{\pi }^2{u}_k\left(t\right)$ and $u_k^{*}\left(t\right){=}_0^{CF}{\mathcal{D}}_t^{\alpha }\left[u_k\right]\left(t\right)$. ◇

As ${}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[u\right](x,t)-u\text{'}{\text{'}}_{xx}(x,t)=s\left(x\right)h\left(t\right)$, and assuming that (3.1) is an absolutely convergent series, the Fourier coefficients of u _k (t) verify

$\begin{array}{rl}{}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[u_k\right]\left(t\right)+k^2{\pi }^2{u}_k\left(t\right)& =h\left(t\right)s_k\forall t\in (0,T)\end{array}$ (3.2)

with s _k the Fourier coefficients of $s\left(x\right):s_k=2{\int }_0^{1}s\left(x\right)\sin \left(k\pi x\right)dx$.

Note that, as ${}_0{}^{CF}{\mathcal{D}}_t^{\alpha }\left[u_k\right]\left(0\right)=0$, it must be $u_k\left(0\right)=\frac{{\displaystyle h\left(0\right)s_k}}{{\displaystyle (k\pi {)}^2}}$.

Now, from (2.5), if ${\delta }_k=M\left(\alpha \right)+k^2{\pi }^2(1-\alpha )$,

$\begin{array}{r}{u}_k\left(t\right)=s_k\{\frac{h\left(0\right)}{k^2{\pi }^2}{e}^{-\frac{k^2{\pi }^2\alpha t}{{\delta }_k}}+\frac{\alpha }{({\delta }_k{)}^2}{\int }_0^{t}h(s)e^{-\frac{k^2{\pi }^2\alpha (t-s)}{{\delta }_k}}ds+\\ +\frac{1-\alpha }{{\delta }_k}\left[h\right(t)-h(0\left)e^{-\frac{k^2{\pi }^2\alpha t}{{\delta }_k}}\right]\}\end{array}$ (3.3)

Remark: From (3.3), u(x,t) in (3.1) is an absolutely and uniformly convergent series. Since $h\in C[0,T],\left|h\right(t\left)\right|\le K$, for certain $K\in {\mathrm{ℝ}}_{>0}$ and for all $t\in [0,T]$, we have

$\begin{array}{l}\left|\sum _{k\in \mathrm{\mathbb{N}}}{u}_k\right(t\left)\sin \right(k\pi x){|}^2\le \\ \\ \le \left(\sum _{k\in \mathrm{\mathbb{N}}}\right|s_k{|}^2\left)\right(\sum _{k\in \mathrm{\mathbb{N}}}\left|\right\{\frac{h\left(0\right)}{k^2{\pi }^2}{e}^{-\frac{k^2{\pi }^2\alpha t}{{\delta }_k}}+\frac{\alpha }{({\delta }_k{)}^2}{\int }_0^{t}h\left(s\right)e^{-\frac{k^2{\pi }^2\alpha (t-s)}{{\delta }_k}}ds+\\ \\ +\frac{1-\alpha }{{\delta }_k}\left[h\right(t)-h(0\left)e^{-\frac{k^2{\pi }^2\alpha t}{{\delta }_k}}\right]\left\}{|}^2\right)\le \\ \\ \le \frac{K^2}{{\pi }^4}\left(\sum _{k\in \mathrm{\mathbb{N}}}\right|s_k{|}^2\left)\right(\sum _{k\in \mathrm{\mathbb{N}}}|\frac{3}{k^2}+\frac{\alpha T}{k^4{\pi }^2(1-\alpha {)}^2}{|}^2).\end{array}$

These last two series are absolutely convergent. ◇

From (3.1) and (3.3), for any position $x_i\in [0,1]$, at time T,

$\begin{array}{l}u(x_i,T)=\sum _{k\in \mathrm{\mathbb{N}}}{s}_k\{\frac{h\left(0\right)}{k^2{\pi }^2}{e}^{-\frac{k^2{\pi }^2\alpha T}{{\delta }_k}}+\frac{\alpha }{({\delta }_k{)}^2}{\int }_0^{T}h(s)e^{-\frac{k^2{\pi }^2\alpha (T-s)}{{\delta }_k}}ds+\\ +\frac{1-\alpha }{{\delta }_k}\left[h\right(T)-h(0\left)e^{-\frac{k^2{\pi }^2\alpha T}{{\delta }_k}}\right]\left\}\sin \right(k\pi x_i)\end{array}$ (3.4)

Suppose we know h(t) and we measure the concentration u(x, T ) at M points of [0, 1], {x ₁ , x ₂ , ..., x _M }; let us denote u _i = u(x _i , T).

We want to obtain an approximation to the function s(x), the spatial component of the source.

From (3.4) we can propose a linear system to compute the Fourier coefficients s _k of s(x), for k = 1, ..., N. Indeed, from

$$\begin{gathered} u_i\cong \sum _{k=1}^N{s}_k\{\frac{{\displaystyle h\left(0\right)}}{{\displaystyle k^2{\pi }^2}}{e}^{-\frac{{\displaystyle k^2{\pi }^2\alpha T}}{{\displaystyle {\delta }_k}}}+\frac{{\displaystyle \alpha }}{{\displaystyle ({\delta }_k{)}^2}}{\int }_0^{T}h(s)e^{-\frac{{\displaystyle k^2{\pi }^2\alpha (T-s)}}{{\displaystyle {\delta }_k}}}ds+ \\ +\frac{{\displaystyle 1-\alpha }}{{\displaystyle {\delta }_k}}[h\left(T\right)-h\left(0\right)e^{-\frac{{\displaystyle k^2{\pi }^2\alpha T}}{{\displaystyle {\delta }_k}}}\left]\right\}\sin \left(k\pi x_i\right), \end{gathered}$$

we have

$\left[\begin{array}{c}{u}_1\\ u_2\\ ...\\ u_M\end{array}\right]=\left[\begin{array}{cccc}{a}_{11}& a_{12}& ...& a_{1N}\\ a_{21}& a_{22}& ...& a_{2N}\\ ...& & & \\ a_{M1}& a_{M2}& ...& a_{MN}\end{array}\right]\left[\begin{array}{c}{s}_1\\ s_2\\ ...\\ s_N\end{array}\right]$ (3.5)

for

$$\begin{gathered} a_{ik}=\{\frac{{\displaystyle h\left(0\right)}}{{\displaystyle k^2{\pi }^2}}{e}^{-\frac{{\displaystyle k^2{\pi }^2\alpha T}}{{\displaystyle {\delta }_k}}}+\frac{{\displaystyle \alpha }}{{\displaystyle ({\delta }_k{)}^2}}{\int }_0^{T}h(s)e^{-\frac{{\displaystyle k^2{\pi }^2\alpha (T-s)}}{{\displaystyle {\delta }_k}}}ds+ \\ +\frac{{\displaystyle 1-\alpha }}{{\displaystyle {\delta }_k}}[h\left(T\right)-h\left(0\right)e^{-\frac{{\displaystyle k^2{\pi }^2\alpha T}}{{\displaystyle {\delta }_k}}}\left]\right\}\sin \left(k\pi x_i\right). \end{gathered}$$

If M = N and $A=\left(a_{ik}\right)\in {\mathrm{ℝ}}^{M\times N}$ is an invertible matrix the solution is straightforward, at least theoretically; if $M\ne N$ the corresponding least squares problem could be posed.

Once the Fourier coefficients s _k are estimated, from (3.3) and (3.1) the function u(x,t) can be approximated too.

4 NUMERICAL EXAMPLES

Clearly, the more Fourier coefficients we can evaluate, the better the approximation of s(x) we get. But, for reasons of computation economy, if the difference between two approximate solutions, calculated with N and $N+1$ Fourier terms respectively, is small enough, within our error expectation, we might decide to stop applying the algorithm.

We chose examples of increasing difficulty in order to evaluate the performance of the numerical scheme.

Problems with known solution u(x,t) were selected. All the tests were performed with simulated noisy measurements at equally spaced points in [0, 1], by assigning a normal distribution to the experimental error, with a variance of 1% around exact values of u(x _i , T).

4.1 Example 1

Let us consider the following “simple” diffusion problem with known solution:

$\left\{\begin{array}{c}\begin{array}{rcl}{}_0{}^{CF}{\mathcal{D}}_t^{0.5}\left[u\right](x,t)-u\text{'}{\text{'}}_{xx}(x,t)& =& t^{3}s\left(x\right),\forall (x,t)\in (0,1)\times (0,1)\end{array}\\ u(0,t)=u(1,t)=0\forall t\in [0,1]\end{array}\right.$ (4.1)

In this case h(t) = t ³ and α = 0.5.

For s(x) = 2 sin(4πx) the exact solution is

$u(x,t)=\frac{{\displaystyle 3(1+8{\pi }^2{)}^2{e}^{\frac{{\displaystyle -8{\pi }^{2}t}}{{\displaystyle 1+8{\pi }^2}}}-3+24{\pi }^2(t-2)+256{\pi }^6{t}^3-96{\pi }^4(t^2-2t+2)}}{{\displaystyle 2048{\pi }^8}}\sin \left(4\pi x\right).$

Assuming neither u(x,t) nor s(x) are known, we simulated values of u(x _i ,t), endowed with some experimental error, at six equally spaced measuring points x _i in [0, 1]. We considered six Fourier coefficients of s(x) (so N = 6) and approximated them by means of (3.5) with M = 6. An invertible matrix (a _ik ) was obtained (note that, in this case, only one Fourier coefficient is needed to obtain the exact known solution).

Figure 1 shows sexact(x) (dashed in red) and sapprox(x) (in green). In Figure 2, uexact(x,t) (in red) and uapprox(x,t) (in green). In Figure 3, a detail of the graph of sexact(x) (dashed in red) and sapprox(x) (in different colors) for various simulations for problem (4.1).

Figure 1:
s_exact(x), dashed in red, and s_approx.(x), in green, for a simulation of the inverse problem (4.1).

Figure 2:
u_exact(x,t), in red, and u_approx.(x,t), in green, for a simulation of the inverse problem (4.1). The scale of the vertical axe must be multiplied by 10−2.

Figure 3:
s_exact(x), dashed in red, and s_approx.(x), in different colors: a detail of the plot for various simulations of problem (4.1).

We calculated the relative quadratic error

$\frac{{\displaystyle \left|\right|s_{\text{exact}}-s_{\text{approx.}}|{|}_{L_2(0,1)}}}{{\displaystyle \left|\right|s_{\text{approx.}}|{|}_{L_2(0,1)}}}$

for 50 different realizations of the experiment. The resulting values are between 0.2% and 2% and averaged 1.1%.

4.2 Example 2

For the following diffusion problem,

$\left\{\begin{array}{ccc}\begin{array}{rc}{}_0{}^{CF}{\mathcal{D}}_t^{0.5}\left[u\right](x,t)-u\text{'}{\text{'}}_{xx}(x,t)& =\end{array}& \sin \left(t\right)s\left(x\right),& \forall (x,t)\in (0,1)\times (0,1)\\ \begin{array}{rc}u(0,t)& =\end{array}& u(1,t)=0& \forall t\in [0,1],\end{array}\right.$ (4.2)

it is h(t) = sin(t) and α = 0.5.

For $s\left(x\right)=\sin \left(2\pi x\right)+\sin \left(3\pi x\right)-2\sin \left(6\pi x\right)$ the exact solution is

$$\begin{gathered} u(x,t)=\sum _{j=2,3}\frac{{\displaystyle e^{-\frac{{\displaystyle j^2{\pi }^{2}t}}{{\displaystyle 2+j^2\pi }}}-\cos \left(t\right)+\sin \left(t\right)+j^2{\pi }^2\sin \left(t\right)}}{{\displaystyle 2+2j^2{\pi }^2+j^4{\pi }^4}}\sin \left(j\pi x\right)- \\ -\frac{{\displaystyle e^{-\frac{{\displaystyle 18{\pi }^{2}t}}{{\displaystyle 1+18\pi }}}-\cos \left(t\right)+\sin \left(t\right)+36{\pi }^2\sin \left(t\right)}}{{\displaystyle 1+36{\pi }^2+648{\pi }^4}}\sin \left(6\pi x\right) \end{gathered}$$

Note that, in this case, regarding the known exact solution, three Fourier coefficients would be enough. However, assuming neither u(x,t) nor s(x) are known, we performed 10 sets of simulations for N = 1, 2, ..., 10 equally spaced measuring points in [0, 1]. For each value of N we simulated 50 repetitions of the experiment and obtained s ₁, s ₂, ..., s _N - the approximate Fourier coefficients of s(x) - by means of (3.5); for each repetition we constructed the approximation $s_{{\text{approx.}}_N}\left(x\right)$ consisting of M = N Fourier terms for the spatial part of the source. Then, by a simple average between the 50 simulations, we obtained the corresponding mean function ${\overline{s}}_{{\text{approx.}}_N}\left(x\right)$ for each N.

During the simulations we observed that, by changing the position of the measuring points, the condition number of the matrix could be improved. This issue deserves a deeper analysis that we will carry out in the future.

Looking for a general criteria to decide when to stop the algorithm, for each pair of consecutive approximate mean solutions we obtained the quadratic distances:

$D_{(N,N+1)}=\frac{{\displaystyle \left|\right|{\overline{s}}_{{\text{aprox.}}_{N+1}}-{\overline{s}}_{{\text{approx.}}_N}|{|}_{L_2(0,1)}}}{{\displaystyle \left|\right|{\overline{s}}_{{\text{approx.}}_{N+1}}|{|}_{L_2(0,1)}}}.$

We show the obtained values in Table 1.

Table 1:
Quadratic differences between simulations

N	1	2	3	4	5	6	7	8	9
D(N,N+1)	1.01	0.65	0.55	0.63	0.81	0.0376	0.0382	0.055	0.062

The smaller distance between approximations was for N = 6. The quadratic errors for the corresponding mean functions ${\overline{s}}_{{\text{approx.}}_6}$ and ${\overline{s}}_{{\text{approx.}}_7}$ were

$\frac{\left|\right|s_{\text{exact}}-{\overline{s}}_{{\text{approx.}}_6}|{|}_{L_2(0,1)}}{\left|\right|{\overline{s}}_{{\text{approx.}}_6}|{|}_{L_2(0,1)}}=0.006$

and

$\frac{\left|\right|s_{\text{exact}}-{\overline{s}}_{{\text{approx.}}_7}|{|}_{L_2(0,1)}}{\left|\right|{\overline{s}}_{{\text{approx.}}_7}|{|}_{L_2(0,1)}}=0.002.$

Figure 4 shows s _exact(x) (dashed in red) and ${\overline{s}}_{{\text{approx}}_6}\left(x\right)$ (in green). Figure 5 shows s _exact(x) (dashed in red) and s _approx(x) for N = 1, 2, 3, ..., 10 points of measurements (in different colors). In Figure 6, u _exact(x,t) (in red) and ${\overline{u}}_{{\text{approx}}_7}(x,t)$ (in green).

Figure 4:
s_exact(x), dashed in red, and ${\overline{s}}_{{\text{approx.}}_6}\left(x\right)$, in green, for 50 simulations of the inverse problem (4.2) with 7 equally spaced measuring points.

Figure 5:
s_exact(x) (dashed in red) and s_approx(x) for N = 1, 2, 3, ..., 10 points of measurements (in different colors) for 10 simulations of the inverse problem (4.2).

Figure 6:
u_exact(x), in red, and ${\overline{u}}_{{\text{approx.}}_7}\left(x\right)$, in green for the inverse problem (4.2).

During the implementation of the simulations we observed that the condition number of the matrix A increases as the number of equally spaced measurement points grows. It is evident that, in the general case, a trade-off must be established between the number of measurement points chosen and the condition number of the matrix.

5 CONCLUSIONS

In this work we propose and approximation scheme to find the solution and the spatial component of the source in a non homogeneous fractional diffusion equation in 1D, from measurements of the concentration at some points. By separating variables, we arrive to a simple fractional differential equation that can be solved explicitly. Numerical examples show the good performance of the proposed scheme if the spatial part of the source can be represented by a few number of Fourier terms of low frequency, but the linear algebraic problem become bad conditioned if the number of terms increases. In the general case, the more Fourier coefficients we obtain, the better the function s(x) can be approximated; nevertheless the increasing condition number of the matrix conspires against this objective. The position of the measuring points is one of the subjects to be analyzed as, in the examples with known solution, we observed that the method becomes unstable if we consider an increasing number of equally spaced points. When carrying out the simulations we also observed that the condition number of the matrix A could be improved by varying the spacing. An optimization procedure could be implemented to determine a better choice for the position of the measurement points; it is pending as a future work.

Acknowledgments

This work was partially supported by UBACyT Project 20020170100350BA.

REFERENCES

• ¹
  R. Agarwal, S. Jain & R.P. Agarwal. Mathematical Modeling and Analysis of Dynamics of Cytosolic Calcium Ion in Astrocytes using Fractional Calculus. Journal of Fractional Calculus and Applications, 9(2) (2018), 1-12.
• ²
  E. Ahmed, A. El-Sayed & H. El-Saka. Equilibrium points, stability and numerical solutions of fractional-order predator-prey and rabies models. J. Math. Anal. Appl., 325 (2007), 542-553.
• ³
  M. Ali, S. Aziz & S. Malik. Inverse source problems for a space-time fractional differential equation. Inverse Prob. Sci. Eng., 28(1) (2020), 47-68.
• ⁴
  A.E. Badia, T. Ha-Duong & A. Hamdi. Identification of a point source in a linear advection-dispersion- reaction equation: Application to a pollution source. Inverse Problems, 21 (2005), 1-17.
• ⁵
  D. Baleanu, K. Diethelm, E. Scalas & J. Trujillo. “Fractional Calculus: Models and Numerical Methods”. Series on Complexity, Nonlinearity and Chaos. World Scientific (2012).
• ⁶
  D. Baleanu , A. Golmankhaneh & M. Baleanu. Fractional Electromagnetic Equations Using Fractional Forms. International Journal of Theoretical Physics, 48 (2009), 3114-3123.
• ⁷
  D. Baleanu , Z. Guvenc & J.T. Machado. “New Trends in Nanotechnology and Fractional Calculus Applications”. Springer (2010).
• ⁸
  D. Baleanu , J.T. Machado & A. Luo. “Fractional Dynamics and Control”. Springer (2012).
• ⁹
  D. Baleanu , G. Wu & S. Zeng. Chaos Analysis and Asymptotic Stability of Generalized Caputo Fractional Differential Equations. Chaos Solitons & Fractals, 102 (2017), 99-105.
• ¹⁰
  A. Boukhouima, K. Hattaf & N.Yousfi. Dynamics of a Fractional Order HIV Infection Model with Specific Functional Response and Cure Rate. International Journal of Differential Equations, 2017 (2017), 1-8.
• ¹¹
  J. Cannon & P.D. Chateau. Structural identification of an unknown source term in a heat equation. Inverse Prob., 14 (1998), 535-551.
• ¹²
  M. Caputo & M. Fabrizio. A new Definition of Fractional Derivative without Singular Kernel. Progr. Fract. Differ. Appl, 1 (2015), 73-85.
• ¹³
  M. Caputo & M. Fabrizio . Applications of New Time and Spatial Fractional Derivatives with Exponential Kernels. Progr. Fract. Differ. Appl , 2 (2016), 1-11.
• ¹⁴
  M. Caputo & M. Fabrizio . On the notion of fractional derivative and applications to the hysteresis phenomena. Meccanica, 52 (2017), 3043-3052.
• ¹⁵
  F. Dib & M. Kirane. An Inverse Source Problem For a Two Terms Time-fractional Diffusion Equation. Bol. Soc. Paran. Mat., 2022(40) (2022), 1-15.
• ¹⁶
  J. Gomez-Aguilar, M. Lopez-Lopez, V.A. Martinez, D. Baleanu & H. Khan. Chaos in a Cancer Model via Fractional Derivatives with Exponential Decay and Mittag Leffler Law. Entropy, 19(12) (2017), 1-19.
• ¹⁷
  E.D. Goufo. Application of the Caputo-Fabrizio Fractional Derivative without Singular Kernel to Korteweg-de Vries-Burgers Equation. Mathematical Modelling and Analysis, 21(2) (2016), 188-198.
• ¹⁸
  A. Hasanov. Simultaneous determination of source terms in a linear parabolic problem from the final overdetermination: weak solution approach. J. Math. Anal. App., 330 (2007), 766-779.
• ¹⁹
  A. Hasanov. An inverse source problem with single Dirichlet type measured output data for a linear parabolic equation. Appl. Math. Lett., 24 (2011), 1269-1273.
• ²⁰
  A. Hasanov, M. Otelbaev & B. Akpayev. Inverse heat conduction problems with boundary and final time measured output data. Inverse Prob. Sci. Eng. , 19(7) (2011), 985-1006.
• ²¹
  A. Hasanov & B. Pektas. Identification of an unknown time-dependent heat source term from over- specified Dirichlet boundary data by conjugate gradient method. Comput. Math. Appl., 65 (2013), 42-570.
• ²²
  B. Jin & W. Rundell. A tutorial on inverse problems for anomalous diffusion processes. Inverse Problems , 31 (2015), 1-40.
• ²³
  T. Johansson & D. Lesnic. Determination of a spacewise dependent heat source. J. Comput. Appl. Math., 209 (2007), 66-80.
• ²⁴
  A.R. Kanth & N. Garg. Computational Simulations for Solving a Class of Fractional Models via Caputo-Fabrizio Fractional Derivative. Procedia Computer Science, 125 (2018), 476-482.
• ²⁵
  M. Kirane & S. Malik. Determination of an unknown source term and the temperature distribution for the linear heat equation involving fractional derivative in time. Applied Mathematics and Computation, 218(1) (2011), 163-170.
• ²⁶
  M. Kirane, S. Malik & M. Al-Gwaizb. An inverse source problem for a two dimensional time fractional difusión equation with nonlocal boundary conditions. Math. Meth. Appl. Sci., 36(9) (2012), 1056-1069.
• ²⁷
  M. Kulbay, B. Mukanova & C. Sebu. Identification of separable sources for advection-diffusion equations with variable diffusion coefficient from boundary measured data. Inverse Prob. Sci. Eng. , 25(2) (2017), 279-308.
• ²⁸
  F. Mainardi & P. Paradisi. Fractional Diffusive Waves. Journal of Computational Acoustics, 9(4) (2001), 1417-1436.
• ²⁹
  M. Kirane, M. Sadybekov & A. Sarsenbi. On an inverse problem of reconstructing a subdiffusion process from nonlocal data. Math. Meth. Appl. Sci. , 42(6) (2019), 2043-2052.
• ³⁰
  O.J. Peter, A.S. Shaikh, M.O. Ibrahim, K.S. Nisar, D. Baleanu , I. Khan & A.I. Abioye. Analysis and Dynamics of Fractional Order Mathematical Model of Covid-19 in Nigeria Using Atangana-Baleanu Operator. Computers, Materials and Continua, 66(2) (2021), 1823-1848.
• ³¹
  A. Prilevko & D. Tkachenko. Well-posedness of the inverse source problem for parabolic systems. Differ. Equ., 40 (2004), 1619-1626.
• ³²
  N.A. Salti, E. Karimov & S. Kerbal. Boundary-value problems for fractional heat equation involving Caputo-Fabrizio derivative. NTMSCI 4, 4 (2016), 79-89.
• ³³
  A. Shaikh, A. Tassaddiq, K.S. Nisar & D. Baleanu . Analysis of differential equations involving Caputo-Fabrizio fractional operator and its applications to reaction-diffusion equations. Advances in Difference Equations, 178 (2019), 188-198.
• ³⁴
  G.S. Teodoro, J.T. Machado & E.C. de Oliveira. A review of definitions of fractional derivatives and other operators. Journal of Computational Physics, 388 (2019), 195-208.
• ³⁵
  M.I. Troparevsky, S. Seminara & M. Fabio. “Chapter 1: A Review on Fractional Differential Equations and a Numerical Method to Solve Linear Boundary Value Problems”. In: W. Legnani and T. Moschandreou Eds., Nonlinear Systems Volume 2, IntechOpen Book Series. IntechOpen (2019).
• ³⁶
  X. Xiong, Q. Zhou & Y. Honb. On an inverse problem of reconstructing a subdiffusion process from nonlocal data. J. Math. Anal. Appl. , 393 (2012), 185-199.
• ³⁷
  M. Yamamoto. Conditional stability in determination of force terms of heat equations in a rectangle. Math. Comput. Model, 18 (1993), 79-88.
• ³⁸
  Y. Yu, P. Perdikaris & G. Karniadakis. Fractional modeling of viscoelasticity in 3D cerebral arteries and aneurysms. Journal of Computational Physics , 323 (2016), 219-242.
• ³⁹
  L. Zada, R. Nawaz, S. Ahsan, K.S. Nisar & D. Baleanu . New iterative approach for the solutions of fractional order inhomogeneous partial differential equations. AIMS Mathematics, 6(2) (DOI: 10.3934/math.2021084), 1348-1365.
  » https://doi.org/10.3934/math.2021084
• ⁴⁰
  P. Zhuang, F. Liu, V. Anh & I. Turner. Numerical methods for the variable order fractional advection diffusion equation with a nonlinear source term. SIAM J. Numer. Anal., 47(3) (2009), 1760-1781.

Publication Dates

• Publication in this collection
  05 Sept 2022
• Date of issue
  Jul-Sep 2022

History

• Received
  31 Aug 2021
• Accepted
  18 Feb 2022