Parameter Identification in the Two-Dimensional Riesz Space Fractional Diffusion Equation

Abstract

This paper presents the application of the swarm intelligence algorithm for solving the inverse problem concerning the parameter identification. The paper examines the two-dimensional Riesz space fractional diffusion equation. Based on the values of the function (for the fixed points of the domain) which is the solution of the described differential equation, the order of the Riesz derivative and the diffusion coefficient are identified. The paper includes numerical examples illustrating the algorithm’s accuracy.

1. Introduction

Models with fractional derivatives have found applications in many fields of science and engineering, such as control theory [1], mechanics [2], image processing [3] and heat conduction. The use of the fractional derivatives in modeling the heat conduction phenomena has been presented in the papers [4,5], wherein the models with the fractional derivative of Riemann–Liouville and Caputo type was considered. In case of the heat conduction problems the equations with the fractional derivatives are particularly useful for modeling the phenomena that take place in the porous materials and composites due to the fact that the anomalous diffusion process occurs there. For such materials the models with the fractional derivatives give better results than the models based on the classic derivatives, as has been shown, for example, in papers [6,7]. In the articles [7,8] the heat conduction process was considered in composites and porous aluminum. More about the applications of the fractional derivatives can be found in [9,10,11,12].

In this paper the model described by the differential equation with the fractional derivative of the Riesz type [13] is being considered. The definition of this derivative is based on the left and right-hand sided Riemann–Liouville derivatives. The algorithm for solving the inverse problem for model with the partial Riesz derivative with respect to the spacial variable is being presented. The inverse problems in general are a class of problems related to the proper selection of the input parameters to obtain the expected values as an output. To solve the inverse problem unambiguously, additional information is required. In this article the problem of reconstructing the diffusion coefficient and the order of the fractional derivative is being solved and the additional information is a set of values of a function which is a solution of the considered differential equation for the chosen points of the domain.

Solving the inverse problem often requires solving the direct problem multiple times for the chosen values of the reconstructed parameters; in our case these are the diffusion coefficient and the order of the fractional derivative. For solving the direct problem, the alternating direction implicit method [13], which allows one to solve the multidimensional problem by iteratively solving a number of one-dimensional problems, was used. The functional representing the error of the approximate result was constructed using the obtained solution of the direct problem and then it was minimized by application of the real ant colony optimization [14] algorithm. That is the probabilistic artificial intelligence algorithm inspired by the behavior of the swarm of ants. In this paper, a numerical example presenting the accuracy of the method for various numbers of measurement points and various input data disturbance is also presented.

2. Formulation of the Problem

We consider a two-dimensional differential equation with the fractional derivative of a Riesz type with respect to the spacial variable:

$$c\varrho \frac{\partial u(x,y,t)}{\partial t}={\lambda }_x\frac{{\partial }^{{\alpha }_1}u(x,y,t)}{{\partial \left|x\right|}^{{\alpha }_1}}+{\lambda }_y\frac{{\partial }^{{\alpha }_2}u(x,y,t)}{{\partial \left|y\right|}^{{\alpha }_2}}+d(u,x,y,t),$$

(1)

defined on the domain $D=\{(x,y,t):x\in [0,L_x],y\in [0,L_y],t\in [0,T],L_x,L_y,T\in {\mathbb{R}}_{+}\}$. Such an equation is called a two-dimensional Riesz space fractional diffusion equation [13], so the ${\lambda }_x$ and ${\lambda }_y$ coefficients are called diffusion coefficients from now on. In this equation the fractional derivatives of order ${\alpha }_1,{\alpha }_2\in (1,2)$ are the derivatives of a Riesz type defined as follows:

$$\begin{array}{c}\hfill \frac{{\partial }^{{\alpha }_1}u(x,y,t)}{{\partial \left|x\right|}^{{\alpha }_1}}=-c_{{\alpha }_1}\left(\frac{{\partial }^{{\alpha }_1}u(x,y,t)}{\partial x^{{\alpha }_1}}+\frac{{\partial }^{{\alpha }_1}u(x,y,t)}{\partial {(-x)}^{{\alpha }_1}}\right),\\ \hfill \frac{{\partial }^{{\alpha }_2}u(x,y,t)}{{\partial \left|y\right|}^{{\alpha }_2}}=-c_{{\alpha }_2}\left(\frac{{\partial }^{{\alpha }_2}u(x,y,t)}{\partial y^{{\alpha }_2}}+\frac{{\partial }^{{\alpha }_2}u(x,y,t)}{\partial {(-y)}^{{\alpha }_2}}\right),\end{array}$$

