Ulam stabilities of nonlinear coupled system of fractional differential equations including generalized Caputo fractional derivative

1.   Introduction

The fractional diffusion equation (FDE), which is obtained by replacing the first-order time derivative and/or second-order space derivative in the standard diffusion equation by a generalized derivative of fractional order respectively, were successfully used for modelling relevant physical processes, see ^{[1,5,14,15,21]}. Recently the research on inverse problems connected with fractional derivatives becomes more and more popular. Since Cheng in ^[2] studied an inverse problem on fractional diffusion equation, many topics are well discussed ^{[11,12,23,26,27]}. For the problem of fractional numerical differentiation, in ^[17,18], the authors give different regularization methods. In ^[13,26,30], some inverse source problems for fractional diffusion equations are considered. In ^[25,28], Liu et al consider the backward time-fractional diffusion problem. In ^[4,8,10,16], many results on inverse coefficient problems are established. For more reference on inverse problems for fractional diffusion equations, please consult the survey paper ^[9]. However, in some situation of anomalous diffusion, the diffusion indexes and the diffusion coefficient are unknown. This leads to determination of coefficients which is a classical inverse problem. It should be mentioned that most of the existing literature investigate the determination of only one unknown parameter or functions. However, in many practical situations, one wishes to simultaneously reconstruct more than one physical parameters. To the authors' knowledge, there are few works on this aspect. For examples, in ^[4], the fractional order $\alpha$ in ${}_{t}D_{*}^{\alpha}u(x, t) = \triangle u(x, t)$ is determined by an analytic method. In ^[24], the fractional orders $\alpha, \beta$ in ${}_{t}D_{*}^{\alpha}u(x, t) = -(-\triangle)^{\frac{\beta}{2}}u(x, t)$ are reconstructed by the classical Levenberg-Marquardt method based on disrete least squares functional. The simultaneous inversion for the fractional order $\alpha$ and the space-dependent diffusion coefficient has been considered in ^[2,10]. In this paper, we reconstruct three important parameters in the time-space fractional diffusion equation from only one boundary measurement.

Let us consider the time and space-symmetric fractional diffusion equation in one-dimensional space

subject to homogeneous Neumann boundary conditions

and the initial condition

where $u$ is a solute concentration, $\kappa > 0$ represents the diffusion coefficient. ${}_{t}D_{*}^\alpha$ is the Caputo time fractional derivative of order $\alpha$ ($0 < \alpha < 1$) with the starting point at $t = 0$ defined as follows ^[19]:

The symmetric-space fractional derivative $(-\triangle)^{\frac{\beta}{2}}$ of order $\beta$ ($1 < \beta\leq 2$) is defined by ^[3,6,7]. For readability, we reproduce the following definition for $(-\triangle)^{\frac{\beta}{2}}$, $1 < \beta\leq 2$:

Definition 1. ^[29] Suppose the Laplace operator $-\triangle$ has a complete set of orthonormal eigenfunctions $\varphi_n$ corresponding to eigenvalues $\lambda_n^2$ on a bounded domain $D$, i.e., $(-\triangle)\varphi_n = \lambda_n^2 \varphi_n$ on a bounded domain $D$, ${\bf{B}}(\varphi) = 0$ on $\partial D$ is one of the standard three homogeneous boundary conditions. Let

then for any $f\in \mathbb{G}_\gamma$: $(-\triangle)^{\frac{\beta}{2}}$ is defined by

In the case of $\alpha = 1, \, \beta = 2$, Eq (1.1) reduces to the classical diffusion equation. For $0 < \alpha < 1, \, \beta = 2$, Eq (1.1) models subdiffusion due to particles having long-tailed resting times. For $\alpha = 1, \, 1 < \beta < 2$, Eq (1.1) corresponds to the Lévy process. Hence the solution of (1.1) is important for describing the competition between these two anomalous diffusion processes.

If $\alpha, \beta, \kappa, f(x)$ are given, the solution for the direct problem of Eqs (1.1)–(1.3) can be obtained analytically by the method of separation of variables: Setting $u(x, t) = X(x)T(t)$ and substituting into (1.1) yields

