Next Article in Journal
Design Optimization of Improved Fractional-Order Cascaded Frequency Controllers for Electric Vehicles and Electrical Power Grids Utilizing Renewable Energy Sources
Next Article in Special Issue
Adaptive Control Design for Euler–Lagrange Systems Using Fixed-Time Fractional Integral Sliding Mode Scheme
Previous Article in Journal
Daughter Coloured Noises: The Legacy of Their Mother White Noises Drawn from Different Probability Distributions
Previous Article in Special Issue
Finite-Time Stabilization of Unstable Orbits in the Fractional Difference Logistic Map
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Analysis of Direct and Inverse Problems for a Fractional Parabolic Integro-Differential Equation

by
Miglena N. Koleva
1,* and
Lubin G. Vulkov
2
1
Department of Mathematics, Faculty of Natural Sciences and Education, University of Ruse “Angel Kanchev”, 8 Studentska Street, 7017 Ruse, Bulgaria
2
Department of Applied Mathematics and Statistics, Faculty of Natural Sciences and Education, University of Ruse “Angel Kanchev”, 8 Studentska Street, 7017 Ruse, Bulgaria
*
Author to whom correspondence should be addressed.
Fractal Fract. 2023, 7(8), 601; https://doi.org/10.3390/fractalfract7080601
Submission received: 27 June 2023 / Revised: 13 July 2023 / Accepted: 2 August 2023 / Published: 4 August 2023

Abstract

:
A mathematical model consisting of weakly coupled time fractional one parabolic PDE and one ODE equations describing dynamical processes in porous media is our physical motivation. As is often performed, by solving analytically the ODE equation, such a system is reduced to an integro-parabolic equation. We focus on the numerical reconstruction of a diffusion coefficient at finite number space-points measurements. The well-posedness of the direct problem is investigated and energy estimates of their solutions are derived. The second order in time and space finite difference approximation of the direct problem is analyzed. The approach of Lagrangian multiplier adjoint equations is utilized to compute the Fréchet derivative of the least-square cost functional. A numerical solution based on the conjugate gradient method (CGM) of the inverse problem is studied. A number of computational examples are discussed.

1. Introduction

The second-order hyperbolic and parabolic systems of partial differential equations, as well as coupled PDE-ODE systems, are a substantial theoretical foundation for many mathematical and engineering problems [1,2,3,4,5,6,7,8,9,10]. In the present paper, we concentrate on the features of a parabolic PDE-ODE system.
Recently, because of the the many applications, fractional integro-differential equations are the focus of a large body of results. The classical integer order derivative is a local operator, which is not adequate for a description of many processes in physics, mechanics, industrial finance, etc. The fractional derivative is a nonlocal operator, which is often used to model phenomena in heat mass transfer [4,5,11], medicine [12], viscoelastic materials [8,13], porous media [1,5,7,10,14,15], mathematical biology and in particularly the honeybee population [16], mathematical finance and economics—see [17,18] and references therein—and atmosphere pollution [6,19].
In [4], analytical exact solutions to the neutron fractional parabolic PDE-ODE system are derived. The inverse problems (IPs) of differential equations with time fractional derivatives concern the recovery of a variety unknown parameters, such as source term, fractional order, coefficients, etc. For instance, Georgiev and Vulkov [16] present a numerical study of the inverse problem for coefficient determination in honeybee population models with Caputo and Caputo–Fabrizio time fractional derivative, and in [13] a parameter identification inverse problem is studied.
In this work, we numerically recover the unknown diffusion coefficient for a time fractional integro-diffusion convection reaction equation, obtained from a model system of differential equations.

A System of a Parabolic and a Ordinary Differential Equation

As a simple physical motivation, we consider a model of non-Darcian flow in transport media used to describe different real processes, see, e.g., [3,5,7,8,10,20].
Suppose that diffusion and solute transport initiate within the mobile phase, while the change in solute concentration within the immobile phase is a dynamic process influenced by the porous media’s low permeability and significant heterogeneity. The fractional calculus has been introduced as an efficient tool for modeling. A general form of a differential equation system that describes such processes reads as follows [3,5,7,20]:
α u t α x a ( x ) u x v ( x ) u + r λ ( x ) ( u w ) = 0 in Q T = Ω × ( 0 , T ) , Ω = ( 0 , l ) ,
α 1 w t α 1 + λ ( x ) w = λ ( x ) u in Q T .
with the initial condition and appropriate boundary conditions. Here, 0 < α , α 1 1 and
α v t α = 1 Γ ( 1 α ) 0 t ( t s ) α v s ( x , s ) d s
is the pointwise Caputo derivative [13,21,22,23]. With u and w are denoted the concentration of a contaminant in the relevant phases. The authors of [24] using Equation (2) write an integral formula for w ( x , t ) and placed it in (1) to obtain an integro-differential equation for u ( x , t ) , ( x , t ) Q T . On this basis, we consider a more general equation:
α u t α = x a ( x , t ) u x + b ( x , t ) u x + c ( x , t ) u + 0 t ρ ( x , t , s ) u ( x , s ) d s + f ( x , t ) = Z u + f .
The remainder of the paper is constructed as follows. In Section 2, we give some notations and discuss the well-posedness of the direct problem. Section 3 is devoted to the numerical solution of the direct problem. In Section 4, we introduce the IP and study its quasisolution. In Section 5, we introduce a finite difference method, combined with the conjugate gradient method (CGM) to solve the IP. Section 6 presents computational results that validate the analysis and demonstrate the efficiency of the proposed algorithm. Some comments and conclusions finalized this paper.

2. Well-Posedness of the Direct Problem

In the IP, it is assumed in advance that the corresponding direct problem is well-posed. The main purpose of this section is the discussion of this question. In this section, our main focus is on the discussion of this question. We consider the Equation (4), 0 < a 0 a ( x , t ) a 1 , ( x , t ) Q ¯ T and with no loss of generality, zero Dirichlet boundary conditions
u ( 0 , t ) = 0 , u ( l , t ) = 0 ,
and initial condition
u ( x , 0 ) = u 0 ( x ) .

2.1. Notations and Preliminaries