(2)

where $c_{\alpha }=\frac{1}{2cos\left(\frac{\pi \alpha }{2}\right)}$ and:

$$\begin{array}{c}\hfill \frac{{\partial }^{{\alpha }_1}u(x,y,t)}{\partial x^{{\alpha }_1}}=\frac{1}{\mathsf{\Gamma }(2-{\alpha }_1)}\frac{{\partial }^2}{\partial x^2}\underset{0}{\overset{x}{\int }}\frac{u(s,y,t)}{{(x-s)}^{{\alpha }_1-1}}ds,\\ \hfill \frac{{\partial }^{{\alpha }_1}u(x,y,t)}{\partial {(-x)}^{{\alpha }_1}}=\frac{{(-1)}^2}{\mathsf{\Gamma }(2-{\alpha }_1)}\frac{{\partial }^2}{\partial x^2}\underset{x}{\overset{L_x}{\int }}\frac{u(s,y,t)}{{(s-x)}^{{\alpha }_1-1}}ds,\\ \hfill \frac{{\partial }^{{\alpha }_2}u(x,y,t)}{\partial y^{{\alpha }_2}}=\frac{1}{\mathsf{\Gamma }(2-{\alpha }_2)}\frac{{\partial }^2}{\partial y^2}\underset{0}{\overset{y}{\int }}\frac{u(x,s,t)}{{(y-s)}^{{\alpha }_2-1}}ds,\\ \hfill \frac{{\partial }^{{\alpha }_2}u(x,y,t)}{\partial {(-y)}^{{\alpha }_2}}=\frac{{(-1)}^2}{\mathsf{\Gamma }(2-{\alpha }_2)}\frac{{\partial }^2}{\partial y^2}\underset{y}{\overset{L_y}{\int }}\frac{u(x,s,t)}{{(s-y)}^{{\alpha }_2-1}}ds,\end{array}$$

(3)

where $\mathsf{\Gamma }$ is the Gamma function. To the Equation (1) we add the initial condition:

$$u(x,y,0)=f(x,y),x\in [0,L_x],y\in [0,L_y],$$

(4)

and homogeneous Dirichlet boundary conditions:

$$\begin{array}{c}\hfill u(0,y,t)=u(L_x,y,t)=0,\\ \hfill u(x,0,t)=u(x,L_y,t)=0.\end{array}$$

(5)

The differential equations of the fractional order have found applications in modeling various kinds of phenomena, such as the diffusion or heat conduction in porous materials [4,7]. In the model presented above we assume that ${\alpha }_1={\alpha }_2=\alpha$ and ${\lambda }_x={\lambda }_y=\lambda$. The considered inverse problem is to reconstruct a value of an order of the derivative $\alpha$ and the diffusion coefficient $\lambda$ using the values of the u function for the fixed points of the domain. We call these values the input data for the inverse problem and denote them as:

$${\widehat{U}}_{i,j}^k=u(x_i,y_j,t_k),i,j=1,2,\dots ,N_1,k=1,2,\dots ,N_2,$$

where $N_1$ is the number of measurement points and $N_2$ is the number of measurements taken for each point. By solving the differential Equation (1) (numerically) for the fixed values of the $\alpha$ and $\lambda$ parameters we can get the values of the function $U_{i,j}^k(\alpha ,\lambda )$ at the measurement points. By taking these values, in each point of the mesh, and comparing them with the input data, we can build a function that defines the error of the approximate solution:

$$J(\alpha ,\lambda )=\sum _{i,j=1}^{N_1}\sum _{k=1}^{N_2}{\left(U_{i,j}^k(\alpha ,\lambda )-{\widehat{U}}_{i,j}^k\right)}^2.$$