Letting $\omega$ is the separation constant, we obtain two fractional ordinary linear differential equations for $X(x)$ and $T(t)$ as

respectively. Following the definition of fractional Laplacian $(-\triangle)^{\frac{\beta}{2}}$ defined on a bounded domain, Eq (1.5) can be expressed as

Hence under homogeneous Neumann conditions, the eigenvalues are $\omega_n = \lambda_n^\beta = (\frac{n\pi}{L})^\beta$($1 < \beta\leq 2$) for $n = 0, 1, 2, \cdots$ and the corresponding eigenfunctions are $x_n = \cos(\frac{n\pi x}{L})$. Finally, the analytical solution of Eqs (1.1)–(1.3) is

where

Here we have used the result $E_{\alpha, 1}(0) = 1$, where $E_{\alpha, \beta}(z)$ is the Mittag-Leffler function defined by

By a similar method to ^[22], we can prove that (1.8) certainly gives the weak solution to Eqs (1.1)–(1.3).

Consider the following inverse problem:

Given $g(t): = u(0, t)$ (or $\tilde{g}(t): = u(L, t)$) with unknown $\alpha, \beta, \, \kappa$, we want to recover the orders $\alpha, \beta$ and the coefficient $\kappa$ from the data $g(t)$ (or $\tilde{g}(t)$).

Usually $g(t)$ is measured and only available data on $g(t)$ is its perturbation $g^\delta(t)$, we assume that there exists a known noise level $\delta$ such that

where the norm $\|\cdot\|$ denotes $L^2$-norm.

In this paper, our main work is to give the uniqueness result on determination of $\alpha, \beta, \, \kappa$ from the data $g(t)$ and two numerical methods for solving the inverse problems. Although in the paper ^[24], the authors give a uniqueness result on a similar problem, the result holds only for $0 < \alpha < 1/2$ and sufficiently large $T$. Our result do not require this restriction and holds for $0 < \alpha < 1$ and a finite $T$. This is done by adding some more smoothness assumption on the initial data $f(x)$.

Throughout this paper, sometimes we denote the solution of the problem as $u(x, t) = u(\alpha, \beta, \kappa, x, t)$ to show its dependence on $\alpha, \beta, \kappa$.

2.   Inverse problem and its uniqueness

Now from (1.8), we have the relationship:

The uniqueness of the inverse problem is stated as follows:

Theorem 1. Suppose that $u_1(\alpha_1, \beta_1, \kappa_1; x, t)$ and $u_2(\alpha_2, \beta_2, \kappa_2; x, t)$ represent the solutions of the inverse problem with $\alpha = \alpha_1, \, \beta = \beta_1, \, \kappa = \kappa_1$ and $\alpha = \alpha_2, \, \beta = \beta_2, \, \kappa = \kappa_2$ respectively. We assume that the initial data satisfies

If $u_1(\alpha_1, \beta_1, \kappa_1; 0, t) = u_2(\alpha_2, \beta_2, \kappa_2; 0, t)$ $(0 < t < T)$, then

Proof. From Eq (2.1), we have

By assumption $u_1(\alpha_1, \beta_1, \kappa_1;0, t) = u_2(\alpha_2, \beta_2, \kappa_2;0, t)$, it yields that

for $0 < t < T$.

As an initial step, we will prove that $\alpha_1 = \alpha_2$.

By the definition of Mittag-Leffler function, we have

And by the same argument, we obtain

From Eq (2.4), the following result

holds true. Now, we need to estimate the terms:

and

According to the inequality

where $\delta_0, C_0$ are some positive constants, we get

Furthermore, we recall that $f(x)\in H^4(0, L)$, then

where $C_1$ is a constant independent on $n$.

For the term on the right-side hand of Eq (2.9),

holds, when $0 < \epsilon_0 < 3/2$ and $0 < \varepsilon < \alpha_1$. Since

then

On the other hand, by Cauchy-Schwarz inequality