In the following, we assume classical interpretations of fractional calculus and equations and will refer to monographs [22,23]. Further, let
L p ( 0 , T ) : = v : 0 T | v ( t ) | p d t <
with p 1 and
W 1 , 1 ( 0 , T ) : = v L 1 ( 0 , T ) : d v d t L 1 ( 0 , T ) ,
where the norms are conventionally established by
v L p ( 0 , T ) : = 0 T | v ( t ) | p d t 1 / p , v W 1 , 1 ( 0 , T ) : = v L 1 ( 0 , T ) + d v d t L 1 ( 0 , T ) .
In view of the Young inequality on the convolution, it could be directly established that the classical Caputo derivative (3) could be well-defined for
v W 1 , 1 ( 0 , T ) and d α v d t α L 1 ( 0 , T ) .
The Riemann–Liouville left and right integrals of fractional order α are defined as follows [22,23]:
I 0 + α u ( t ) = 1 Γ ( α ) 0 t u ( s ) ( t s ) 1 α d s , I T α u ( t ) = 1 Γ ( α ) t T u ( s ) ( s t ) 1 α d s ,
and the corresponding Riemann–Liouville left and right fractional derivatives for 0 < α < 1 are defined as
D 0 + α u ( t ) = 1 Γ ( 1 α ) d d t 0 t u ( s ) ( t s ) α d s = d d t I 0 + 1 α u ( t ) , D T α u ( t ) = 1 Γ ( 1 α ) d d t t T u ( s ) ( s t ) α d s = d d t I T 1 α u ( t ) ,
Further, the notations I 0 + α u and D 0 + α u will be used instead of I 0 + α u ( t ) and D 0 + α u ( t ) , respectively, etc.
Let us recall the integration by parts formula, see, e.g., [25] (Proposition 2.19). Let 0 < α < 1 and assume that w ( t ) A C [ 0 , T ] , w ( t ) C 1 α [ 0 , T ] , and v ( t ) A C 1 α T [ 0 , T ] , where
C α [ 0 , T ] : = y ( t ) : y ( t ) t α C [ 0 , T ] A C α T [ 0 , T ] : = y ( t ) : y ( t ) ( T t ) α A C [ 0 , T ] ,
with C [ 0 , T ] and A C [ 0 , T ] being the classes of continuous and absolutely continuous functions in [ 0 , T ] [22,23].
Then, we have the integration by parts formula
0 T α w t α ( t ) v ( t ) d t = w ( T ) I T 1 α v ( T ) w ( 0 ) I T 1 α v ( 0 ) + 0 T w ( t ) D T α v ( t ) d t .
The subsequent Gronwall-type inequality is established in [26].
Lemma 1.
Let the absolutely continuous function y ( t ) 0 satisfy the inequality
d α y ( t ) t α c 1 y ( t ) + c 2 ( t ) , 0 < α 1
for almost all t in [ 0 , T ] , where c 1 is a positive constant and c 2 ( t ) 0 is an integrable function on [ 0 , T ] . Then,
y ( t ) y ( 0 ) E α ( c 1 t α ) + Γ ( α ) E α , α ( c 1 t α ) D 0 t α c 2 ( t ) ,
where
E α ( z ) = k = 0 z k / Γ ( α k + 1 ) and E α , μ ( z ) = k = 0 z k / Γ ( α k + μ )
are the Mittag–Leffler functions.
Further, we also use the relation between Caputo and Rieamann–Liouville derivatives
d 0 α y ( t ) d t α = D 0 + α y ( t ) y ( 0 ) Γ ( 1 α ) t α , t ( 0 , T ] ,

2.2. Solution of the Direct Problem

Now, following the results of Section 3 in [24] we show first the forward (direct) problem (4)–(6) and present a priory estimates for the corresponding solutions.
Theorem 1.
Suppose that the smoothness conditions hold:
0 < a 0 < a ( x ) a 1 , max | b ( x , t ) | , b x ( x , t ) , | c ( x , t ) | , | ρ ( x , t ) | C ,
and
a ( x , t ) C 1 , 0 ( Q T ) , b ( x , t ) , c ( x , t ) , f ( x , t ) C ( Q T ) .
Then,
u ( x , t ) C 2 , 0 ( Q T ) C 1 , 0 ( Q ¯ T ) , α u ( x , t ) , t α C ( Q T ) .
and the a priori estimate holds
u 0 2 + I 0 + α u x 0 2 C u 0 0 2 + I 0 + α f 0 2 ,
where C > 0 is a constant, which does not depend on u 0 and f.
Here, we use standard notations [23,27]
( f , g ) L 2 = 0 l f g d x , ( f , f ) L 2 = f 0 2 ,
as well as for the Hölder spaces C 1 , 0 , C 2 , 0 .
The proof of Theorem 1 resembles the one of Theorem 1 in [24], namely, first we take the scalar product of Equation (4) with u. Then, we use Lemma 3 from [26] (continuous version of Lemma 3 in the current study). Next, we apply several times the ε -Cauchy inequality to estimate the terms of the scalar product and to obtain a basic inequality. Further, we act with integral operator I 0 + α to both sides of this inequality and obtain (8).
It follows from (8) that the uniqueness and stability of the solution of (4)–(6) in the sense of the Sobolev norm
u 0 2 + I 0 + α u x 0 2 ,
i.e., the well-posedness of the direct problem.

3. Numerical Solution of the Direct Problem

Here we construct a monotone finite difference discretization for (4)–(6) and investigate the stability and convergence.

3.1. Difference Scheme