(6)

By minimizing the above function, we get the approximate values of the $\alpha$ and $\lambda$ parameters.

3. Methods of Solution

The algorithm for solving the considered inverse problem consists of two major parts: solving the direct problem—determining the solution of the differential Equation (1) with the given conditions; and finding the minimum of the objective function (6).

3.1. Solution of the Direct Problem

In order to solve the differential Equation (1) in the considered area we introduce a grid:

$$S=\{(i\Delta x,j\Delta y,k\Delta t):\Delta x=\frac{L_x}{N},\Delta y=\frac{L_y}{M},\Delta t=\frac{T}{K},i=0,1,\dots N,j=0,1,\dots ,M,k=0,1,\dots ,K\}.$$

Let us take a following notation: $U_{i,j}^k\approx u(x_i,y_j,t_k)=u(i\Delta x,j\Delta y,k\Delta t)$. The initial condition may be denoted as $U_{i,j}^0=f(x_i,y_j)=f_{i,j}$, while for the source component d we use notation $d_{i,j}^k=d(u_{i,j}^k,x_i,y_j,t_k)$. We approximate the first derivative with respect to time in the following way:

$$\frac{\partial u(x_i,y_j,t_k)}{\partial t}\approx \frac{U_{i,j}^k-U_{i,j}^{k-1}}{\Delta t}.$$

The fractional derivative (3) with respect to the x variable is approximated by the shifted Grünwald-Letnikov formula [15,16]:

$$\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{\partial x^{\alpha }}\approx \frac{1}{{(\Delta x)}^{\alpha }}\sum _{r=0}^{i+1}{g}_{\alpha }^{\left(r\right)}{U}_{i-r+1,j}^k,$$

(7)

$$\frac{{\partial }^{\alpha }u(x_i,y_j,t_k)}{\partial {(-x)}^{\alpha }}\approx \frac{1}{{(\Delta x)}^{\alpha }}\sum _{r=0}^{N-i+1}{g}_{\alpha }^{\left(r\right)}{U}_{i+r-1,j}^k,$$

(8)

where:

$$g_{\alpha }^{\left(r\right)}={(-1)}^r\frac{\mathsf{\Gamma }(\alpha +1)}{\mathsf{\Gamma }(\alpha -r+1)\mathsf{\Gamma }(r+1)}.$$

We use similar method to approximate the derivative with respect to the y variable. By taking into account all of the above dependencies, Equation (1) may be written in the form:

$$\begin{array}{c}\hfill \frac{U_{i,j}^k-U_{i,j}^{k-1}}{\Delta t}=-\frac{\lambda c_{\alpha }}{{(\Delta x)}^{\alpha }}\left[\sum _{r=0}^{i+1}{g}_{\alpha }^{\left(r\right)}{U}_{i-r+1,j}^k+\sum _{r=0}^{N-i+1}{g}_{\alpha }^{\left(r\right)}{U}_{i+r-1,j}^k\right]\\ \hfill -\frac{\lambda c_{\alpha }}{{(\Delta y)}^{\alpha }}\left[\sum _{r=0}^{j+1}{g}_{\alpha }^{\left(r\right)}{U}_{i,j-r+1}^k+\sum _{r=0}^{M-j+1}{g}_{\alpha }^{\left(r\right)}{U}_{i,j+r-1}^k\right]+f_{i,j}^{k-1}.\end{array}$$

(9)

Next we introduce the operator ${\delta }_x$ defined as follows:

$${\delta }_x{u}_{i,j}^k=-\frac{\lambda c_{\alpha }}{{(\Delta x)}^{\alpha }}\left[\sum _{r=0}^{i+1}{g}_{\alpha }^{\left(r\right)}{U}_{i-r+1,j}^k+\sum _{r=0}^{N-i+1}{g}_{\alpha }^{\left(r\right)}{U}_{i+r-1,j}^k\right]$$

We also define the operator ${\delta }_y$ similarly. By substituting into (9), we get the following equation:

$$(1-\Delta t{\delta }_x-\Delta t{\delta }_y)U_{i,j}^k=U_{i,j}^{k-1}+\Delta t{f}_{i,j}^{k-1}.$$