Noting that $\epsilon_0 < 3/2$ and $1 < \beta_1 < 2$, then $2(1+\epsilon_0-4) < 2(\beta_1+\epsilon_0-4) < -1$. Combining (2.10), we can conclude that

Therefore,

with small $\varepsilon > 0$. Similarly, we can verify that

with small $\varepsilon > 0$.

Inserting the above two equalities into (2.7), we have

If $\alpha_2 < \alpha_1$, then

holds for $t\to 0.$ That means

for $t\to 0.$ However, this is impossible. Because if $\sum_{n = 1}^{\infty}\kappa_2 (\frac{n\pi}{L})^{\beta_2}f_n = 0$, then due to $f_n > 0, 1 < \beta_2\leq 2$,

This is a contradiction. Hence $\alpha_2 < \alpha_1$ does not hold. Similarly, $\alpha_2 > \alpha_1$ does not hold, too. As a result, we get $\alpha_1 = \alpha_2$.

Next, we will prove that $\beta_1 = \beta_2, \, \kappa_1 = \kappa_2$.

Now we have the following equation from (2.4)

for $\alpha_1 = \alpha_2 = \alpha$. By the analytic proposition of $E_{\alpha, 1}(t)$, the above equation holds for $t > 0$. First for $0 < \alpha < 1$, we have the Laplace transform of Mittag-Leffler function:

Taking Laplace transform on both sides of Eq (2.20), we have

i.e.,

That is

We can analytically continue both sides of (2.24) in $\eta$. So (2.24) holds for

Now we deduce that $\kappa_1(\frac{n\pi}{L})^{\beta_1} = \kappa_2(\frac{n\pi}{L})^{\beta_2}$. From (2.24), W. L. O. G., we assume that $\kappa_1(\frac{n\pi}{L})^{\beta_1} < \kappa_2(\frac{n\pi}{L})^{\beta_2}$. Then we can take a suitable disc $R_1$ which includes $-\kappa_1(\frac{n\pi}{L})^{\beta_1}|_{n = 1}$ but does not include $\{-\kappa_1(\frac{n\pi}{L})^{\beta_1}\}_{n\geq 2}\cup\{-\kappa_2(\frac{n\pi}{L})^{\beta_2}\}_{n\ge 1}$. According to (2.24) and the analyticity of both sides of (2.24), we have

and hence, the Cauchy integral formula and Cauchy integral theorem yield

However, $f_1\neq 0$, i.e., $2\pi if_1\neq 0$. Therefore $\kappa_1(\frac{n\pi}{L})^{\beta_1} < \kappa_2(\frac{n\pi}{L})^{\beta_2}$ does not hold.

By the same argument, $\kappa_1(\frac{n\pi}{L})^{\beta_1} > \kappa_2(\frac{n\pi}{L})^{\beta_2}$ does not hold, either. Therefore, there holds $\kappa_1(\frac{n\pi}{L})^{\beta_1} = \kappa_2(\frac{n\pi}{L})^{\beta_2}$. Now, we are in the position to prove $\kappa_1 = \kappa_2$, $\beta_1 = \beta_2$ from $\kappa_1(\frac{n\pi}{L})^{\beta_1} = \kappa_2(\frac{n\pi}{L})^{\beta_2}$. It is easy to see that if $\kappa_{1}\neq\kappa_{2}$ or $\beta_{1}\neq\beta_{2}$ then there must exist at least a constant $n$ such that $\kappa_1(\frac{n\pi}{L})^{\beta_1}\neq\kappa_2(\frac{n\pi}{L})^{\beta_2}$. By analysis of its contrapositivity, we easily have $\kappa_1 = \kappa_2$, $\beta_1 = \beta_2$.

In general, the conditions (2.2) is not easy for one to verify. Therefore we give the weak conditions for the uniqueness:

Theorem 2. Suppose that $u_1(\alpha_1, \beta_1, \kappa_1; x, t)$ and $u_2(\alpha_2, \beta_2, \kappa_2; x, t)$ represent the solutions of the inverse problem with $\alpha = \alpha_1, \, \beta = \beta_1, \, \kappa = \kappa_1$ and $\alpha = \alpha_2, \, \beta = \beta_2, \, \kappa = \kappa_2$ respectively. We assume that the initial data satisfies