We construct a second-order monotone Samarskii-type [28] finite difference scheme for the discretization in space and L 2 1 σ scheme [29] of order 3 α for the temporal approximation.
On the rectangle Q ¯ T , we build the uniform mesh
w ¯ h τ = w ¯ h × w ¯ τ = { ( x i , t n ) , x i w ¯ h , t n w ¯ τ } ,
where
w ¯ h = { x i = i h , i = 0 , 1 , , N , N h = l } , w ¯ τ = { t n = n τ , n = 0 , 1 , , M , M τ = T } .
Further, we set u i = u i ( t ) = u ( x i , t ) , U i n = U ( x i , t n ) u ( x i , t n ) .
To build a monotone scheme for (4)–(6), that satisfies the maximum principle for arbitrary mesh step sizes, we consider the perturbation of Equation (4), see [28]
α u t α = L ˜ u + 0 t ρ ( x , t , s ) u ( x , s ) d s + f ( x , t ) , L ˜ u = a ˜ ( x ) x a ( x , t ) u x + b ( x , t ) u x , a ˜ ( x , t ) = 1 1 + R ( x , t ) , R ( x , t ) = h | b ( x , t ) | 2 a ( x , t ) .
We write b ( x , t ) as a sum
b = b + + b , b + = 1 2 ( b + | b | ) 0 , b = 1 2 ( b | b | ) 0 ,
and approximate b u x in space as follows:
b u x i = b a a u x i q i + p i + 1 ( a ) u x , i + q i p i ( a ) u x ¯ , i ,
where
q i + = b + ( x i , t ) a ( x i , t ) 0 , q i = b ( x i , t ) a ( x i , t ) 0 , p i ( a ) = 1 2 ( a ( x i , t ) + a ( x i 1 , t ) ) or p i ( a ) = a ( x i 1 / 2 , t ) , u x , i = u i + 1 u i h , u x ¯ , i = u i u i 1 h .
Thus, the resulting semidiscretization
α u t α = a ˜ i p i u x ¯ , i x , i + q + i p i + 1 u x , i + q i p i u x ¯ , i + f i ( t ) ,
achieves a second-order of convergence in space, as is proved in [28], since
q ± = b ± + O ( h 2 ) , a q ± = b ± , b + b = | b | , ( p ( a ) ) ( + 1 ) u x = a u x + 0.5 h x p ( a ) u x + O ( h 2 ) , ( p ( a ) ) ( + 1 ) = p i + 1 ( a ) , p ( a ) u x ¯ = a u x 0.5 h x a u x + O ( h 2 ) , p u x ¯ x = x k u x + O ( h 2 ) .
The temporal derivative is discretized by L 2 1 σ formula of order O ( τ 3 α ) at ( x i , t n + σ ) , t n + σ = ( n + σ ) τ , σ = 1 α / 2 , n = 0 , 1 , , M 1 , derived in [29].
α u t α D α U t , i : = τ α Γ ( 2 α ) j = 0 n c j ( n , α , σ ) u ( x i , t n j + 1 ) u ( x i , t n j ) ,
where, c 0 ( 0 , α , σ ) = d 0 ( α , σ ) for n = 0 and for n 1
c j ( n , α , σ ) = d 0 ( α , σ ) + b 1 ( α , σ ) , j = 0 , d j ( α , σ ) + b j + 1 ( α , σ ) b j ( α , σ ) , 1 j n 1 , d n ( α , σ ) b n ( α , σ ) , j = n ,
d m ( α , σ ) = σ 1 α , m = 0 , ( l + σ ) 1 α ( m 1 + σ ) 1 α , m 1 , b m ( α , σ ) = 1 2 α ( m + σ ) 2 α ( m 1 + σ ) 2 α 1 2 ( m + σ ) 1 α + ( m 1 + σ ) 1 α ,   m 1 .
Using (9), the σ —weighted finite difference discretization of the Equation (4) is
D c α U t , i = σ Λ ˜ U i n + 1 + ( 1 σ ) Λ ˜ U i n + c i n + σ σ U i n + 1 + ( 1 σ ) U i n + I i n + σ + f i n + σ ,
i = 1 , 2 , , N 1 , where
Λ ˜ U i n + 1 = a ˜ i n + σ p i n + σ U x ¯ , i n + 1 x , i + q + i n + σ p i + 1 n + σ U x , i n + 1 + q i n + σ p i n + σ U x ¯ , i n + 1 , Λ ˜ U i n = 1 h a i + 1 / 2 n U x , i n a i 1 / 2 n U x ¯ , i n + b i n + σ U x , i n + U x ¯ , i n 2 , p i n + σ = a i n + σ + a i 1 n + σ 2 or p i n + σ = a i 1 / 2 n + σ , a i n + σ = a ( x i , t n + σ ) , a i ± 1 / 2 n = a ( x i ± h / 2 , t n )
and similarly for b i n + σ , ( q ± ) i n + σ and f i n + σ .
By I i n in (10), we denote the approximation of the integral in the Equation (4), namely,
I i n + σ 0 t n + σ ρ ( x i , t n + σ , s ) u ( x i , s ) d s = j = 1 n t j 1 t j ρ ( x i , t n + σ , s ) u ( x i , s ) d s + t n t n + σ ρ ( x i , t n + σ , s ) u ( x i , s ) d s j = 1 n u ( x i , t j 1 / 2 ) t j 1 t j ρ ( x , t n + σ , s ) d s + u ( x i , t n ) t n t n + σ ρ ( x , t n + σ , s ) d s j = 1 n U i j 1 + U i j 2 t j 1 t j ρ ( x i , t n + σ , s ) d s + U i n t n t n + σ ρ ( x i , t n + σ , s ) d s .
The approximation (11) is obtained in the same fashion as the discretization of the fractional integral operator.
Further, the integrals in (11) can be computed by the trapezoidal rule, midpoint rule or exactly, depending on the function ρ .
The Dirichlet boundary conditions (5) and initial condition (6) are incorporated into the numerical scheme straightforwardly:
σ U 0 n + 1 + ( 1 σ ) U 0 n = 0 , σ U N n + 1 + ( 1 σ ) U N n = 0 , n = 1 , 2 , , M , U i 0 = u 0 ( x i ) , i = 0 , 1 , , N .
Lemma 2
([30]). Suppose that the non-negative U n , φ n , n = 0 , 1 , fulfill the inequalities
D t j + τ α U n ε 1 U n + 1 + ε 2 U n + φ n , n 1 ,
where ε 1 > ε 2 0 are constants. Then, there exists τ * , such that: for τ τ * , we have
U n 2 U 0 + t n α Γ ( 1 + α ) max 0 n ˜ n φ n ˜ E α ( 2 ϵ t n α ) , 1 n M ,
where ϵ = ϵ 1 + ϵ 2 ( 2 2 1 α ) 1 .
Further, we use the difference scheme inequality as well.
Lemma 3
([29]). For every mesh function U ( t ) , defined on w ¯ τ , the inequality holds:
U ( σ ) D t j + σ α U 1 2 D t j + σ α ( U 2 ) ,
where U ( σ ) = σ U n + 1 + ( 1 σ ) U , where U = U n .

3.2. Solvability and Convergence

In the following theorem, we present the energy estimate of the difference scheme (10)–(12).
Theorem 2.
Let the conditions Theorem 1 be satisfied. Then, there exists τ * such that, if τ τ * , for the problem (10), (12) the a priori estimate holds:
U n + 1 0 C ( U 0 0 2 + max F n ˜ 0 2 ) ,
where the positive constant C is independent on h and τ.
Proof. 
We apply the energy method of Samarskii [28] for the case d = I 0 . With this aim, we introduce the scalar product along with its corresponding norms:
( U , V ) = i = 1 N 1 U i V i h , ( U , V ] = i = 1 N U i V i h , ( U , U ) = ( 1 , U 2 ) = U 0 2 .
We obtain the scalar product of U ( σ ) with (10):
D t n + σ α U , U ( σ ) = a ˜ ( p ( a ) U x ¯ ( σ ) ) x , U ( σ ) + q + ( p ( a ) ) ( + 1 ) U x ( σ ) , U ( σ ) + q p ( a ) U x ( σ ) , U ( σ ) + m = 0 n + σ ρ ˜ i , m n U i m τ ˜ , U ( σ ) + F , U ( σ ) .
We estimate each term of the equality (14) using Lemma 3:
D t n + σ α U , U ( σ ) 1 2 1 , D t n + σ α U 2 1 2 D t n + σ α U 0 2 , a ˜ p ( a ) U x ¯ ( σ ) x , U ( σ ) = a ˜ p ( a ) U x ( σ ) U ( σ ) | x = x 0 x = x N p ( a ) U x ¯ ( σ ) , ( a ˜ U ( σ ) ) x ¯ p ( a ) a ˜ x ¯ U x ¯ ( σ ) U ( σ ) p ( a ) a ˜ ( 1 ) , ( U x ¯ ( σ ) ) 2 1 1 + h C 1 p ( a ) a ˜ , ( U x ¯ ( σ ) ) 2 + ε U x ¯ ( σ ) ] | 0 2 + C 2 ( ε ) U ( σ ) 0 2 , q p ( a ) U x ¯ ( σ ) , U ( σ ) + q + ( p ( a ) ) ( + 1 ) U x ( σ ) , U ( σ ) ε U x ¯ ( σ ) ] | 0 2 + C 3 ( ε ) U ( σ ) 0 2 , m = 0 n + σ ρ ˜ i , m n U i m τ ˜ , U ( σ ) 1 2 U ( σ ) 2 + m = 0 n + σ ρ ˜ i , m n U i m τ ¯ 2 1 2 U ( σ ) 0 2 + 1 2 m = 0 n + σ ρ ˜ i , m n m = 0 n + σ U m 2 τ ¯ 1 2 U ( σ ) 0 2 + C 4 m = 0 n + σ ( 1 , U m 2 ) τ ˜ 1 2 U ( σ ) 0 2 + C 4 m = 0 n + σ U m 0 2 τ ¯ , F , U ( σ ) 1 2 U ( σ ) 0 2 + 1 2 F 0 2 .
Finally, inserting the last estimates in (14), we obtain
D t n + σ α U 0 2 + U x ¯ ( σ ) 0 2 ε C 5 U x ¯ ( σ ) ] | 0 2 + C 6 ( ε ) U ( σ ) 0 2 + C 7 m = 0 n + σ U m 0 2 τ ¯ + C 8 F 0 2 .
Then, taking ε = 1 / ( 2 C 5 ) , we obtain
D t n + σ α U 0 2 + 1 2 U x ¯ ( σ ) 0 2 C 9 U ( σ ) 0 2 + C 10 m = 0 n + σ U m 0 2 τ ¯ + C 11 F 0 2 .
In view of
m = 0 n + σ U m 0 2 τ ¯ = m = 0 n U m 0 2 τ ¯ + 0.5 τ U n 0 2 ,
we rewrite (15) as follows:
D t n + σ α U 0 2 C 12 ( σ ) U n + 1 0 2 + C 13 ( σ ) U n 0 2 + C 14 G n ,
where
G n = m = 0 n U m 0 2 τ ¯ + F 0 2 .
Now, the application of Lemma 1 to the last inequality results in
U j + 1 0 2 C 15 U 0 0 2 + t n α Γ ( 1 + α ) max 0 n ˜ n G n ˜ ,
where the constant C 15 > 0 does not depend on h and τ .
In view of the expression for G n , we restate the last inequality in the following:
U n + 1 0 2 C 15 U 0 0 2 + max 0 n ˜ n m = 0 n ˜ U n 0 2 τ ¯ + F n ˜ 0 2 .
Letting z n + 1 = max 0 n ˜ n U n ˜ 0 2 we obtain from the last inequality
z n + 1 C 16 m = 0 n z m τ + C 17 G 1 n ,
where
G 1 n = U 0 0 2 + t n α Γ ( 1 + α ) max F n ˜ 0 2 .
The application of the Gronwall lemma from [S, Ch.3] to (16) gives the a priori estimate (13). □
The inequality (13) provides the unicity and stability of the difference scheme (10), (12) with respect to the right-hand side and initial data.
Further, again using Theorem 1, we will obtain an accuracy estimate for the difference problem (10), (12). Let V i n = U i n u i n , where u i n = u ( x i , t n ) and u ( x , t ) is the solution of the differential problem (4)–(6), while U ( x i , t n ) = U i n ˙ is the solution of the difference problem (10), (12). Then, inserting U = V + u in (10), (12), for the mesh function V, we derive the following discrete problem:
D t n + σ α V = L ˜ V ( σ ) + m = 1 n + σ ρ ˜ i m n V i m + φ i n , ( x , t ) w h , τ , V 0 ( σ ) = V N ( σ ) = 0 , V ( x , 0 ) = 0 ,
with approximation error φ = O ( h 2 + τ 2 ) .
Applying Theorem 1 to the problem (17), we find
V n + 1 0 2 C max 0 n ˜ n φ n ˜ ,
where the constant C > 0 is independent on h and τ . From this inequality follows the convergence of the discretization (10), (12) towards the solution of the differential problem (4)–(6) on the each time level, so that there exists τ * s.t. for τ τ 0 the estimate holds U n + 1 u n + 1 0 C ( h 2 + τ 2 ) .