(10)

To solve the direct problem we use the alternating direction implicit method (ADIM). This method reduces the multidimensional problem to a several independent one-dimensional problems to solve. We write the Equation (10) in a following form:

$$(1-\Delta t{\delta }_x)(1-\Delta t{\delta }_y)U_{i,j}^k=U_{i,j}^{k-1}+\Delta t{f}_{i,j}^{k-1}.$$

(11)

ADIM is a iterative method and in the considered case we may define it in two steps:

• First we solve the problem in the direction of the $Ox$ axis (for the fixed $y_j$):

  $$(1-\Delta t{\delta }_x)U_{i,j}^{*}=U_{i,j}^{k-1}+\Delta t{f}_{i,j}^{k-1}.$$

  (12)
  This way we get the intermediate solution $U_{i,j}^{*}$.
• In the next step we solve the problem in the direction of $Oy$ axis (for the fixed $x_i$):

  $$(1-\Delta t{\delta }_y)U_{i,j}^k=U_{i,j}^{*}.$$

  (13)

In the considered problem the Dirichlet boundary conditions are zeros. Due to that, we may easily determine the boundary values for the solution $U^{*}$ which are:

$$U_{0,j}^{*}=(1-\Delta t{\delta }_y)U_{0,j}^k=0,U_{N,j}^{*}=(1-\Delta t{\delta }_y)U_{N,j}^k=0.$$

We must write the differential schemes (12)–(13) in the matrix form, thereby obtaining the system of equations, solutions of which are the values of the function u at the grid points. More about the ADIM method, its stability and convergence may be found in the [13].

3.2. Minimalization of the Objective Function

The second necessary step to solving the considered inverse problem is the minimization of the functional (6). For this we use the real ant colony optimization (RealACO) algorithm [14], which is the swarm artificial intelligence algorithm. We describe it briefly below, and for that, use the following notation:

J—minimalized function (objective function), n—size of the problem,
nT—number of threads, M = nT · p—number of ants in the population,
I—number of iterations, L—number of pheromone spots, q, ξ—algorithm parameters.

The pheromone spots mean the potential solutions. The ants modify these solutions while looking for better results. We now present the consecutive steps of the algorithm.

Start of the algorithm

1.

Set the algorithm’s input parameters: $L,M,I,nT,q,\xi$.

2.

Generate L initial solutions serving as pheromone spots. Assign them to the $T_0$ set (which is the initial archive).

3.

Compute the values of the objective function for all of the pheromone spots (parallel computing) and sort the archive T0 from the best solution to the worst.

Iteration process

4.

Assign the probabilities to the pheromone spots according to the formula:

$$p_l=\frac{{\omega }_l}{{\displaystyle \sum _{l=1}^L}{\omega }_l}l=1,2,\dots ,L,$$

where weights ${\omega }_l$ are connected to the l-th solution and are expressed by the formula:

$${\omega }_l=\frac{1}{qL\sqrt{2\pi }}·e^{\frac{-{(l-1)}^2}{2q^2{L}^2}}.$$

5.

The ant must randomly chose the l-te solution with the probability $p_l$.

6.

The ant must transform the j-th coordinate ($j=1,2,\dots ,n$) of the l-th solution $s_j^l$ by sampling the neighborhood using the probability density function (in this case the Gauss function):

$$g(x,\mu ,\sigma )=\frac{1}{\sigma \sqrt{2\pi }}·e^{\frac{-{(x-\mu )}^2}{2{\sigma }^2}},$$

where $\mu =s_j^l,\sigma =\frac{\xi }{L-1}{\displaystyle \sum _{p=1}^L}|s_j^p-s_j^l|$.

7.

The steps 5 and 6 must be repeated by every ant. This way, we get M new solutions (pheromone spots).

8.

Partition the new solutions into $nT$ groups. Compute the value of the objective function for each of the new solutions (parallel computing).

9.

Add the new solutions to the archive $T_i$; sort them by the quality and reject the M worst solutions.

10.

The steps 3–9 must be repeated I times.