If $u_1(\alpha_1, \beta_1, \kappa_1; 0, t) = u_2(\alpha_2, \beta_2, \kappa_2; 0, t)$ for $0 < t < T$, then

$\alpha_1 = \alpha_2, \, \beta_1 = \beta_2, \, \kappa_1 = \kappa_2$.

Proof. W. L. O. G, let $L = 1, \varphi_{n}(x) = \sqrt{2}\cos n\pi x$ for $n\in \bf{N}$ and $\varphi_{0}(x) = 1$. Then we consider

Then

Setting $f_{n} = (f, \varphi_n)\sqrt{2}$ for $n\in \bf{N}$ and $f_{0} = (f, 1)$, we have

Assume

Then due to Podlubny ^[19] (p.34),

Since $f\in D((-\triangle)^{{\beta_k}/{2}})$ and $(-\triangle)^{{\beta_k}/{2}}$ is self-adjoint for $k = 1, 2$, we have

Hence we have

Now we can prove that $\alpha_1 = \alpha_2$.

Indeed we assume that $\alpha_1 < \alpha_2$. Dividing (2.29) by $t^{\alpha_1}$, we obtain

By $\alpha_1 < \alpha_2$, letting $t \rightarrow 0$, we obtain $\frac{k_1}{\Gamma(\alpha_1+1)}(-\triangle)^{{\beta_1}/{2}}f(0) = 0$, which contradicts (2.29).

The uniqueness proof on $\beta$ and $\kappa$ is the same as that in Theorem 1.

Remark on (2.29). We can satisfy (2.29) by a generous condition by the comparison principle of $(-\triangle)^{{\beta}/{2}}$.

3.   Two numerical methods

In this section, we propose two numerical methods for solving this problem based least squares functional. The first is based on Tikhonov method in the function space. The second method is based on the classical Levenberg-Marquardt optimization method in the discrete Euclid space.

3.1.   Tikhonov method combined with the gradient flow

Denote $a = \left(\begin{array}{c} \alpha \\ \beta \\ \kappa \end{array}\right)\in \mathbb{R}^3$, let $u(x, t;a): = u(\alpha, \beta, \kappa)(x, t)$ be the unique solution of forward problem. A feasible way to numerical computation for the unknown $a$ is to solve the following minimization problem.

The gradient $\nabla J(a)$ of the functional $J(a)$ is given by

If let $\nabla J(a)$ = 0, then we can get the Euler equation for the minimizer. However, it is a nonlinear equation and is not easily be solved directly. Here we turn to the approximate solution by the iterative method. Using the gradient flow method with an initial value $a_0$, we get

where $t$ is the artificial time. Using a simple method, i.e. the explicit Euler method, we arrives the following iteration schemes with a time step size $\tau$:

i.e.,

3.2.   Levenberg-Marquardt method

Because in most of the practical applications, the data are measured at discrete times. Assume the measured data is given by $g^\delta(t_i), \, i = 0, 1, \cdots, q$. Let us consider the minimization problem in discrete case:

where $u(a; 0, t_i)$ is the computed data from the forward problem with a given $a$, which is used to fit the measured data. A standard method for solving this least squares problem is the Levenberg-Marquardt method with a damped parameter $\tilde{\lambda}$ which plays the same role as the regularization parameter $\lambda$ in Tikhonov method. For readability, we give the details of this algorithm:

A updated sequences is given by

where $\Delta a^j$ is the updated stepsize of $a^j$ in each iteration step $j$. We consider the the minimization problem about $\Delta a^j$ at each iteration step $j$:

Make the Taylor expansion for $u(a^j+\Delta a^j; 0, t_i)$ at $a^j$ and take a linear approximation, we have

Plus this into (3.6), we get

However, this least square problem is ill-posed due to the original problem, therefore we consider the Tikhonov method:

where $\nabla_a^{tr}u(a; 0, t_i)\cdot a^j$ is computed by finite difference method and is given by $\nabla_a^{tr}u(a; 0, t_i)\cdot a^j\approx \sum_{k = 1}^3 \frac{u(a_k^j+h; 0, t_i)-u(a_k^j; 0, t_i)}{h}\Delta a_k^j$. Now the minimization problem (3.8) is a linear problem and can be easily solved for the updated stepsize $\Delta a^j$ with a regularization $\tilde{\lambda}$.

3.3.   Numerical test

In this section, we consider a simple example to show the effectiveness of the aforementioned two algorithms, i.e., Tikhonov method and Levenberg-Marquardt method. We want to determine the parameters $(\alpha, \beta, \kappa)$ in the following problem

The exact solution is given by

Now the input data $g(t): = u(0, t)$ is obtained and the noisy data $g^\delta(t)$ is generated in the following way:

where $\sigma(t) = \theta r(t)$ and $r(t)$ is a random number between $[0, 1]$ and $\theta = \max\{g(t)\}*\eta\%$ is the noise level. In the numerical experiment, we fix the parameters $T = 1, \eta = 1$. The algorithm of calculating Mittag effler function is given in ^[20]. First we plot the solution $u(x, t = 1)$ of direct problem when $\alpha = 0.5, \beta = 1.2, \kappa = 0.10$, which is shown in Figure 1. The input exact data $g(t)$ is displayed in Figure 2.

The exact fractional orders and diffusion coefficient are ${ {{ α}} = {\bf{0.5}}, {{ β}} = {\bf{1.2}}, {{ κ}} ={\bf{ 0.10}}}$.

In the numerical test for the Tikhonov method, the parameters and their values in the computation are listed below:

1. Summing idex of $u(x, t)$ in (3.12): $n = 20$ is taken.

2. The number of $t_i = i/q\in[0, 1], (i = 0, \cdots, q)$ and $q = 20$ in two methods ($q$ is the total points for trapezoidal rule of the numerical integral in (3.4)).

3. Stepsizes $\tau_{s} (s = \alpha, \beta, \kappa)$ in (3.4) for coputing the derivatives $\frac{\partial u}{\partial \alpha}$, $\frac{\partial u}{\partial \beta}$ and $\frac{\partial u}{\partial \kappa}$ by finite difference method: $\tau_\alpha = \tau_\beta = \tau_\kappa = 0.01$.

4. The initial guess for $(\alpha, \beta, \kappa) = (0.4, 1.0, 0.05)$.

5. The regularization parameter $\lambda = 0.001$.

6. The iterative step size is $\tau = (0.07, 0.01, 0.02)$.

6. Stop criterion: when $\|a^{j+1}-a^j\|\leq 0.01$.

Finally we obtain the approximate value $(\alpha, \beta, \kappa) = (0.5289, 1.1836, 0.0943)$. Keep the same parameters with the damped parameter $\tilde{\lambda} = 0.001$, we use the Levenberg-Marquardt method to get the approximate value where we don't use the Matlab optimization toolbox on Levenberg-Marquardt method.

This result shows the numerical methods are effective. Here we list more results using the above parameters for the Levenberg-Marquardt method. First we fix the initial guess for $\alpha = 0.1, \kappa = 0.05$ and let $\beta$ range from 1 to 1.4. The numerical results are displayed in Table 1.

The numerical results show that the methods are stable.

4.   Conclusions

In this paper, we give the proof of uniqueness for determining three parameters in a time-space fractional diffusion equation by means of observation data from accessible boundary. By our uniqueness result a Tikhonov method and the Levenberg-Marquardt method are tested preliminarily. Some further research on the stability of the proposed methods will be investigated in the future.

Acknowledgments

The authors are grateful to the anonymous referees for their insightful comments and suggestions which contributed to greatly improve the original version of the manuscript. The authors would like to thank for the Prof. Masahiro Yamamoto's important comments on Theorem 2.

This work is partially supported by the Natural Science Foundation of China (No. 11661072), the Natural Science Foundation of Northwest Normal University, China (No. NWNU-LKQN-17-5).

Conflict of interest

There is no conflict of interest.