4. Quasi Solution of the IP

For the ease of simplicity we further consider space space-dependent coefficient a = a ( x ) in the Equation (4). Then, the IP is to reconstruct a ( x ) under the measurements
u ( x i * , t ) = g i ( t ) , i = 1 , , I , x i * ( 0 , l ) .
We employ the Lagrange multiplier technique [2,27,31] to obtain the Fréchet gradient of the least-squares discrepancy functional associated with the quasisolution of the IP (4)–(6), (18). Because of the instability of their solutions with respect to (18), these IPs are ill-posed [2,27,31].
Moreover, in our case, the solution’s uniqueness of the IP is dependent on the placement of the points x i * in (18) along with the initial and boundary conditions.
Consider the operator form A ( a ) = g of the problem (4)–(6), (18). Here, A : A G denotes an injective operator, a A = { a 0 a ( x ) a 1 , a C ( Ω ¯ ) } represents the admissible set and g = ( g 1 ( t ) , , g I ( t ) ) G , where G is Euclidian space. Our IP is reformulated as a minimization problem:
a * = argmin a A J ( a ) , J ( a ) = i = 1 I 0 T ( u ( x i * , t ; a ) g i ( t ) ) 2 d t ,
with u ( x * , t ; a ) being the solution to the problem (4)–(6).
In order to numerically investigate the recovering of the diffusion coefficient a ( x ) from the observations at interior points, we employ the conjugate gradient method (CGM). It is grounded in the utilization of the Fréchet gradient of the objective functional (19). The increment of the functional δ J ( a ) = J ( a + δ a ) J ( a ) could be express as follows [2,27,31,32]:
δ J ( a ) = U , δ a H ,
where with · , · we denote the scalar product in H = L 2 ( Ω ) . Therefore, J ( a ) = U . We use the adjoint operator
Z * ψ : = x a ( x , t ) ψ x b ( x , t ) ψ + c ( x , t ) b ( x , t ) x ψ + t T ρ ( x , s , t ) ψ ( x , s ) d s .
Theorem 3.
The functional J ( a ) in (19) is Fréchet differentiable and its gradient is defined by
J ( a ) = 0 T ψ x u x d t ,
where ψ ( x , t ) is the solution of the backward adjoint problem to the forward problem (4)–(6), (18), i.e.,
D T α ψ = Z * ψ + i = 1 I 2 ( u g i ) δ ( x x i * ) i n Q T ,
I T 1 α ψ ( x , T ) = 0 , x Ω ,
and
ψ ( 0 , t ) = 0 , ψ ( l , t ) = 0 , t ( 0 , T ) .
The proof resembles the one of Theorem 3 in [24]. A key role has the so-called sensitivity problem. By u ( x , t ; a ) : = u ( x , t ) we denote that u ( x , t ) depends on the parameter a ( x ) . Assume that for a A and a + δ a A , δ u ( x , t ; δ a ) : = u ( x , t ; a + δ a ) u ( x , t ; a ) .
Hence, the deviation δ u : = δ u ( x , t ; δ a ) fulfills the next IBVP with an accuracy up to terms of order o ( | δ a | ) 2 :
α δ u t α = Z ( δ u ) + x δ a u x i n Q T ,
δ u ( x , 0 ) = 0 , x Ω ,
δ u ( 0 , t ) = 0 , δ u ( l , t ) = 0 , t ( 0 , T ) .
Next, for the minimization of the functional (19), we apply the Lagrange multiplier method [2,27,31] and follow the procedure outlined in [24] to obtain the gradient (20).

5. Conjugate Gradient Steepest Descent Method

We explore the CGM [2,24,27,31,33] in developing a numerical algorithm to solve the IP for recovering the diffusion coefficient from interior concentration observations in porous media. As is common knowledge, the Fréchet gradient (20) plays a key role.
The iterative procedure is implemented by applying the CGM (see, e.g., [2,27,31,33])
a k + 1 ( x ) = a k ( x ) + β k d k ( x ) , k = 0 , 1 , ,
where k denotes the number of iteration, a k and a k + 1 correspond to consecutive approximations of the minimizer, and d k ( x ) and β k refer to the search direction and the search step size, respectively.
At each step, the search direction d k ( x ) is computed as a linear combination of the steepest descent direction at the current approximation and the search direction from previous iterations, namely,
d 0 ( x ) = J [ a 0 ( x ) ] , d k ( x ) = J [ a k ( x ) ] + γ k d k 1 ( x ) , i = 1 , 2 , ,
with γ k being the conjugation coefficient. An often-used choice of γ k is as follows [2,27,32]:
γ 0 = 0 , γ k = 0 l | J [ a k ( x ) ] | 2 d x 0 l | J [ a k 1 ( x ) ] | 2 d x , k = 1 , 2 ,
first proposed by Fletcher and Reeves [34]. The coefficient β k is determined through the use of line search, which involves minimizing along the specified search direction, β k = arg min β J ( a k ( x ) + β d k ( x ) ) . Applying the sensitivity problem (24)–(26), a method similar to [35] can be employed to linearize the objective functional J ( a k ( x ) + β d k ( x ) ) with respect to β and obtain the search step size
β k = i = 1 I 0 t ( u k x i * , t ; a k ) g i ( t ) ) δ u k x i * , t ; a k ) i = 1 I 0 t [ δ u k ( x i * , t ; a k ) ] 2 d t ,
where u k and δ u k are the corresponding direct and sensitivity problems solutions and a ( x ) = a k ( x ) , δ a = d k ( x ) .