Knowing the values of the parameters $L,M$ and I, we can determine the number of the objective function calls during the algorithm execution. This number is $L+M·I$. More about the method and its application for the inverse problems (particularly in the field of heat conductivity) can be found in [5,17,18,19].

4. Numerical Example

We now present the method described in the previous section by using it in a numerical example. We consider the Equation (1) with the homogenous Dirichlet boundary conditions. We also take under consideration the following data:

$$c=900.0,\varrho =2700.0,L_x=L_y=1.0,T=300.0,f(x,y)=\mathrm{150,000}{(x-1)}^4{x}^2{(y-1)}^4{y}^2,$$

$$\begin{array}{c}d(u,x,y,t)=-\mathrm{1,215,000,000}{(x-1)}^4{x}^2{(y-1)}^4{y}^2\hfill \\ +\mathrm{191,492}(t-300)(x^{0.5}-8x^{1.5}+19.2x^{2.5}-18.2857x^{3.5}+6.09524x^{4.5}){(y-1)}^4{y}^2\\ +(t-300)(\mathrm{65,297.4}{(1-x)}^{4.5}-\mathrm{587,677}{(1-x)}^{3.5}(0.0833333+x)\\ +\mathrm{514,217}{(1-x)}^{2.5}(0.025+0.166667x+x^2)){(y-1)}^4{y}^2+\mathrm{191,492}(t-300){(x-1)}^4\\ x^2(1.y^{0.5}-8.y^{1.5}+19.2y^{2.5}-18.2857y^{3.5}+6.09524y^{4.5})+\mathrm{32,648.7}(-300+t){(x-1)}^4{x}^2\\ \hfill \left(2{(1-y)}^{4.5}-18{(1-y)}^{3.5}(0.0833333+y)+15.75{(1-y)}^{2.5}(0.025+0.166667y+y^2)\right).\end{array}$$

The exact solution of such a problem for $\alpha =1.5$ and $\lambda =240$ is a function:

$$u(x,y,t)=500(300-t)x^2{y}^2{(x-1)}^4{(y-1)}^4.$$

As was described in the Section 2 the reconstructed parameters of the considered model are order of the derivative $\alpha$ and the diffusion coefficient $\lambda$. The aim of this numerical example is to illustrate the method’s accuracy. Therefore, the input data for the inverse problem were generated for the exact values $\alpha =1.5$, $\lambda =240.0$ using the grid of size $150\times 150\times 600$. For the algorithm for solving the inverse problem, the grid of size $100\times 100\times 300$ was used. Using a different grid size for generating the input data for the inverse problem and a different size for solving the inverse problem was necessary to avoid the occurrence of the so-called inverse crimes [20,21]. Additionally, in order to examine the stability of the algorithm, the input data were disturbed by the pseudo-random normal distribution errors of size $1\%$ and $3\%$. Running time of the algorithm on a standard PC was about 14 h.

In this example we considered two cases of measurement point locations:

• Four points placed in the “corners” of the domain: $A(0.2,0.8),B(0.8,0.8),C(0.2,0.2),D(0.8,0.2)$; then we have $N_1=4$.
• Two points placed in the opposite “corners” of the domain: $B(0.8,0.8),C(0.2,0.2)$; then we have $N_1=2$.

The locations of the measurement points are presented in the Figure 1. For each configuration of points, we simulated the measurements every 1 s, 2 s and 4 s ($N_2=300,N_2=150,N_2=75$).

Table 1. Results of calculations in the case of four measurement points (${\overline{\alpha }}_i$—reconstructed value of order derivative, ${\overline{\lambda }}_i$—reconstructed value of diffusion coefficient, $\delta$—the relative error of reconstruction, $\sigma$—standard deviation, J—the value of functional).

Measurements Intervals	Noise	$\overline{\mathit{\alpha }}$	${\mathit{\delta }}_{\overline{\mathit{\alpha }}}[\%]$	${\mathit{\sigma }}_{\overline{\mathit{\alpha }}}$	$\overline{\mathit{\lambda }}$	${\mathit{\delta }}_{\overline{\mathit{\lambda }}}[\%]$	${\mathit{\sigma }}_{\overline{\mathit{\lambda }}}$	J
1 s
	$0\%$	$1.4944$	$0.38$	$0.0677$	$241.66$	$0.69$	$33.19$	$0.0135$
	$1\%$	$1.5114$	$0.76$	$0.0502$	$233.87$	$2.55$	$21.77$	$1.7948$
	$3\%$	$1.5083$	$0.56$	$0.0967$	$235.54$	$1.86$	$46.92$	$16.2992$
2 s
	$0\%$	$1.5181$	$1.21$	$0.0474$	$231.30$	$3.62$	$21.56$	$0.0243$
	$1\%$	$1.5129$	$0.86$	$0.0588$	$232.90$	$2.95$	$27.43$	$0.8427$
	$3\%$	$1.5094$	$0.63$	$0.0246$	$234.68$	$2.21$	$11.53$	$7.9230$
4 s
	$0\%$	$1.4937$	$0.41$	$0.0444$	$241.84$	$0.76$	$21.25$	$0.0021$
	$1\%$	$1.4959$	$0.27$	$0.0626$	$240.91$	$0.38$	$29.82$	$0.4088$
	$3\%$	$1.5480$	$3.21$	$0.0511$	$219.71$	$8.45$	$21.82$	$3.5613$

and

Table 2. Results of calculations in the case of two measurement points (${\overline{\alpha }}_i$—reconstructed value of order derivative, ${\overline{\lambda }}_i$—reconstructed value of diffusion coefficient, $\delta$—the relative error of reconstruction, $\sigma$—standard deviation, J—the value of functional).

Measurements Intervals	Noise	$\overline{\mathit{\alpha }}$	${\mathit{\delta }}_{\overline{\mathit{\alpha }}}[\%]$	${\mathit{\sigma }}_{\overline{\mathit{\alpha }}}$	$\overline{\mathit{\lambda }}$	${\mathit{\delta }}_{\overline{\mathit{\lambda }}}[\%]$	${\mathit{\sigma }}_{\overline{\mathit{\lambda }}}$	J
1 s
	$0\%$	$1.5115$	$0.77$	$0.1666$	$233.47$	$2.71$	$72.28$	$0.0974$
	$1\%$	$1.4901$	$0.65$	$0.0793$	$243.52$	$1.46$	$38.56$	$1.7712$
	$3\%$	$1.5084$	$0.56$	$0.1738$	$235.38$	$1.92$	$50.94$	$16.1719$
2 s
	$0\%$	$1.4492$	$3.38$	$0.1637$	$261.82$	$9.09$	$21.76$	$0.0191$
	$1\%$	$1.4922$	$0.51$	$0.0502$	$242.45$	$1.02$	$21.35$	$0.8085$
	$3\%$	$1.5066$	$0.44$	$0.0813$	$235.71$	$1.78$	$35.82$	$7.8638$
4 s
	$0\%$	$1.4773$	$1.51$	$0.0666$	$249.11$	$3.79$	$32.52$	$0.0007$
	$1\%$	$1.4664$	$2.23$	$0.1769$	$254.01$	$5.83$	$43.95$	$0.4032$
	$3\%$	$1.5796$	$5.31$	$0.0821$	$207.71$	$13.45$	$35.06$	$3.4968$

present the obtained results of the $\alpha$ and $\lambda$ parameter reconstructions for various sizes of disturbance of the input data and different numbers of measurements taken. The presented algorithm is probabilistic, and therefore for each case it was executed five times. The tables include the values of standard deviation calculated for these executions.