6. Numerically Solving the Inverse Problem

In this section, we will develop numerical discretizations for solving the IP. To this aim, we need to approximate adjoint problem (21)–(23), sensitivity problem (24)–(26), gradient of the functional (20), etc.

6.1. Discretization of the Adjoint Problem

Now, we consider (21)–(23). Applying the temporal inversion t * = T t in (21), we obtain a forward problem with an initial instead of terminal condition and ψ ( x , t ) = ψ ( x , T t * ) = ψ * ( x , t * ) , a ( x , T t ) = a * ( x , t * ) , b ( x , T t ) = b * ( x , t * ) , c ( x , T t ) = c * ( x , t * ) , u ( x , T t ) = u * ( x , t * ) . With the fractional derivative D T α ψ , we deal as in [14], applying also the variable change s * = T s , namely,
D T α ψ ( x i , t ) = 1 Γ ( 1 α ) d d t t T ψ ( x i , s ) ( s t ) α i d s = 1 Γ ( 1 α i ) d d t * t * 0 ψ * ( x i , s * ) ( T s * ( T t * ) ) α i d s * = 1 Γ ( 1 α i ) d d t * 0 t * ψ * ( x i , s * ) ( t * s * ) α i d s * = D 0 + α ψ * ( x i , t * ) .
Next, in view of (7), considering ψ i ( x , T ) = ψ i * ( x , 0 ) , we obtain
D 0 + α i ψ * ( x i , t * ) = α ψ * ( x i , t * ) ( t * ) α + ψ * ( x i , 0 ) Γ ( 1 α ) ( t * ) α .
We take ψ * ( x i , 0 ) = 0 , since in this case (22), (23) are satisfied for t * = T t .
We treat the time integral in (21) similarly
t T ρ ( x , t , s ) u ( x , s ) d s = 0 t * ρ * ( x , t * , s * ) u * ( x , s * ) d s * ,
where ρ * ( x , t * , s * ) = ρ ( x , T t * , T s * ) , u * ( x , s * ) = u ( x , T s * ) .
Then, we approximate the adjoint problem in the same fashion as the direct problem and derive
D c α ( Ψ * ) t * , i = σ Λ ˜ * ( Ψ * ) i n + 1 + ( 1 σ ) Λ ˜ * ( Ψ * ) i n + ( c ˜ * ) i n + σ σ ( Ψ * ) i n + 1 + ( 1 σ ) ( Ψ * ) i n + ( I * ) i n + σ + ( f * ) i n + σ , i = 1 , 2 , , N 1 ,
where ( y * ) i n + σ = y * ( x i , t n + σ * ) , c ˜ * = c * 2 b * x , ( Ψ * ) i n = ψ * ( x i , t n ) and
Λ ˜ * ( Ψ * ) i n + 1 = ( a ˜ * ) i n + σ ( p * ) i n + σ ( Ψ * ) x ¯ , i n + 1 x , i ( q * ) + i n + σ ( p * ) i + 1 n + σ ( Ψ * ) x , i n + 1 ( q * ) i n + σ ( p * ) i n + σ ( Ψ * ) x ¯ , i n + 1 , Λ ˜ * ( Ψ * ) i n = 1 h ( a * ) i + 1 / 2 n ( Ψ * ) x , i n ( a * ) i 1 / 2 n ( Ψ * ) x ¯ , i n ( b * ) i n + σ ( Ψ * ) x , i n + ( Ψ * ) x ¯ , i n 2 , ( I * ) i n + σ = j = 1 n ( Ψ * ) i j 1 + ( Ψ * ) i j 2 t j 1 * t j * ρ ( x , t n + σ * , s * ) d s * + ( Ψ * ) i j t n * t n + σ * ρ * ( x , t n + σ * , s * ) d s * .
For Dirichlet boundary conditions (23), we have
σ ( Ψ * ) 0 n + 1 + ( 1 σ ) ( Ψ * ) 0 n = 0 , σ ( Ψ * ) N n + 1 + ( 1 σ ) ( Ψ * ) N n = 0 .

6.2. Discretization of the Sensitivity Problem

To obtain numerical the solution δ U i n of the sensitivity problem (24)–(26), we proceed in the same manner as for the direct problem:
D c α δ U t , i = σ Λ ˜ δ U i n + 1 + ( 1 σ ) Λ ˜ δ U i n + ( c ) i n + σ σ δ U i n + 1 + ( 1 σ ) δ U i n + δ a i + 1 / 2 h σ δ U x , i n + 1 + ( 1 σ ) δ U x , i n δ a i 1 / 2 h σ δ U x ¯ , i n + 1 + ( 1 σ ) δ U x ¯ , i n , σ δ U 0 n + 1 + ( 1 σ ) δ U 0 n = 0 , σ δ U N n + 1 + ( 1 σ ) δ U N n = 0 , δ U i 0 = 0 , i = 0 , 1 , , N ,
where δ a i 1 / 2 = 0.5 ( a i 1 + a i ) .

6.3. Discretization of the Gradient of the Functional and γ k

For the approximation of the gradient of the functional at spatial grid node x i , i = 1 , 2 , , N 1 , we apply central second order discretization and the trapezoidal rule:
J ( a k ) ( x i ) = j = 1 M t j 1 t j ψ ( x i , s ) x u ( x i , s ) x d s 1 4 h 2 j = 1 M t j 1 t j Ψ ˜ i ( s ) U ˜ i ( s ) d s τ 8 h 2 Ψ ˜ i 0 U ˜ i 0 + 2 j = 1 M Ψ ˜ i j U ˜ i j + Ψ ˜ i M 1 U ˜ i M ,
where Ψ ˜ i = Ψ i + 1 Ψ i 1 and U ˜ i = U i + 1 U i 1 .
In order to obtain the conjugate coefficient at iteration number k > 0 , we need to solve the integral
0 l J ( a k ) = 0 l 0 T ψ x u x d s d x τ 16 h l = 0 M 1 i = 0 N 1 ( Ψ ˜ i j U ˜ i j + Ψ ˜ i + 1 j U ˜ i + 1 j + Ψ ˜ i j + 1 U ˜ i j + 1 + Ψ ˜ i + 1 j + 1 U ˜ i + 1 j + 1 ) .

6.4. Realization

The CGM for solving the IP is realized by the Algorithm 1.
Algorithm 1 Inverse problem
  • Require:  a 0 , ε
  • Ensure:  a * , U ( x i , t n ; a * )
  • k 0
     Calculate J ( a k ) , using (19)
  • while J ( a k ) J ( a k 1 ) > ε J  do
  •    U ( x i , t n ; a k ) direct problem (10), (12)
  •    Ψ ( x i , t n ; a k ) adjoint problem (30), (31)
      Calculate J ( a k ) , using (33)
      Calculate d k , using (28)
  •    δ a k + 1 d k
  •    δ U ( x i , t n ; a k ) sensitivity problem (32)
      Calculate β k , using (29)
  •    a k + 1 a k + β k d k in view of (27)
  •    k k + 1
  • end while
  • a * a k + 1
  • U ( x i , t n ; a * ) (10), (12)

7. Computational Results

Here, we provide numerical results from computations for both direct and inverse problems in order to illustrate the effectiveness and accuracy of the suggested numerical method. The computational domain is [ 0 , 1 ] × [ 0 , 1 ] .
We consider (4)–(6), with r.h.s f ( x , t ) , initial and boundary conditions, such that u ( x , t ) = t 4 + α sin ( π x ) , is the exact solution. The other model parameters are
a ( x ) = 2 x 2 + 3 , b ( x , t ) = x + t , c ( x , t ) = 3 x + t 2 , ρ ( x , t ) = 2 x + t 2 , u 0 ( t ) = u 1 ( t ) = 0 .
Example 1 
(Direct problem). In this test, we check the accuracy of the numerical scheme (10), (12) for the direct problem. We give errors and convergence rates in maximal norm for N = M , namely,
E N = max 0 i N max 0 n N | u ( x i , t n ) U i n | , CR h = log 2 E N E 2 N .
In Table 1, we illustrate the order of convergence in space and time. As a benchmark test, we also include the case of integer time derivative, namely α = 1 . The computations are conducted for τ = h . Results show that the order of convergence is O ( τ 2 + h 2 ) .
Example 2 
(Inverse problem). With this test, we demonstrate the efficiency in recovering the diffusion coefficient a ( x ) and solution u of the IP. The initial guess for the diffusion coefficient is generated by [36]
a 0 = a ( x ) + 2 ρ a a ( x ) ( ϱ ( x ) 0.5 ) ,
where ϱ ( x ) is a uniformly distributed random function in the interval [ 0 , 1 ] and a ( x ) is the exact diffusion coefficient. In the same fashion, we generate noisy data
g i ( t n ) = u ( x i , t n ) + 2 ρ u u ( x i , t n ) ( ϱ ( x ) 0.5 ) , i = 1 , 2 , , I , n = 1 , 2 , , M ,
where u ( x i , t n ) is the exact solution at points of measurements ( x i , t n ) .
The measurements are located at grid nodes; namely, we take seven spatial nodes I = 7 , at x i * = [ 0.1 , 0.2 , 0.5 , 0.6 , 0.7 , 0.8 , 0.9 ] and all time layers. The initial data are regularized by smoothing the function a 0 using fifth-degree polynomial curve fitting. The stopping accuracy is ε J = 10 3 . Actually, for all presented experiments at the first 2–3 iterations, the functional J decreases more rapidly and then, the decreasing is very small. On Figure 1 and Figure 2, we plot the true (exact) and restored diffusion coefficient a ( x ) and error | u U | of the computed by the inverse problem solution U for N = M = 80 , α = 0.5 and α = 0.9 , ρ a = 0.05 , ρ u = 0.01 .
In Figure 3, we depict the exact and restored diffusion a ( x ) for N = M = 80 , α = 0.25 and α = 0.95 , ρ a = 0.1 and ρ u = 0.05 . We observe that for both values of the fractional order, the accuracy of the restored function is very close.
The results show that the main factor for the accuracy of the solution is the deviation of the measurements, rather than the fractional order.
Next, we take initial function, right-hand side
u 0 ( x ) = 2 sin ( π x ) , f ( x , t ) = 5 Γ ( 1 + α 1 ) + Γ ( 1 + 2 α 1 ) Γ ( 1 + α 1 ) t α 1 + Γ ( 1 + 3 α 1 ) Γ ( 1 + 2 α 1 ) t 2 α 1 sin ( π x )
and generate perturbed data from (34), (35), replacing the exact solution u with a numerical solution of the direct problem, computed for the exact diffusion coefficient. Then, the solution of the discrete inverse problem U will be compared with the corresponding numerical solution of the direct problem. The computations are performed for N = M = 80 and the same points of measurements as in previous test. On Figure 4 (left), we illustrate the impact of the fractional order on the accuracy for ρ a = 0.1 , ρ u = 0.08 . As before, we deduce that it is not so significant. On Figure 4 (right), we examine the impact of the perturbation on the accuracy for α = 0.5 . We observe that the deviation has a larger impact on accuracy than the fractional order. On Figure 5, we plot the error of the solution for α = 0.5 and two sets of perturbation ρ a = ρ u = 0.01 and ρ a = ρ u = 0.1 . In spite of the noise affecting the accuracy of the solution, it is evident that the precision of the solution is quite satisfactory when employing the recovered diffusion coefficient.
The computational simulations, presented in this section to test the proposed approach are implemented by MATLAB® R2022a.

8. Conclusions

In this paper, we have solved a single fractional parabolic integro-differential equation, derived by reducing a coupled fractional parabolic PDE-ODE system. First, we discussed the well-posedness of the Dirichlet problem and then we performed an analysis for stability and convergence of Alikhanov’s scheme [29]. We have developed an approach for numerical recovering of the space-dependent diffusion coefficient using a finite number of space points observations. We have experimentally illustrated the efficacy of the suggested method. The order of the fractional derivative does not significantly affect the results. The numerical tests with perturbed data showed the accurate recovery of the diffusion coefficient with a satisfactory error.
We intend to expand the investigation presented in this paper to other inverse problems of type [37,38,39,40,41,42] for one and multidimensional systems of fractional differential equations.

Author Contributions

Conceptualization, L.G.V.; methodology, M.N.K. and L.G.V.; investigation, M.N.K. and L.G.V.; resources, M.N.K. and L.G.V.; writing—original draft preparation, M.N.K. and L.G.V.; writing—review and editing, L.G.V.; visualization, M.N.K. All authors have read and agreed to the published version of the manuscript.

Funding

This research was supported by the Bulgarian National Science Fund under the Project KP-06-N 62/3 “Numerical methods for inverse problems in evolutionary differential equations with applications to mathematical finance, heat-mass transfer, honeybee population and environmental pollution”, 2022.

Data Availability Statement

Not applicable.

Acknowledgments