In the case of four measurement points, the best solutions (the smallest errors) were received for measurements taken every 1 s. The errors of the order of the derivative $\alpha$ were no more than $1\%$ and the biggest value of the error of reconstruction of the $\lambda$ coefficient was gotten for the input data disturbed by $1\%$ error. In each case the lower the value of the functional, the lower the disturbance of the input data. In the case of four measurement points, the worst solution we got was for the input data disturbed by the $3\%$ error and measurements taken every 4 s. That was the consequence of, first of all, the biggest error of the input data, and second of all, the smallest number of measurements taken ($N_2=75$). The values of the standard deviation of the obtained results are more or less at a similar level. The value of the standard deviation was affected by the RealACO algorithm parameters and number of the algorithm calls. During some of the executions (when the measurements were taken every 4 s and the disturbance was $3\%$) the solution ended falling into a local minimum, which significantly affected the value of the standard deviation. The proper selection of the RealACO parameters and increasing the number of algorithm calls caused a decrease of the value of the standard deviation. In case of the two measurement points, we can observe that in most cases (but not in all of them) the obtained errors of the solutions are bigger than in the case of the four measurement points. The values of the functional are, however, at a similar level in both cases. By comparing the two tables, we can also observe that the values of the standard deviations are bigger in most cases for the four measurement points. In this case, increasing the number of measurement points should positively affect the obtained results. For the obtained approximate values of the $\alpha$ and $\lambda$ parameters, the absolute errors of the function u reconstruction were calculated for the measurement points and are presented in the

Table 3. The errors of the function u reconstruction for the measurement points in the case of four measurement points (${\Delta }_{avg}$—average absolute error, ${\Delta }_{max}$—maximum average error).

Noise	$0\%$	$1\%$	$3\%$	$0\%$	$1\%$	$3\%$	$0\%$	$1\%$	$3\%$
	co 1 s			co 2 s			co 4 s
${\Delta }_{avg}$	$0.0021$	$0.0024$	$0.0026$	$0.0042$	$0.0038$	$0.0022$	$0.0021$	$0.0018$	$0.0141$
${\Delta }_{max}$	$0.0091$	$0.0269$	$0.0148$	$0.0223$	$0.0367$	$0.0258$	$0.0065$	$0.0074$	$0.0485$

and

Table 4. The errors of the function u reconstruction for the measurement points in the case of two measurement points (${\Delta }_{avg}$—average absolute value, ${\Delta }_{max}$—maximum absolute value).

Noise	$0\%$	$1\%$	$3\%$	$0\%$	$1\%$	$3\%$	$0\%$	$1\%$	$3\%$
	co 1 s			co 2 s			co 4 s
${\Delta }_{avg}$	$0.0071$	$0.0027$	$0.0037$	$0.0064$	$0.0018$	$0.0053$	$0.0019$	$0.0035$	$0.0245$
${\Delta }_{max}$	$0.0365$	$0.0073$	$0.0176$	$0.0163$	$0.0071$	$0.0299$	$0.0041$	$0.0087$	$0.0759$

. In each case the errors were very small. The average values of the relative error in case of four measurement points for the $3\%$ disturbance and measurements intervals 1 s, 2 s and 4 s were $0.31$%, $0.51$% and $6.03$% accordingly. In case of two measurement points the errors were $0.34$%, $1.77$% and $10.14$%, accordingly, for those measurements intervals. The higher values of the error were the consequence of the fact that the solution on the boundary of the considered area was equal to zero. The absolute errors in those cases were small and for the four measurement points did not exceed the values of $0.0141$ on average and $0.0485$ at maximum. For the two measurement points, in turn, these values do not exceed the $0.0245$ average and $0.0759$ maximum.

Figure 2 presents the absolute errors of the solution u reconstruction at the measurement points B and C in cases of two measurement points, the smallest number of measurements taken and various sizes of the input data disturbance. In case of the accurate input data and 1% disturbance, the errors of the solution reconstruction for the measurement points were minimal. The errors were slightly higher with 3% disturbance.

5. Conclusions

The aim of the paper was to reconstruct the order of the derivative and the diffusion coefficient while the values of the solution in the selected points of the domain were known. The obtained results show that for the proper number of the measurement points and a reasonable size of the input data disturbance, the presented algorithm allows one to reconstruct the parameters of interest well. Moreover, the absolute values of the considered solution reconstruction are small. Consecutive executions of the algorithm gave similar results, which was proven by the values of the standard deviation obtained.

In the future we plan to try to improve the algorithm by combining the heuristic algorithm used with the deterministic algorithm in order to speed up the procedure. The second research area will be the addition of the regulatory element, the purpose of which will be the improvement of the algorithm’s operation in cases of larger input data disturbances and smaller numbers of measurements taken.