The authors express their gratitude to the anonymous reviewers, whose valuable comments and suggestions have improved the quality of the paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Alaimo, G.; Piccolo, V.; Cutolo, A.; Deseri, L.; Fraldi, M.; Zingales, M. A fractional order theory of poroelasticity. Mech. Res. Commun. 2019, 100, 103395. [Google Scholar] [CrossRef]
  2. Alifanov, O.M.; Artioukhine, E.A.; Rumyantsev, S.V. Extreme Methods for Solving Ill-Posed Problems with Applications to Inverse Heat Transfer Problems; Begell House: Danbury, CT, USA, 1995. [Google Scholar]
  3. Armstrong, J.E.; Frind, E.O.; McClellan, R.D. Nonequilibrium mass transfer between the vapor, aqueous, and solid phases in unsaturated soils during vapor extraction. Water Resour. Res. 1994, 30, 355–368. [Google Scholar] [CrossRef]
  4. Burqan, A.; Shqair, M.; El-Ajou, A.; Ismaeel, S.M.E.; AlZhour, Z. Analytical solutions to the coupled fractional neutron diffusion equations with delayed neutrons system using Laplace transform method. AIMS Math. 2023, 8, 19297–19312. [Google Scholar] [CrossRef]
  5. Gao, G.; Zhan, H.; Feng, S.; Fu, B.; Ma, Y.; Huang, G. A new mobile-immobile model for reactive solute transport with scale-dependent dispersion. Water Resour. Res. 2010, 46, 1–16. [Google Scholar] [CrossRef]
  6. Goulart, A.G.; Lazo, M.J.; Suarez, J.M.S. A new parameterization for the concentration flux using the fractional calculus to model the dispersion of contaminants in the planetary boundary layer. Phys. A Stat. Mech. Appl. 2019, 518, 38–49. [Google Scholar] [CrossRef] [Green Version]
  7. Li, X.; Wen, Z.; Zhu, Q.; Jakada, H. A mobile-immobile model for reactive solute transport in a radial two-zone confined aquifer. J. Hydrol. 2020, 580, 124347. [Google Scholar] [CrossRef]
  8. Mahiudin, M.; Godhani, D.; Feng, L.; Liu, F.; Langrish, T.; Karim, M. Application of Caputo fractional rheological model to determine the viscoelastic and mechanical properties of fruit and vegetables. J. Postharvest Biol. Technol. 2020, 163, 111147. [Google Scholar] [CrossRef]
  9. Vougalter, V. Inverse problems for some systems of parabolic equations with coefficient depending on time. In Mathematical Methods in Modern Complexity Science; Springer International Publishing: Berlin/Heidelberg, Germany, 2021; pp. 189–197. [Google Scholar]
  10. Zhou, H.; Yang, S.; Zhang, S. Modeling non-darcian flow and solute transport in porous media with the caputo–fabrizio derivative. Appl. Math. Model. 2019, 68, 603–615. [Google Scholar] [CrossRef]
  11. Gyulov, T.B.; Vulkov, L.G. Reconstruction of the lumped water-to-air mass transfer coefficient from final time or time-averaged concentration measurement in a model porous media. In Proceedings of the 16th Annual Meeting of the Bulgarian Section of SIAM, Sofia, Bulgaria, 21–23 December 2021. [Google Scholar]
  12. Georgiev, S.; Vulkov, L. Numerical coefficient reconstruction of time-depending integer- and fractional-order SIR models for economic analysis of COVID-19. Mathematics 2022, 10, 4247. [Google Scholar] [CrossRef]
  13. Caputo, M.; Plastino, W. Diffusion in porous layers with memory. Geophys. J. Int. 2004, 158, 385–396. [Google Scholar] [CrossRef] [Green Version]
  14. Koleva, M.N.; Vulkov, L.G. Parameters Estimation in a Time-Fractiona Parabolic System of Porous Media. Fractal Fract. 2023, 7, 443. [Google Scholar] [CrossRef]
  15. Shen, S.; Weizhong Dai, W.; Cheng, J. Fractional parabolic two-step model and its accurate numerical scheme for nanoscale heat conduction. J. Comput. Appl. Math. 2020, 375, 112812. [Google Scholar] [CrossRef]
  16. Georgiev, S.; Vulkov, L. Parameters identification and numerical simulation for a fractional model of honeybee population dynamics. Fractal Fract. 2023, 7, 311. [Google Scholar] [CrossRef]
  17. Koleva, M.N.; Vulkov, L.G. Fast Positivity Preserving Numerical Method for Time-Fractional Regime-Switching Option Pricing Problem. In Advanced Computing in Industrial Mathematics. BGSIAM 2020; Studies 295 in Computational Intelligence; Georgiev, I., Kostadinov, H., Lilkova, E., Eds.; Springer: Cham, Switzerland, 2023; Volume 1076. [Google Scholar]
  18. Tarasov, V.E. (Ed.) Special Issue “Fractional Calculus in Economics and Finance”. Fractal Fract. 2018. Available online: https://www.mdpi.com/journal/fractalfract/special_issues/Fractional_Calculus_Economics_Finance (accessed on 26 June 2023).
  19. Kandilarov, J.; Vulkov, L. Determination of concentration source in a fractional derivative model of atmospheric pollution. AIP Conf. Proc. 2021, 2333, 090014. [Google Scholar] [CrossRef]
  20. Liu, W.; Li, G.; Jia, X. Numerical simulation for a fractal MIM model for solute transport in porous media. J. Math. Res. 2021, 13, 31–44. [Google Scholar] [CrossRef]
  21. Jovanović, B.S.; Vulkov, L.G.; Delić, A. Boundary value problems for fractional PDE and their numerical approximation. In Numerical Analysis and Its Applications; Lecture Notes in Computer Science; Springer: Berlin/Heidelberg, Germany, 2013; Volume 8236, pp. 38–49. [Google Scholar]
  22. Kilbas, A.A.; Srivastava, H.M.; Trujillo, J.J. Theory and Applications of Fractional Differential Equations; North-Holland Mathematics Studies; Elsevier: Amsterdam, The Netherlands, 2006; Volume 204. [Google Scholar]
  23. Podlubny, I. Fractional Differential Rquations; Academic Theory and Applications of Fractional Differential Equations; Elsevier: Amsterdam, The Netherlands, 1998. [Google Scholar]
  24. Gyulov, T.; Vulkov, L. Gradient Optimization in Reconstruction of the Diffusion Coefficient in a Time Fractional Integro-Differential Equation of Pollution in Porous Media. In Modelling and Development of Intelligent Systems MDIS; Communications in Computer and Information Science; Simian, D., Stoica, L.F., Eds.; Springer: Cham, Switzerland, 2022; Volume 1761. [Google Scholar]
  25. Wei, T.; Xian, J. Variational method for a backward problem for a time-fractional diffusion equation. ESAIM M2AN 2019, 53, 1223–1244. [Google Scholar] [CrossRef] [Green Version]
  26. Alikhanov, A.A. A priori estimates for solutions of boundary value problems for fractional-order equations. Differ. Equ. 2010, 46, 660–666. [Google Scholar] [CrossRef] [Green Version]
  27. Hasanov, A.H.; Romanov, V.G. Introduction to Inverse Problems for Differential Equations; Springer: Cham, Switzerland, 2017. [Google Scholar] [CrossRef]
  28. Samarskii, A.A. The Theory of Difference Schemes; Marcel Dekker: New York, NY, USA, 2001. [Google Scholar]
  29. Alikhanov, A.A. A new difference scheme for the time fractional diffusion equation. J. Comput. Phys. 2015, 280, 424–438. [Google Scholar] [CrossRef] [Green Version]
  30. Li, D.; Liao, H.-L.; Sun, W.; Wang, J.; Zhang, J. Analysis of L1 Galerkin FEMs for time-fractional nonlinear parabolic problems. Commun. Comput. Phys. 2018, 24, 86–103. [Google Scholar] [CrossRef] [Green Version]
  31. Chavent, G. Nonlinear Least Squares for Inverse Problems: Theoretical Foundations and Step-by-Step Guide for Applications; Scientific Computation; Springer: Dordrecht, The Netherlands, 2010. [Google Scholar] [CrossRef] [Green Version]
  32. Hanke, M. Conjugate Gradient Type Methods for Ill-Posed Problems, 1st ed.; Chapman and Hall/CRC: New York, NY, USA, 1995. [Google Scholar]
  33. Vabishchevich, P.N.; Denisenko, A.Y. Numerical methods for solving the coefficient inverse problem. Comput. Math. Model. 1992, 3, 261–267. [Google Scholar] [CrossRef]
  34. Fletcher, R.; Reeves, C.M. Function minimization by conjugate gradients. Comput. J. 1964, 7, 149–154. [Google Scholar] [CrossRef] [Green Version]
  35. Fakhraie, M.; Shidfar, A.; Garshasbi, M. A computational procedure for estimation of an unknown coefficient in an inverse boundary value problem. Appl. Math. Comput. 2007, 187, 1120–1125. [Google Scholar] [CrossRef]
  36. Samarskii, A.A.; Vabishchevich, P.N. Numerical Methods for Solving Inverse Problems in Mathematical Physics; De Gruyter: Berlin, Germany, 2007; 438p. [Google Scholar]
  37. Durdiev, D.K.; Zhumaev, Z.Z. One-dimensional inverse problems of finding the kernel of integrodifferential heat equation in a bounded domain. Ukr. Math. J. 2022, 73, 1723–1740. [Google Scholar] [CrossRef]
  38. Li, Z.; Yamamoto, M. Uniqueness for inverse problems of determining orders of multi-term time-fractional derivatives of diffusion equations. Appl. Anal. 2015, 94, 570–579. [Google Scholar] [CrossRef] [Green Version]
  39. Li, Z.; Huang, X.; Yamamoto, M. A stability result for the determination of order in time-fractional diffusion equations. J. Inverse -Ill-Posed Probl. 2020, 28, 379–388. [Google Scholar] [CrossRef]
  40. Sun, L.; Wei, T. Identification of the zeroth-order coefficient in a time fractional diffusion equation. Appl. Numer. Math. 2017, 111, 160–180. [Google Scholar] [CrossRef]
  41. Tikhonov, A.; Arsenin, V. Solutions of Ill-Posed Problems; V. H. Winston: Winston, WA, USA, 1977. [Google Scholar]
  42. Li, Z.; Yamamoto, M. Inverse Problems of Determining Coefficients of the Fractional Partial Differential Equations; Kochubei, A., Luchko, Y., Eds.; Volume 2 Fractional Differential Equations; De Gruyter: Berlin, Germany; Boston, MA, USA, 2019; pp. 443–464. [Google Scholar]
Figure 1. Exact (solid line) and restored (line with circles) diffusion coefficient a ( x ) (left) and solution u (right), α = 0.5 , ρ a = 0.05 , ρ u = 0.01 .
Figure 1. Exact (solid line) and restored (line with circles) diffusion coefficient a ( x ) (left) and solution u (right), α = 0.5 , ρ a = 0.05 , ρ u = 0.01 .
Fractalfract 07 00601 g001
Figure 2. Exact (solid line) and restored (line with circles) diffusion coefficient a ( x ) (left) and solution u (right), α = 0.9 , ρ a = 0.05 , ρ u = 0.01 .
Figure 2. Exact (solid line) and restored (line with circles) diffusion coefficient a ( x ) (left) and solution u (right), α = 0.9 , ρ a = 0.05 , ρ u = 0.01 .
Fractalfract 07 00601 g002
Figure 3. Exact (solid line) and restored (line with circles) diffusion coefficient a ( x ) , α = 0.25 (left) and α = 0.95 (right), ρ a = 0.1 , ρ u = 0.03 .
Figure 3. Exact (solid line) and restored (line with circles) diffusion coefficient a ( x ) , α = 0.25 (left) and α = 0.95 (right), ρ a = 0.1 , ρ u = 0.03 .
Fractalfract 07 00601 g003
Figure 4. Dependence on the fractional order (left): exact (solid line) and restored diffusion coefficient a ( x ) , α = 0.25 (line with triangles), α = 0.5 (line with squares), α = 0.95 (line with circles), ρ a = 0.1 , ρ u = 0.08 ; dependence on the perturbation (right): Exact (solid line) and restored diffusion coefficient a ( x ) , ρ a = 0.01 , ρ u = 0.01 (line with triangles), ρ a = 0.05 , ρ u = 0.01 (line with circles), ρ a = 0.01 , ρ u = 0.05 (line with squares), ρ a = 0.1 , ρ u = 0.1 (line with stars).
Figure 4. Dependence on the fractional order (left): exact (solid line) and restored diffusion coefficient a ( x ) , α = 0.25 (line with triangles), α = 0.5 (line with squares), α = 0.95 (line with circles), ρ a = 0.1 , ρ u = 0.08 ; dependence on the perturbation (right): Exact (solid line) and restored diffusion coefficient a ( x ) , ρ a = 0.01 , ρ u = 0.01 (line with triangles), ρ a = 0.05 , ρ u = 0.01 (line with circles), ρ a = 0.01 , ρ u = 0.05 (line with squares), ρ a = 0.1 , ρ u = 0.1 (line with stars).
Fractalfract 07 00601 g004
Figure 5. Error of the solution for ρ a = ρ u = 0.01 (left) and ρ a = ρ u = 0.1 (right).
Figure 5. Error of the solution for ρ a = ρ u = 0.01 (left) and ρ a = ρ u = 0.1 (right).
Fractalfract 07 00601 g005
Table 1. Errors and spatial convergence rate of the solution of the direct problem.
Table 1. Errors and spatial convergence rate of the solution of the direct problem.
α N = 40 N = 80 N = 160 N = 320 N = 640 N = 1280
0.25 E 1.536 × 10 4 3.916 × 10 5 9.886 × 10 6 2.4838 × 10 6 6.221 × 10 7 1.557 × 10 7
CR h 1.9721.9861.9931.9971.999
0.5 E 5.351 × 10 4 1.398 × 10 4 3.507 × 10 5 8.739 × 10 6 2.191 × 10 6 5.489 × 10 7
CR h 1.9361.9952.0051.9961.997
0.85 E 9.656 × 10 4 2.559 × 10 4 6.393 × 10 5 1.586 × 10 5 3.974 × 10 6 9.969 × 10 7
CR h 1.9162.0012.0121.9971.997
1 E 1.086 × 10 3 2.891 × 10 4 7.226 × 10 5 1.792 × 10 5 4.498 × 10 6 1.129 × 10 6
CR h 1.9092.0002.0111.9951.996
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Koleva, M.N.; Vulkov, L.G. Numerical Analysis of Direct and Inverse Problems for a Fractional Parabolic Integro-Differential Equation. Fractal Fract. 2023, 7, 601. https://doi.org/10.3390/fractalfract7080601

AMA Style

Koleva MN, Vulkov LG. Numerical Analysis of Direct and Inverse Problems for a Fractional Parabolic Integro-Differential Equation. Fractal and Fractional. 2023; 7(8):601. https://doi.org/10.3390/fractalfract7080601

Chicago/Turabian Style

Koleva, Miglena N., and Lubin G. Vulkov. 2023. "Numerical Analysis of Direct and Inverse Problems for a Fractional Parabolic Integro-Differential Equation" Fractal and Fractional 7, no. 8: 601. https://doi.org/10.3390/fractalfract7080601

APA Style

Koleva, M. N., & Vulkov, L. G. (2023). Numerical Analysis of Direct and Inverse Problems for a Fractional Parabolic Integro-Differential Equation. Fractal and Fractional, 7(8), 601. https://doi.org/10.3390/fractalfract7080601

Article Metrics

Back to TopTop