A Quasilinearization Approach for Identification Control Vectors in Fractional-Order Nonlinear Systems

Abstract

This paper is concerned with solving the problem of identifying the control vector problem for a fractional multi-order system of nonlinear ordinary differential equations (ODEs). We describe a quasilinearization approach, based on minimization of a quadratic functional, to compute the values of the unknown parameter vector. Numerical algorithm combining the method with appropriate fractional derivative approximation on graded mesh is applied to SIS and SEIR problems to illustrate the efficiency and accuracy. Tikhonov regularization is implemented to improve the convergence. Results from computations, both with noisy-free and noisy data, are provided and discussed. Simulations with real data are also performed.

1. Introduction

Many real processes are delayed in time, and the resulting memory effects require the use of fractional calculus to model them mathematically. For example, the problems, based on time-fractional derivative and a set of coupled nonlinear ODEs, describes the concentrations of species in so-called compartment problems [1]. The model posits that each compartment exhibits homogeneity within a well-mixed population. The coupling terms in the ODEs depict interactions among populations in distinct compartments. These terms could entail straightforward constant rate removal processes, or they might signify reactions among various populations. For certain compartmental models, understanding the time of entry of an individual into a compartment is crucial. Fractional-order systems of ODEs are one of the most effective approaches for studying transmission dynamics of epidemic [2,3,4,5,6,7,8,9,10], fractional pharmacokinetic processes [11], honeybee population dynamics [12,13], mechanics [14], etc. Recently, there has also been an increase in the number of studies based on COVID-19 fractional epidemiological modeling, see e.g., refs. [15,16,17]. One class of these models, describing the spread of an infectious disease, are the SIR (susceptible-infected-removed) compartment model, SEIR (susceptible-exposed-infected-recovered), etc.

A numerical solution of epidemic fractional SIR model of measles with and without demographic effects by using Laplace Adomian decomposition method is presented in [2]. The effect of the fractional order was investigated. The homotopy analysis method is proposed in [3] for the SIR model that tracks the temporal dynamics of childhood illness in the face of a preventive vaccine. Authors of [8] provide some remarkable and helpful properties of the SIR epidemic model with nonlinear incidence and Caputo fractional derivative.

The SEIR fractional differential model with Caputo operator is considered in [7]. Authors construct a numerical method in order to explain and comprehend influenza A(H1N1) epidemics. They conclude that the nonlinear fractional-order epidemic model is appropriate for the study of influenza and fits the real data.

In [10], the author proposes a SEIR measles epidemiological model with a Caputo fractional order operator. Existence, uniqueness and positivity of the solution was proved. The problem is solved by the Adams technique. It is shown that the incorporation of memory features into dynamical systems designed for modeling infectious diseases reveals concealed dynamics of the infection that classical derivatives are unable to detect. Existence and uniqueness results for the Caputo fractional nonlinear ODE system, describing transmission dynamics of COVID-19, are obtained in [15,17]. In [13], the author presents a existence, uniqueness and stability analysis for a fractional multi-order honeybee colony population nonlinear ODE system with Caputo and Caputo–Fabrizio operators.

One of the most popular approaches for numerically solving nonlinear ODEs and systems of ODEs with Caputo derivative is the fractional Adam method [18,19,20]. This method is applied, for example in [10] for solving the Caputo fractional SEIR measles epidemiological problem. In [17], authors a construct numerical method based on the Adams–Bashforth–Moulton approximation for solving the nonlinear ODEs system with a Caputo operator.

The Newton–Kantorovich method for solving both classical and fractional nonlinear initial-boundary values and initial value problems is employed in many papers, see e.g., refs. [21,22,23,24]. In [25], the quasi-Newton’s method, based on the Fréchet derivative, is developed for solving nonlinear fractional order differential equations.

The original concept of the quasilinearization method was developed by Bellman and Kalaba [26] and recently has been widely investigated and applied by Ben-Romdhane et.al., to several highly nonlinear problems [27] combined with other methods by Feng et al. [28] and by Sinha et al. [29]. By combining the technique of lower and upper solutions with the quasilinearization method and incorporating the Newton–Fourier idea, it becomes possible to simultaneously construct lower and upper bounding sequences that converge quadratically to the solution of the given problem [30]. This approach, referred to as generalized quasilinearization, has been further refined and theoretically extended to address a broad class of problems.

The quasilinearization method is a realization of Bellman–Kalaba quasilinearization, representing a generalization of the Newton–Raphson method. It is employed for solving either scalar or systems of nonlinear ordinary differential equations or nonlinear partial differential equations. One can also interpret quasilinearization as Newton’s method for the solution of a nonlinear differential operator equation.

In [31,32,33] existence, uniqueness and convergence results for the generalized quasilinearization method for the Caputo fractional nonlinear differential equation, represented as the Volterra fractional integral equation, are obtained.

Inverse or identification problem solutions have a key role in many practical applications [34,35,36,37]. It is a process of reconstruction from a set of observations/measurements, unknown parameters/coefficients, and initial conditions or source, i.e., the quantities that cannot be directly measured. Such problems have attracted a lot of attention from many researchers during the last decade. In spite of many results on the existence, uniqueness and stability of the solution for linear ODEs and PDEs, the nonlinearity of the differential equation renders it challenging to solve. In general, inverse problems are highly ill-posed [34,35,36,37,38], and this makes them difficult to solve, even numerically. Small fluctuations (noisy measurements) in the input data can lead to significant changes in the outcome. Therefore, to achieve meaningful results, the inversion or reconstruction process requires stabilization, commonly known as regularization.

Numerical methods for solving inverse problems for fractional-order ODEs are proposed in [14,39,40]. The shooting method is employed in [14] for numerically solving terminal value problems for Caputo fractional differential equations. In [39], the authors develop a numerical method for recovering a right-hand side function of a nonlinear fractional-order ODE, using Chelyshkov wavelets. Numerical discretizations of inverse problems for Caputo fractional ODEs under interval uncertainty are constructed in [40].

Recently, a lot of research on inverse problems concern fractional-order nonlinear ODEs systems. In [41] there is proposed a parameter and initial condition identification block pulse function method for fractional order systems. A time–domain determination inverse problem for fractional-order chaotic systems is solved in [42] by a hybrid particle swarm optimization–genetic algorithm method. An adaptive recursive least-square method is applied in [43] for parameter estimation in fractional-order chaotic systems.

In [44], the authors developed a new numerical approach, based on the Adams–Bashforth–Moulton method, for solving parameter estimation inverse problem for fractional-order nonlinear dynamical biological systems. Time fractional-order biological nonlinear ODE systems in Caputo sense with uncertain parameters are considered in [4,5]. In particular, the solution for HIV-I infection and Hepatitis E virus infection models is investigated. In [45], the author investigate a fractional-order delay nonlinear ODE system to study tuberculosis transmission. A numerical method for recovering dynamic parameters of tuberculosis is developed.

In [16], authors use the Monte Carlo method based on the genetic algorithm in order to estimate unknown parameters in fractional SIR models, describing the spread of the COVID-19 disease. The authors of [46] solve the numerical parameter identification inverse problem for the Caputo fractional SIR model, in order to estimate the economic damage of COVID-19 pandemic. They minimize a least-squares cost functional.

The parameter identification inverse problem for the nonlinear ODE system with Caputo differential operator describing honeybee colony collapse was solved numerically in [12]. The authors applied a gradient optimization method to minimize a quadratic cost functional.

The aim of this paper is to present an extension to the fractional multi-order (${\alpha }_i$) ODE system, the quasilinearization parameter identification method proposed in [47,48] for the classical ODEs system. We develop a robust and convergent numerical method of order $2-\underset{i}{max}{\alpha }_i$. The approach successfully recovers the unknown parameters even for noisy data. The Tikhonov regularization technique is applied, which increases the range of convergence with respect to the initial guess.

The remaining part of the paper is organized as follows. In the next section, we formulate direct and inverse problems on the base of finite number solution observations. In Section 3, the fractional-order quasilinearization inverse method is presented. The implementation of the method is described in Section 4. The usefulness of our approach is illustrated by numerical simulations for SIR and SEIR epidemiological models in Section 5. The paper is finalized with some conclusions.

2. Direct and Inverse Problems

First, we introduce the direct (forward) problem for the nonlinear fractional differential equation system with initial condition:

(1)$\frac{d^{\mathit{\alpha }}\mathbf{x}\left(t\right)}{d{t}^{\mathit{\alpha }}}=\mathbf{f}(t,\mathbf{x}\left(t\right);\mathbf{p}),t_0\le t\le T,\mathbf{x}\left(t_0\right)={\mathbf{x}}_0,$

where

$\mathit{\alpha }=({\alpha }_1,{\alpha }_2,\dots ,{\alpha }_I),\mathbf{x}\left(t\right)=(x^1\left(t\right),\dots ,x^I\left(t\right)),\mathbf{p}=(p^1,\dots ,p^L)\in {\mathbb{R}}^L,$

${\mathbb{R}}^L$ is a real L-dimensional Euclidean space, and $d^{\alpha }/d{t}^{\alpha }$ is the left Caputo fractional derivative of order $0<\alpha \le 1$ for arbitrary function $y\in \mathcal{A}C[t_0,T]$, i.e., y is absolutely continuous on $[t_0,T]$, and

(2)$\frac{d^{\alpha }y}{d{t}^{\alpha }}=\left\{\begin{array}{cc}\frac{1}{\mathsf{\Gamma }(1 - \alpha )}{\int }_{t_0}^t\frac{1}{{(t - \varsigma )}^{\alpha }}\frac{dy}{dt}d\varsigma \hfill & ,0<\alpha <1;\hfill \\ \frac{dy}{dt}\hfill & ,\alpha =1,\hfill \end{array}\right.$

The vector function

$\mathbf{f}(t,\mathbf{x};\mathbf{p})=\left(f^1(t,\mathbf{x};\mathbf{p}),\dots ,f^I(t,\mathbf{x};\mathbf{p})\right):[t_0,T]\times {\mathbb{R}}^I\times {\mathbb{R}}^L\to {\mathbb{R}}^I,$

where the symbol ‘;’ stands to underline the dependence on parameter $\mathbf{p}$, is called vector field. We will use also the notation ${\mathbf{x}}^{\mathit{\alpha }}=\left({\left(x^1\right)}^{{\alpha }_1},{\left(x^2\right)}^{{\alpha }_2},\dots ,{\left(x^I\right)}^{{\alpha }_I}\right)$.

The direct problem is to determine the solution $\left(x^1\left(t\right),x^2\left(t\right),\dots ,x^I\left(t\right)\right)$, when the parameters $\mathbf{p}$, right-hand side $\mathbf{f}$ and initial data ${\mathbf{x}}_0$ in (1) are given.

We denote by $C^n(t_0,T)$ the set of functions $\phi$ which are n-times differentiable on the interval $[t_0,T]$ endowed by the norm

${\parallel \phi \parallel }_{C^n}=\underset{t\in [t_0,T]}{max}\sum _{k=1}^n\left|\frac{d^k\phi }{d{t}^k}\right|,n={0,1,2,\dots ,\parallel \phi \parallel }_{C^0}={\parallel \phi \parallel }_C.$

Here, $|·|$ is the Euclidean norm when the argument is the vector and the corresponding operator norm if the argument is matrix.

Further, we will suppose that for the vector-valued function, $\mathbf{f}$ from (1) will have continuous partial differential derivatives up to second order

(3)$\begin{array}{ccc}& & \frac{\partial \mathbf{f}}{\partial \mathbf{x}}=\left\{\frac{\partial f^i}{\partial x^j}\right\},i,j=1,\dots ,I,\hfill \end{array}$

(4)$\begin{array}{ccc}& & \frac{\partial \mathbf{f}}{\partial \mathbf{p}}=\left\{\frac{\partial f^i}{\partial p^l}\right\},\frac{{\partial }^2\mathbf{f}}{\partial \mathbf{p}\partial \mathbf{x}}=\left\{\frac{{\partial }^2{\mathbf{f}}^i}{\partial {\mathbf{p}}^l\partial x^j}\right\},i,j=1,\dots ,I,l=1,\dots ,L,\hfill \end{array}$

in a bounded convex region $D\subset {\mathbb{R}}^n$ and

$\mathcal{P}=\{p\in {\mathbb{R}}^L:|p^l| <\tilde{M},l=1,\dots ,L\},\tilde{M}>0.$

We will use the following result of existence and uniqueness of the solution to the problem (1).

Theorem 1. 

Let ${\alpha }_i\in (0,1)$, $i=1,2,\dots ,I$ and the function $\mathbf{f}(t,\mathbf{x})$ satisfies the conditions (3) and (4). Then, the problem (1) has one solution $\mathbf{x}\in C^2[t_0,T]$, which is the solution of nonlinear Volterra integral equation

(5) $\mathbf{x}\left(t\right)={\mathbf{x}}_0+\frac{1}{\mathsf{\Gamma }\left(\alpha \right)}{\int }_{t_0}^t{(t-\varsigma )}^{\alpha -1}\mathbf{f}(\varsigma ,\mathbf{x}\left(\varsigma \right);\mathbf{p})d\varsigma$

and vice versa.

Proof. 

For proof, we use the results of [49], as well of [19,50].

The condition (3) is sufficient for f to be Lipschitz continuous for all $t\in [t_0,T]$ and $p\in \mathcal{P}$. Also, $\left|\frac{\partial f^i}{\partial x^j}(t,x\left(t\right);p)\right|\le \mathrm{const}.$, $i,j=1,2,\dots ,I$, $p\in \mathcal{P}$, $t\in [t_0,T]$. Then, Theorem 2 from [49] and Theorem 3.1 from [19] assure global existence of unique solution $x\in C^2[t_0,T]$ to the problem (1), as well as the representation (5).    □

Further, since the system (1) depends on the vector parameters $\mathbf{p}=(p^1,p^2,\dots ,p)$, we state the following assertion.

Corollary 1. 

Let the conditions of Theorem 1 and (4) be fulfilled. Then, the solution $x(t,\mathbf{p})$ has continuous partial derivatives with respect to $\mathbf{p}=(p^{‘}1,p^2,\dots ,p^L)$.

Proof. 

Let us consider the auxiliary problem:

(6)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{d^{\mathit{\alpha }}\mathbf{x}\left(t\right)}{d{t}^{\mathit{\alpha }}}=\mathbf{f}(t,\mathbf{x}\left(t\right);\mathbf{p}),t_0\le t\le T,\mathbf{x}\left(t\right)={\mathbf{x}}_0,\hfill \\ & \frac{d^{\mathit{\alpha }}\mathbf{p}}{d{t}^{\mathit{\alpha }}}=0,0\le t\le T.\hfill \end{array}\end{array}$

Then, the analysis of the system (6), similar to that of the Theorem 1 provides the proof.    □

The inverse problem in this paper consists of identifying the parameter vector $\mathbf{p}\in \mathcal{P}$, where $\mathcal{P}$ is an admissible set upon the observed behavior of the solution of dynamical system (1).

The idea of the parameter identification inverse problem is to minimize the difference between a measured state ${\mathbf{x}}_{\sigma }(t;\mathbf{p})$ of the system and the ones calculated with the mathematical model. However, in the real practice, the measurements are available at a final number of times. So, let us denote ${\mathbf{x}}_{\sigma }(t_m;\mathbf{p})$, $m=1,2\dots ,\mathcal{M}$ as the measured state of the dynamical system (1) and $\mathbf{x}(t_m;\mathbf{p})$, $m=1,2\dots ,\mathcal{M}$ as those obtained by solving the fractional-order differential problem (1). We gather the residuals for the measurements to be calibrated

(7)${\mathbf{r}}_m\left(\mathbf{p}\right)=\mathbf{x}(t_m;\mathbf{p})-{\mathbf{x}}_{\sigma }(t_m;\mathbf{p}),m=1,2,\dots ,\mathcal{M},$

in the residual vector $\mathbf{r}\left(\mathbf{p}\right)\in {\mathbb{R}}^{\mathcal{M}}$, i.e.

(8)$\mathbf{r}\left(\mathbf{p}\right)={[r_1\left(\mathbf{p}\right),r_2\left(\mathbf{p}\right),\dots ,r_{\mathcal{M}}\left(\mathbf{p}\right)]}^{T_r},$

where $T_r$ means the transpose operation.

We treat the inverse problem in the nonlinear least-square form

(9)$\underset{p\in \mathcal{P}}{min}F\left(\mathbf{p}\right),F\left(\mathbf{p}\right):=\frac{1}{2}{\parallel \mathbf{r}\left(\mathbf{p}\right)\parallel }^2=\frac{1}{2}{\mathbf{r}}^{T_r}\left(\mathbf{p}\right)\mathbf{r}\left(\mathbf{p}\right).$

This conception is developed in the next section, where a functional of type $\underset{t_0}{\overset{T}{\int }}F\left(\delta \mathbf{p}\right)$ with $\delta \mathbf{p}$ being residual on iterations of quazilinearization of the problem (1).

3. Quasilinearization Optimization Approach to (1)

For clarity of the exposition, we describe the method for the simple case of two equations with three unknown parameters, i.e., $I=2$, $L=3$.

Let $\left(x_0^1\left(t\right),x_0^2\left(t\right)\right)$ be an initial approximation. Applying quasilinearization [26,51] to (1), at each iteration $k=1,2,\dots ,$ we arrive at the linear Cauchy problem:

(10)$\begin{array}{ccc}& & \begin{array}{cc}\hfill \frac{d^{{\alpha }_i}{x}_k^i}{d{t}^{{\alpha }_i}}=& f^i(t,x_{k-1}^1,x_{k-1}^2;\mathbf{p})+\frac{\partial f^i}{\partial x^i}(t,x_{k-1}^1,x_{k-1}^2;\mathbf{p})(x_k^i-x_{k-1}^i)\hfill \\ & +\frac{\partial f^i}{\partial x^{3-i}}(t,x_{k-1}^1,x_{k-1}^2;\mathbf{p})(x_k^{3-i}-x_{k-1}^{3-i}),i=1,2,\hfill \end{array}\hfill \end{array}$

(11)$\begin{array}{ccc}& & x^i\left(t_0\right)=x_0^i,i=1,2.\hfill \end{array}$

Following [26], for a known parameter vector $\mathbf{p}=(p^1,p^2,p^3)$, solving the problem (10) and (11) one can construct a sequence $\left\{{\mathbf{x}}_k\left(t\right)\right\}$, which converges to $\mathbf{x}(t;\mathbf{p})$ with a quadratic rate of convergence. However, we do not know the exact value of parameter $\mathbf{p}\in \mathcal{P}$, so that starting from an an initial guess ${\mathbf{p}}_0\in \mathcal{P}$, after each quasilinearization iteration we specify the parameter value ${\mathbf{p}}_k={\mathbf{p}}_{k-1}+\delta {\mathbf{p}}_k$ from the linear system (10) and (11) by the following sensitivity scheme.

In line with the theory of sensitivity functions [52], we present the differences $\{x^i\left(t\right)-x_k^i(t;{\mathbf{p}}_k)\}$, $i=1,2$, with second-order accuracy as follows:

(12)$\begin{array}{c}\hfill r^i=x^i(t;\mathbf{p})-x_k^i(t;{\mathbf{p}}_k)=\underset{\delta x_k^i}{\underbrace{x^i(t;\mathbf{p})-[x_k^i(t;{\mathbf{p}}_{k-1})}}+{\displaystyle \sum _{l=1}^3}\frac{\partial x_k^i}{\partial p^l}(t;p_{k-1})(p_k^l-p_{k-1}^l)]\\ \hfill +o\left(\right|\delta \mathbf{p}\left|\right)as\left|\delta \mathbf{p}\right|\to 0,i=1,2.\end{array}$

The Jacobian $P_k$ with elements ${\left(P_k^{il}\right)}_{i,l=1,1}^{2,3}$, called sensitivity matrix,

(13)$P_k\equiv {\left(P_k^{il}\right)}_{i,l=1,1}^{2,3}=\left(\begin{array}{ccc}\frac{\partial x^1}{\partial p^1}(t;{\mathbf{p}}_{k-1})& \frac{\partial x^1}{\partial p^2}(t;{\mathbf{p}}_{k-1})& \frac{\partial x^1}{\partial p^3}(t,{\mathbf{p}}_{k-1})\\ \frac{\partial x^2}{\partial p^1}(t;{\mathbf{p}}_{k-1})& \frac{\partial x^2}{\partial p^2}(t;{\mathbf{p}}_{k-1})& \frac{\partial x^2}{\partial p^3}(t,{\mathbf{p}}_{k-1})\end{array}\right)=\left[\begin{array}{c}{P}_k^1,P_k^2,P_k^3\end{array}\right],$

can be found after solving the system (10) and (11). Here, $P_k^s$, $s=1,2,3$ are the corresponding columns of the matrix $P_k$.

We suppose the existence of continuous mixed derivatives $\frac{{\partial }^2{f}^i}{\partial p^l\partial x^j}$, $i=1,2$, $l=1,2,3$, $j=1,2$. Then, the sensitivity matrix $P_k$ is a solution of the following matrix equation, obtained by differentiation of the system (10) with respect to the vector $\mathbf{p}$:

(14)$\frac{d^{{\alpha }_i}{P}_k}{d{t}^{{\alpha }_i}}=\mathcal{J}(t,x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1})P_k+Q_k+\frac{\partial f(t,x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1})}{\partial \mathbf{p}},P_k\left(0\right)=0,$

where the Jacobian $\mathcal{J}(t,x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1})$ is defined as follows

$\begin{array}{ccc}& & \mathcal{J}(t,x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1})={\left[\begin{array}{c}{\mathcal{J}}^1(t,x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1}),{\mathcal{J}}^2(t,x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1})\end{array}\right]}^{T_r}={\left[\begin{array}{c}{\mathcal{J}}^1,{\mathcal{J}}^2\end{array}\right]}^{T_r},\hfill \\ & & {\mathcal{J}}^i=\left[\begin{array}{c}\frac{\partial f^i}{\partial x^i}(t,x_{k-1}^1,x_{k-1}^2;\mathbf{p})(x_k^i-x_{k-1}^i),\frac{\partial f^i}{\partial x^{3-i}}(t,x_{k-1}^1,x_{k-1}^2;\mathbf{p})(x_k^{3-i}-x_{k-1}^{3-i})\end{array}\right],\hfill \end{array}$

and $Q_k$ is $2\times 3$ matrix with elements

$q_k^{il}=\sum _{j=1}^2\frac{{\partial }^2{f}^i(t,x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1})}{\partial p^l\partial x^j}(x_k^j-x_{k-1}^j),i=1,2,l=1,2,3.$

The increment $\delta {\mathbf{p}}_k$ is obtained by minimization of the functional [47,48]:

(15)$J_k\left(\delta \mathbf{p}\right)={\int }_{t_0}^T{\left(\left(\begin{array}{c}\delta x_k^1\\ \delta x_k^2\end{array}\right)-P_k\delta \mathbf{p}\right)}^{T_r}\left(\left(\begin{array}{c}\delta x_k^1\\ \delta x_k^2\end{array}\right)-P_k\delta \mathbf{p}\right)dt.$

It is easy to see that the functional $J_k\left(\delta \mathbf{p}\right)$ is twice continuously differentiable and we have the extremum necessary condition

(16)$J_k^{\prime }\left(\delta \mathbf{p}\right)=2{\int }_{t_0}^T\left(P_k^{T_r}{P}_k\delta \mathbf{p}-P_k^{T_r}\left(\begin{array}{c}\delta x_k^1\\ \delta x_k^2\end{array}\right)\right)dt=0$

and

$J_k^{\prime \prime }\left(\delta p\right)=2{\int }_{t_0}^T{P}_k^{T_r}{P}_{k}dt>0.$

It follows from (16), that the minimum $\delta {\mathbf{p}}_k$ satisfies the system of algebraic equations

(17)$C_k\delta \mathbf{p}=D_k,C_k={\left\{c_k^{lj}\right\}}_{l=1,j=1}^{3,3}={\int }_{t_0}^T{P}_k^{T_r}{P}_{k}dt.$

Therefore, $C_k$ is $3\times 3$ symmetric matrix with elements

$c_k^{lj}={\int }_{t_0}^T{\left(\begin{array}{c}\frac{\partial x_k^1(t;{\mathbf{p}}_{k-1})}{\partial p^l}\\ \frac{\partial x_k^2(t;{\mathbf{p}}_{k-1})}{\partial p^l}\end{array}\right)}^{T_r}\left(\begin{array}{c}\frac{\partial x_k^1(t;{\mathbf{p}}_{k-1})}{\partial p^j}\\ \frac{\partial x^2(t;{\mathbf{p}}_{k-1})}{\partial p^j}\end{array}\right)dt$

and

$D_k={\left\{d_k^l\right\}}_{l=1}^3={\int }_{t_0}^T{P}_k^T\left(\begin{array}{c}\delta x_k^1\\ \delta x_k^2\end{array}\right)dt,$

is a 3-component vector:

$d_k^l=\underset{t_0}{\overset{T}{\int }}\left(\begin{array}{c}\frac{\partial x_k^1(t;{\mathbf{p}}_{k-1})}{\partial p^l}\\ \frac{\partial x_k^2(t;{\mathbf{p}}_{k-1})}{\partial p^l}\end{array}\right)\left(\begin{array}{c}{x}^1(t;\mathbf{p})-x^1(t;{\mathbf{p}}_{k-1})\\ x^2(t;\mathbf{p})-x^2(t;{\mathbf{p}}_{k-1})\end{array}\right)dt.$

The matrix $C_k$ is a Gram matrix of vectors $p_k^j$ and

$c_k^{lj}={(P_k^l,P_k^j)}_{L_2(t_0,T)}.$

It is shown in [47,48,53] that

$det\left(C_k\right)=\mathsf{\Gamma }(P_k^1,P_k^2,P_k^3)\ge 0$

and it is strongly positive, that is to say, $C_k$ is non-singular if and only if the vectors $P_k^l$, $l=1,2,3$ are linearly independent. In this case, we summarize the approach in the following steps.

Algorithm

Step 1. Choose initial guess ${\mathbf{p}}_0$, $(x_0^1\left(t\right),x_0^2\left(t\right))$ and set $k=1$.

Step 2. Solve the linear problem (10) and (11) to find $\left(x_k^1(t;{\mathbf{p}}_{k-1}),x^2(t;{\mathbf{p}}_{k-1})\right)$, and the system (14) to find $P_k$.

Step 3. Solve the linear algebraic Equations (17) to find the new parameter value ${\mathbf{p}}_k={\mathbf{p}}_{k-1}+\delta {\mathbf{p}}_k$.

Step 4. One of the expressions can be used as a criterion to stop the iteration process

$\parallel {\mathbf{p}}_k-{\mathbf{p}}_{k-1}\parallel ,J_k\left(\delta \mathbf{p}\right),$

${\int }_0^T{\left(\begin{array}{c}{x}^1(t;\mathbf{p})-x_k^1(t;\mathbf{p})\\ x^2(t;\mathbf{p})-x_k^2(t;\mathbf{p})\end{array}\right)}^{T_r}\left(\begin{array}{c}{x}_k^1(t;\mathbf{p})-x_k^1(t;{\mathbf{p}}_m)\\ x_k^2(t;\mathbf{p})-x_k^2(t;{\mathbf{p}}_m)\end{array}\right)dt,$

when it reaches a sufficiently small value. Otherwise, $k:=k+1$ and go to Step 2.

Now, following results in [47,48], more specifically Section 3 in [47], we discuss the convergence of the quasilinearization method, based on Algorithm for general case $i=1,2,\dots ,I$, $l=1,2,\dots ,L$.

We suppose that the solution of the problem (1) satisfies

${\parallel \mathbf{x}\parallel }_C\le X=\mathrm{const}.,\underset{{\parallel \mathbf{x}\parallel }_C\le X}{max}\underset{\mathbf{p}\in \mathcal{P}}{sup}max\left({\parallel \mathbf{f}\parallel }_C,{\parallel \mathcal{J}\left(\mathbf{x}\right)\parallel }_C\right)\le Y=\mathrm{const}.$

It follows from (5), (10) and (11) the integral representation

(18)${\mathbf{x}}_k\left(t\right)={\mathbf{x}}_0\left(t\right)+\frac{1}{\mathsf{\Gamma }\left(\mathit{\alpha }\right)}\underset{t_0}{\overset{T}{\int }}{(t-\varsigma )}^{1-\mathit{\alpha }}\left[f(\tau ,{\mathbf{x}}_{k-1};{\mathbf{p}}_{k-1})+\mathcal{J}(\tau ,{\mathbf{x}}_{k-1};{\mathbf{p}}_{k-1})({\mathbf{x}}_k-{\mathbf{x}}_{k-1})\right]d\varsigma ,$

which gives for $k=1$:

$\parallel {\mathbf{x}}_1{\parallel }_C \le \parallel {\mathbf{x}}_0{\parallel }_C+Y(X+1){(T-t_0)}^{\mathit{\alpha }}+Y{(T-t_0)}^{\mathit{\alpha }}{\parallel {\mathbf{x}}_1\parallel }_C.$

If we take $T-t_0$, such that

(19)${(T-t_0)}^{{\alpha }_i}\le \frac{X - \parallel {\mathbf{x}}_0{\parallel }_C}{Y(1+2X)}\mathrm{for}\mathrm{all}i=1,2,\dots ,I,$

then $\parallel {\mathbf{x}}_1{\parallel }_C \le X$. Further, when (19) holds, by induction follows that $\parallel {\mathbf{x}}_k{\parallel }_C \le X$ for $k\ge 1$.

Next, we assume that

$\underset{{\parallel \mathbf{x}\parallel }_C\le X}{max}\underset{\mathbf{p}\in \mathcal{P}}{sup}\underset{i,j,l}{max}{∥\frac{{\partial }^2{f}^i}{\partial p^l\partial x^j}∥}_C\le Y.$

Also, we suppose that $\parallel P_k{\parallel }_C \le X$, for sufficiently small $T-t_0$ and any $k\ge 1$.

Further, using the representation (18) and following the line of consideration in (Section 3 of [47]), we derive

(20)$\parallel {\mathbf{x}}_{k+1}(t;{\mathbf{p}}_k)-{\mathbf{x}}_{k+1}(t;{\mathbf{p}}_k){\parallel }_C \le \frac{{\epsilon }^{2^k}}{X_1},$

where $X_1=\underset{1\le i\le I}{max}\left\{Y{(T-t_0)}^{{\alpha }_i}/\left(1-Y{(T-t_0)}^{{\alpha }_i}\right)\right\}$ and $\epsilon =X_1{\parallel {\mathrm{x}}_1(t;{\mathbf{p}}_k)-{\mathbf{x}}_0\left(t\right)\parallel }_C$.

Additionally, as in [47], one can deduce for the functional $J_k=J_k\left(\delta {\mathbf{p}}_k\right)$, $\delta {\mathbf{p}}_k={\mathbf{p}}_k-{\mathbf{p}}_{k-1}$ of (12), the inequalities

(21)$0\le J_{k+1}\le J_k+\mathrm{const}.{\epsilon }^{2^k}.$

If we require $\epsilon <1$, then there exists a limit $J_{*}=\underset{k\to \infty }{lim}{J}_k$.

Furthermore, if the functional $J_k\left(\delta \mathbf{p}\right)$ is a strong convex function with convex constant $0.5\varrho$, see [54], then, following the results of Section 3 in [47], one can show that

(22)$\frac{1}{2}\varrho {|\delta \mathbf{p}-\delta {\mathbf{p}}_k|}^2\le J_k\left(\delta \mathbf{p}\right)-J_k\left(\delta {\mathbf{p}}_k\right)$

and

(23)$0\le \frac{1}{2}\varrho {|\mathbf{p}-{\mathbf{p}}_k|}^2\le \mathrm{const}.{\varrho }^{2^k},0\le J_k\le \mathrm{const}.{\varrho }^{2^k}.$

Therefore, the next assertion holds.

Theorem 2. 

Assume that the vector-function $\mathbf{f}$ and its second-order derivatives (3) and (4) are continuous and bounded as $(t,\mathbf{x};\mathbf{p})\in (t_0,t_1)\times D\times \mathcal{P}$. Let $x^1(t;\mathbf{p}),x^1(t;\mathbf{p}),\dots ,x^I(t;\mathbf{p})$, $p\in \mathcal{P}$ be solution to the problem (1) and assume that the conditions (20)–(22) are fulfilled. Then, the Algorithm for the general case $i=1,2,\dots ,I$, $l=1,2,\dots ,L$ is convergent with quadratic rate of convergence and the estimate (23) holds.

4. Numerical Implementation of the Method

In this section, we present the numerical approach and discuss the realization.

4.1. Numerical Discretization

For the numerical discretization of the fractional derivative, we apply $L1$ formula on non-uniform mesh [55,56]. We define the nonuniform temporal mesh ${\overline{\omega }}_{\tau }$: $t_0<t_1<\cdots <t_M=T$ with step size ${\tau }_{n+1}=t_{n+1}-t_n$, $n=0,1,\dots ,M-1$, ${\tau }_n<{\tau }_{n+1}$ and $\tau =\underset{0\le n\le M}{max}{\tau }_n$ and denote ${\left(x^i\right)}^n=x^i\left(t_n\right)$, $i=1,2$. Assume that $M=\mathcal{M}$ and we have measurements at each grid node $t_n$.

The Caputo fractional derivative of the function ${\left(x^i\right)}^{n+1}$ is approximated by

$\begin{array}{ccc}\hfill \frac{d^{{\alpha }_i}{\left(x^i\right)}^{n+1}}{d{t}^{{\alpha }_i}}& \approx & \frac{1}{\mathsf{\Gamma }(1-{\alpha }_i)}{\displaystyle \sum _{s=0}^n}\frac{{\left(x^i\right)}^{s+1} - {\left(x^i\right)}^s}{{\tau }_{s+1}}\underset{t_s}{\overset{t_{s+1}}{\int }}{(t_{n+1}-\eta )}^{-{\alpha }_i}d\eta \hfill \\ & =& {\displaystyle \sum _{s=0}^n}\left({\left(x^i\right)}^{s+1}-{\left(x^i\right)}^s\right){\rho }_{n,s}^i,i=1,2,\hfill \end{array}$

where

${\rho }_{n,s}^i=\frac{{(t_{n+1}-t_s)}^{1-{\alpha }_i}-{(t_{n+1}-t_{s+1})}^{1-{\alpha }_i}}{\mathsf{\Gamma }(2-{\alpha }_i){\tau }_{s+1}}\mathrm{and}{\rho }_{n,n}^i=\frac{{\tau }_{n+1}^{-{\alpha }_i}}{\mathsf{\Gamma }(2-{\alpha }_i)}.$

Further, we involve the notation

$G_n\left(x^i\right):=\sum _{s=1}^n\left({\rho }_{n,s}^i-{\rho }_{n,s-1}^i\right){\left(x^i\right)}^s+{\rho }_{n,0}^i{\left(x^i\right)}^0,n=0,1,\dots ,M-1,i=1,2.$

Therefore, the discretized problem (10) and (11) become

(24)$\begin{array}{ccc}& & \begin{array}{cc}\hfill {\rho }_{n,n}^i{x}_k^i& - f^i(t_{n+1},x_{k-1}^1,x_{k-1}^2;\mathbf{p})-\frac{\partial f^i}{\partial x^i}(t_{n+1},x_{k-1}^1,x_{k-1}^2;\mathbf{p})(x_k^i-x_{k-1}^i)\hfill \\ & -\frac{\partial f^i}{\partial x^{3-i}}(t_{n+1},x_{k-1}^1,x_{k-1}^2;\mathbf{p})(x_k^{3-i}-x_{k-1}^{3-i})=G_n\left(x^i\right),i=1,2,\hfill \end{array}\hfill \end{array}$

(25)$\begin{array}{ccc}& & x^i\left(t_0\right)=x_0^i,i=1,2.\hfill \end{array}$

and the corresponding to (14) discrete sensitivity problem in vector-matrix form is

(26)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & E_{\rho }{P}_k-\mathcal{J}(t_{n+1},x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1})P_k-Q_k-\frac{\partial f(t_{n+1},x_{k-1}^1,x_{k-1}^2,{\mathbf{p}}_{k-1})}{\partial \mathbf{p}}=G_n\left(P\right),\hfill \\ & P_0\left(0\right)=0,l=1,2,3,i=1,2.\hfill \end{array}\end{array}$

where $q_k^{il}$ are computed at time layer $t_{n+1}$, $E_{\rho }=\left(\begin{array}{cc}{\rho }_{n,n}^1& 0\\ 0& {\rho }_{n,n}^2\end{array}\right)$ and $G_n\left(P\right)$ is a $2\times 3$ matrix with elements

$G_n^{il}=G_n\left(P^{il}\right)=\sum _{s=1}^n\left({\rho }_{n,s}^i-{\rho }_{n,s-1}^i\right){\left(P^{il}\right)}^s+{\rho }_{n,0}^i{\left(P^{il}\right)}^0,i=1,2,l=1,2,3.$

Here $x_k^i$ and $P_k^{il}$ are solutions $x^i$ and $P^{il}$, respectively at k-th iteration of the new time layer $t_{n+1}$, i.e., $x_k^i=x_k^i\left(t_{n+1}\right)$ and $P_k^{il}=P_k^{il}\left(t_{n+1}\right)$.

The regularization technique stabilizes the solution in case of the ill-conditioning and sensitivity of the solution to noise and perturbation. If the sensitivity vectors are not linear independent, then the matrix $C_k$ may be ill-conditioned. In this case, to solve the system (17), we use Tikhonov regularization [57]. One regularization replaces the functional $J\left(\delta \mathbf{p}\right)$ with

${∥\left(\begin{array}{c}\delta x_k^1\\ \delta x_k^2\end{array}\right)-P_k\delta \mathbf{p}∥}_{L_2(t_0,T)}^2+ϵ{\parallel \delta \mathbf{p}\parallel }^2.$

This leads to the linear system

(27)$(C_k+ϵE)\delta \mathbf{p}=D_k,$

instead of (17), where E is the identity matrix, and $ϵ$ is a regularization parameter.

Another regularization approach is to replace the functional $J\left(\delta \mathbf{p}\right)$ with

${∥\left(\begin{array}{c}\delta x_k^1\\ \delta x_k^2\end{array}\right)-P_k\delta \mathbf{p}∥}_{L_2(t_0,T)}^2+ϵ{\parallel P_{k-1}+\delta \mathbf{p}-\tilde{\mathbf{p}}\parallel }^2,$

where $\tilde{\mathbf{p}}$ is a known vector to be close to the true value of the unknown parameter. It is expected that the unknown parameter is searched in a neighborhood of $\mathbf{p}$. This leads to the following linear algebraic system of equations, instead of (17):

(28)$(C_k+ϵE)=D_k+ϵ(\tilde{\mathbf{p}}-{\mathbf{p}}_{k-1}).$

To compute the $C_k$ and $D_k$, the integrals are approximated by trapezoidal rule.

4.2. Realization of the Numerical Discretizations

Now, we summarize the results above in the following Algorithms 1 and 2 for parameter identification inverse problem (PIIP).

Algorithm 1 PIIP with regularization (27)

Require: ${\mathbf{p}}_0$, $(x_0^1\left(t\right),x_0^2\left(t\right))$, $ϵ$, $tol$ Ensure: $\mathbf{p}$, $(x^1\left(t_n\right),x^2\left(t_n\right))$, $n=1,2,\dots ,M$  $k\leftarrow 1$, $Q\leftarrow tol+1$;    while $Q>tol$ do    Find $\left(x_k^1(t_n;{\mathbf{p}}_{k-1}),x_k^2(t_n;{\mathbf{p}}_{k-1})\right)$, $n=1,2,\dots ,M$, solving (24) and (25);      Determine $P_k\left(t_n\right)$, $n=1,2,\dots ,M$ from (26);      Find ${\mathbf{p}}_k={\mathbf{p}}_{k-1}+\delta {\mathbf{p}}_k$, solving (27);    $Q\leftarrow \parallel {\mathbf{p}}_k-{\mathbf{p}}_{k-1}\parallel$;     $k\leftarrow k+1$. end while

Algorithm 2 PIIP with regularization (28)

Require: ${\mathbf{p}}_0$, $(x_0^1\left(t\right),x_0^2\left(t\right))$, $ϵ$, $tol$, $\tilde{\mathbf{p}}$. Ensure: $\mathbf{p}$, $(x^1\left(t_n\right),x^2\left(t_n\right))$, $n=1,2,\dots ,M$  $k\leftarrow 1$, $Q\leftarrow tol+1$;    while $Q>tol$ do    Find $\left(x_k^1(t_n;{\mathbf{p}}_{k-1}),x_k^2(t_n;{\mathbf{p}}_{k-1})\right)$, $n=1,2,\dots ,M$, solving (24) and 25);      Determine $P_k\left(t_n\right)$, $n=1,2,\dots ,M$ from (26);      Find ${\mathbf{p}}_k={\mathbf{p}}_{k-1}+\delta {\mathbf{p}}_k$, solving (28);    $Q\leftarrow \parallel {\mathbf{p}}_k-{\mathbf{p}}_{k-1}\parallel$;  $k\leftarrow k+1$. end while

Here, $\parallel ·\parallel$ is a discrete version of ${\parallel ·\parallel }_C$ norm, namely $\parallel \nu \parallel =\underset{0\le n\le M}{max}\left|{\nu }^n\right|$.

We will clarify two points. The first one concerns the convergence. Although the quasilinearization is the second order of convergence in the sense of Theorem 2, since we apply $O\left(M^{-(2-\alpha )}\right)$ approximation of the fractional derivative, we may expect at most a $\underset{i}{min}\{2-{\alpha }_i\}=2-\underset{i}{max}\left\{{\alpha }_i\right\}$ order of convergence of Algorithms 1 and 2.

The second point is related to the regularity of the solution. If the conditions of Theorem 1 are fulfilled, then we may use uniform temporal mesh. Following the results in [55,58], we can conclude that the order of convergence of the numerical method for the direct problem is $2-\underset{i}{max}\left\{{\alpha }_i\right\}$.

Consider the case $t_0=0$. Because the fractional differential equation is a nonlocal problem and the derivative of the solution usually exhibits singularity at $t_0=0$, it is unfeasible to employ high-order numerical methods with uniform meshes. Hence, it is reasonable to employ non-uniform meshes to capture the singularity near the initial time.

Suppose that there exist a positive constant C, such that

$\left|\frac{d^s{x}^i\left(t\right)}{d{t}^s}\right|\le C(1+t^{{\alpha }_i-s}),i=1,2,s=0,1,2,t\in [0,T],$

and consider the temporal grid defined as follows [55,58]

(29)$t_n=T{\left(\frac{n}{M}\right)}^r,r\ge 1,n=0,1,\dots ,M,$

where for $r=1$, the mesh is uniform, otherwise the grid nodes are concentrated close to the origin $t_0=0$. Then, following the results in [55,58], the convergence rate of the numerical solution of the direct problem in maximal norm is

$\parallel \mathbf{x}\left(t\right)-\mathbf{x}\left(t_n\right)\parallel \le O\left(M^{-min\{2-{\alpha }_1,r{\alpha }_1\}}+M^{-min\{2-{\alpha }_2,r{\alpha }_2\}}\right),$

or in general, the case of I number of equations

$\parallel \mathbf{x}\left(t\right)-\mathbf{x}\left(t_n\right)\parallel \le O\left(\sum _{i=1}^I{M}^{-min\{2-{\alpha }_i,r{\alpha }_i\}}\right)=O\left(M^{-\underset{i}{min}\{min\{2-{\alpha }_i,r{\alpha }_i\}\}}\right).$

Therefore, in order to obtain optimal accuracy, we use graded mesh (29), taking $r=(2-\underline{\alpha })/\underline{\alpha }$, where $\underline{\alpha }=\underset{i}{min}{\alpha }_i$.

5. Numerical Simulations

In this section, we illustrate the efficiency of the developed method for recovering three and four parameters in biological ODE systems with two and four equations, respectively. We give absolute ($E_{p^l}$) and relative (${\mathcal{E}}_{p^l}$), errors of the recovered parameters $\mathbf{p}$ and relative errors (${\mathcal{E}}_{x^i}^{\infty }$, ${\mathcal{E}}_{x^i}^2$) of the solution in maximal and $L_2$ norms, respectively

$\begin{array}{cc}& {\mathcal{E}}_{\nu }^{\infty }={\mathcal{E}}_{\nu }^{\infty }\left(M\right)=\frac{\parallel {\nu }^n-{\nu }_{*}^n\parallel }{\parallel {\nu }_{*}^n\parallel },{\mathcal{E}}_{\nu }^2={\mathcal{E}}_{\nu }^2\left(M\right)=\frac{\parallel {\nu }^n-{\nu }_{*}^n{\parallel }_2}{\parallel {\nu }_{*}^n{\parallel }_2},{\parallel \nu \parallel }_2={\left({\displaystyle \sum _{n=0}^M}{\tau }_n{\left(\nu \right)}^2\right)}^{1/2},\\ & E_{p^l}=|p_k^l-p_{*}^l|,{\mathcal{E}}_{p^l}=\frac{|p_k^l-p_{*}^l|}{p_{*}^l},l=1,2,\dots ,L,\end{array}$

for different fractional orders ${\alpha }_i$. Here, $v=x^i$, $p_{*}^j$ is the exact (true) value of the parameter $p^j$, $p_k^j$ and ${\nu }^n={\left(x^i\right)}^n$ are the recovered parameter and solution at time layer $t_n$, and ${\nu }_{*}^n={\left(x^i\right)}_{*}^n$ is the reference solution of the ODEs system at time layer $t_n$, which can be the exact or numerical solution of the direct problem with exact values of the parameters, i.e., $\mathbf{p}={\mathbf{p}}_{*}$.

In addition, we perform simulations to illustrate the model behavior with respect to parameters.

We will also investigate the case of noisy measurements, generated by

(30)$x_{\sigma }^i=x^i(t;p)+{\sigma }^{x^i}\mathcal{N},$

where the measurements $x^i(t;p)$ are obtained from direct problem, ${\sigma }^{x^i}$ is the noise level, and $\mathcal{N}$ is an M-dimensional random variable with standard normal distribution.

We take $t_0=0$, and the computations are conducted mainly for the more intricate scenario of a weak singular solution in graded mesh (29).

5.1. System of Two Equations with Three Unknown Parameters

First, we consider a simplified prototype of SIS fractional order epidemic model [2,5,8], describing the dynamics of susceptible $S\left(t\right)$ and infectious $I\left(t\right)$ individuals at any time t

(31)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{d^{{\alpha }_1}S}{d{t}^{{\alpha }_1}}=\mathsf{\Lambda }-(\beta I+\mu )S,S\left(0\right)=S_0,\hfill \\ & \frac{d^{{\alpha }_2}I}{d{t}^{{\alpha }_2}}=\beta SI-\gamma I,I\left(0\right)=I_0.\hfill \end{array}\end{array}$

Typically in such models $\beta$, $\mathsf{\Lambda }$, $\mu$ and $\gamma$ are constants, representing transmission, birth, death and recovery rates, respectively.

The inverse problem is formulated for recovering the parameters $\mathbf{p}=(\beta ,\mu ,\gamma )$, and computations are performed by Algorithm 1 with $tol=$ ${10}^{-4}$. We assume that $\mathcal{M}=M$ and the measurements are performed at grid nodes.

Example 1.  (Direct problem). First, we test the accuracy of the numerical discretization (24) and (25) for solving the direct problem (31), $T=1$ in the case of weak singular solution. To this aim, we add residual functions $f_S\left(t\right)$ and $f_I\left(t\right)$ in the right-hand side of the differential equations in (31), such that the exact solution is to be $S_{*}=0.8+E_{{\alpha }_1}(-t^{{\alpha }_1})$, $I_{*}=0.2+E_{{\alpha }_1}(-t^{{\alpha }_2})$, where $E_{{\alpha }_i}$, $i=1,2$ is the Mittag–Leffler function, namely $E_{\alpha }\left(z\right)={\displaystyle \sum _{n=0}^{\infty }}\frac{z^n}{\mathsf{\Gamma }(n\alpha +1)}$.

That corresponding to the (24) and (25) discretization for the modified system (31) with an exact solution is

(32)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\rho }_{n,n}^1{S}_k+(\beta I_{k-1}+\mu )S_k+\beta I_k{S}_{k-1}=G_n\left(S\right)+\mathsf{\Lambda }+\beta I_{k-1}{S}_{k-1}+f_S\left(t_{n+1}\right),\hfill \\ & {\rho }_{n,n}^2{I}_k-\beta I_{k-1}{S}_k-\beta I_k{S}_{k-1}+\gamma I_k=G_n\left(I\right)+f_I\left(t_{n+1}\right)-\beta I_{k-1}{S}_{k-1},\hfill \end{array}\end{array}$

where $S_k\left(0\right)=S_0$ and $I_k\left(0\right)=I_0$.

We compute the solution of the direct problem from (32), initiating iteration process at each time layer with stopping criteria $max\{\parallel S_k-S_{k-1}\parallel ,\parallel I_k-I_{k-1}\parallel \}\le$ ${10}^{-6}$, for exact values of the parameters ${\mathbf{p}}_{*}=({\beta }_{*},{\mu }_{*},{\gamma }_{*})=(0.8,0.40,0.43)$ and $\mathsf{\Lambda }=0.4$ [3].

Note that in inverse problem Algorithms 1 and 2, to find the solution $\mathbf{x}(t;{\mathbf{p}}_{k-1})$, an iterative process with the same termination criteria is used to find the solution $\mathbf{x}(t;{\mathbf{p}}_{k-1})$.

In Table 1 we present errors and the corresponding orders of convergence

${\mathcal{CR}}_{\nu }^{\infty }={log}_2\frac{{\mathcal{E}}_x^{\infty }\left(M\right)}{{\mathcal{E}}_x^{\infty }\left(2M\right)},{\mathcal{CR}}_{\nu }^2={log}_2\frac{{\mathcal{E}}_x^2\left(M\right)}{{\mathcal{E}}_x^2\left(2M\right)},$

for different and equal orders of the fractional derivative. We observe that the accuracy of the solutions S is $O\left(M^{-(2-{\alpha }_1)}\right)$ and that of the solution I is $O\left(M^{-(2-{\alpha }_2)}\right)$, both in maximal and $L_2$ norms. Therefore, the convergence of the numerical solution of the system (31) is $\mathcal{E}\le O\left(M^{-(2-max\{{\alpha }_1,{\alpha }_2\})}\right)$, where $\mathcal{E}=max\{{\mathcal{E}}_S^{\infty },{\mathcal{E}}_I^{\infty }\}$ or $\mathcal{E}=max\{{\mathcal{E}}_S^2,{\mathcal{E}}_I^2\}$.

Example 2.  (Inverse problem: exact measurements). We consider an inverse problem for recovering parameters $\mathbf{p}$ in the ODE system (31) with an exact solution, modified just as in Example 1. We consider noisy-free data ($\sigma =0$), and the measurements $x(t;p)$ in (30) are equal to the numerical solution of the direct problem, computed by (24) and (25) for the true values ${\mathbf{p}}_{*}=(0.8,0.4,0.43)$. The sensitivity problem (26) is

$\begin{array}{ccc}& & S_k^{\beta }\left({\rho }_{n,n}^1+{\beta }_k{I}_{k-1}+{\mu }_k\right)+{\beta }_k{S}_{k-1}{I}_k^{\beta }=G_n\left(S^{\beta }\right)-S_{k-1}{I}_k-I_{k-1}(S_k-S_{k-1}),\hfill \\ & & S_k^{\beta }(-{\beta }_k{I}_{k-1})+\left({\rho }_{n,n}^2-{\beta }_k{S}_{k-1}+{\alpha }_k\right)I_k^{\beta }=G_n\left(I^{\beta }\right)+I_{k-1}{S}_k+S_{k-1}(I_k-I_{k-1}),\hfill \\ & & S_k^{\mu }\left({\rho }_{n,n}^1+{\beta }_k{I}_{k-1}+{\mu }_k\right)+{\beta }_k{S}_{k-1}{I}_k^{\mu }=G_n\left(S^{\mu }\right)-S_k,\hfill \\ & & S_k^{\mu }(-{\beta }_k{I}_{k-1})+\left({\rho }_{n,n}^2-{\beta }_k{S}_{k-1}+{\alpha }_k\right)I_k^{\mu }=G_n\left(I^{\mu }\right),\hfill \\ & & S_k^{\gamma }\left({\rho }_{n,n}^1+{\beta }_k{I}_{k-1}+{\mu }_k\right)+{\beta }_k{S}_{k-1}{I}_k^{\gamma }=G_n\left(S^{\gamma }\right),\hfill \\ & & S_k^{\gamma }(-{\beta }_k{I}_{k-1})+\left({\rho }_{n,n}^2-{\beta }_k{S}_{k-1}+{\alpha }_k\right)I_k^{\gamma }=G_n\left(I^{\gamma }\right)-I_k.\hfill \end{array}$

where $\frac{\partial S_k}{\partial p^l}=S_k^{p^l}$, $\frac{\partial I_k}{\partial p^l}=S_k^{p^l}$, $S_k^{p^l}\left(0\right)=0$, $I_k^{p^l}\left(0\right)=0$, $l=1,2,3$.

The regularization parameter is determined experimentally, by multiple runs, such that to avoid singularity of the matrix $C_k$ (very small $det\left(C_k\right)$) and to obtain convergent iteration process. In Table 2, we give the regularization parameter, recovered values of $\mathbf{p}$, errors of the recovered values of the parameters and solution, and number of iterations k for different initial guesses ${\mathbf{p}}_0$ and fractional orders, $M=320$, $T=1$. The reference value of the recovered solution $(S,I)$ is the exact solutions $(S_{*},I_{*})$. Results show that independently of the fractional order, the recovering of the parameters is almost exact, and the accuracy of the restored solution $(S,I)$ is similar to the accuracy of the discrete solution of the direct problem. Therefore, we conclude that the order of convergence of the numerical method for solving PIIP with noisy-free data has the same order of convergence as the numerical method of the direct problem.

For comparison, we give also the results for integer order derivative. In this case, L1 approximation is replaced by forward time stepping. We observe that the recovering of the parameters is with higher precision, and the obtained values of the solution $(S,I)$ have the same accuracy as the discrete solution of the direct problem. Note that for this example, regularization is not required.

The convergence of the method without regularization is limited in relation to the size of the convergence range for the initial guess ${\mathbf{p}}_0$. Moreover, for this test example, a smaller ${\alpha }_1$ leads to a bigger singularity of $C_k$, and convergence is achieved in more iterations. Nevertheless, applying the regularization approach, the parameters are successfully recovered. It is not a surprise that, if the initial guess ${\mathbf{p}}_0$ approaches the true values ${\mathbf{p}}_{*}$, the regularization parameter can be decreased, the iteration process requires smaller number of iterations, and the accuracy improves.

For completeness, we illustrate the performance of the numerical approach for the smooth solution in the sense of Theorem 1. We repeat the same experiment as in Table 2, but now, we take the exact solution to be $S_{*}=0.8+t^4+2t$, $I_{*}=0.2+t^5+3t$, and the mesh is uniform. The computational results are listed in Table 3. The parameter recovering is more accurate than the case of weak singular solution, and for most sets of initial guesses, the convergence is attained for a smaller number of iterations. Furthermore, in most cases, except for ${\alpha }_1=0.4$, ${\alpha }_2=0.6$ and initial guess ${\mathbf{p}}_0=(2,3,2)$, the recovering without regularization fails.

Computational results in this example showed that if the system is not ill-conditioned, and thus if the regularization is not used, the recovering of the parameters is almost exact. Otherwise, the use of regularization is strongly recommended in order to obtain a convergent iteration process, although this is at the expense of the accuracy of the restored parameters. Nevertheless, the recovery is with good enough precision such that the reconstructed solution is with optimal accuracy.

Further, all experiments are performed for the case of weak singular solution.

Example 3.  (Inverse problem: noisy measurements). The test’s problem is the same as in Example 2, but now, we add a noise, i.e., ${\sigma }^S>0$ and ${\sigma }^I>0$ in (30). We consider three levels of noise: $10\%$ (${\sigma }^S=0.065$, ${\sigma }^I=0.03$), $5\%$ (${\sigma }^S=0.032$, ${\sigma }^I=0.015$) and $1\%$ (${\sigma }^S=0.0065$, ${\sigma }^I=0.003$).

The computations are performed for ${\mathbf{p}}_0=(2,3,2)$, $M=320$ and different fractional orders. In the case of ${\alpha }_1={\alpha }_2=0.5$, we take $ϵ=0.001$, while for ${\alpha }_1=0.7$, ${\alpha }_2=0.3$, we set $ϵ=0.0001$. In Table 4 we give the numerical results. The noise level has a bigger impact on the precision of the recovered parameters in the case of derivatives with different fractional orders of large range than the case of moderate and equal values of the fractional orders, but the convergence is faster. For ${\alpha }_1={\alpha }_2=0.5$, the iteration process requires twice as many iterations than for the case of ${\alpha }_1=0.7$, ${\alpha }_2=0.3$. In both scenarios, the precision of the determined parameters is optimal, such that the solution $(S,I)$ is sufficiently accurate.

In Figure 1, we illustrate the convergence of the iteration process, namely, we represent the values of $|\delta p_k^j|$ at each iteration for ${\alpha }_1={\alpha }_2=0.5$ and ${\alpha }_1=0.7$, ${\alpha }_2=0.3$ and $10\%$ noise. For both set of parameters, we observe the fluctuation at initial iterations, and then, as was mentioned above, for ${\alpha }_1={\alpha }_2=0.5$, the convergence is attained at a small number of iterations compared to the case ${\alpha }_1=0.7$, ${\alpha }_1=0.3$. We believe that the reason is in the smaller regularization parameter and bigger grading of the mesh, which is prescribed by (29).

5.2. System of Four Equations with Four Unknown Parameters

We consider the following dynamical system SEIR model, describing measles epidemic [10].

(33)$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \frac{d^{\alpha }S}{d{t}^{\alpha }}={\mathsf{\Lambda }}^{\alpha }-({\beta }^{\alpha }I+{\mu }^{\alpha })S,S\left(0\right)=S_0,\hfill \\ & \frac{d^{\alpha }E}{d{t}^{\alpha }}={\beta }^{\alpha }SI-({\mu }^{\alpha }+{\lambda }^{\alpha })E,E\left(0\right)=E_0,\hfill \\ & \frac{d^{\alpha }I}{d{t}^{\alpha }}={\lambda }^{\alpha }E-({\mu }^{\alpha }+{\gamma }^{\alpha })I,I\left(0\right)=I_0,\hfill \\ & \frac{d^{\alpha }R}{d{t}^{\alpha }}={\gamma }^{\alpha }I-{\mu }^{\alpha }R,R\left(0\right)=R_0,\hfill \end{array}\end{array}$

where $S\left(t\right)$, $E\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$ are susceptible, exposed, infectious and the recovered individuals at time t, measured in months. The epidemiological model parameters describe recruitment rate $\mathsf{\Lambda }$, transmission coefficient rate $\beta$, natural death rate $\mu$, rate of exposure to the epidemic $\lambda$ and recovery rate $\gamma$. Let the total population be $N\left(t\right)=S\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)$. If the birth ${\mathsf{\Lambda }}^{\alpha }$ rate is equal to death rate ${\mu }^{\alpha }N$, then N is a constant, and the problem (33) reduces to a system or three equations. Namely, the function $R\left(t\right)$ is determined from algebraic relation of the total population conservation.

For the simulations, we use the data from [10] for real confirmed measles incidence cases for May–December, 2018 in Pakistan

(34)$\begin{array}{c}\hfill \begin{array}{c}\hfill \mathsf{\Lambda }=374125,\beta =4.569662\times {10}^{-11},\mu =5.25\times {10}^{-4},\lambda =2,\gamma =1.579,\\ \hfill \alpha =8.368542\times {10}^{-1},S\left(0\right)=204989885,E\left(0\right)=3548,I\left(0\right)=6567,R\left(0\right)=0,\end{array}\end{array}$

incorporated in the normalized model (33).

The inverse problem is to recover the vector $\mathbf{p}=(\beta ,\mu ,\lambda ,\gamma )$. We use the numerical solution of the direct problem (33) and (34), as a reference (exact) solution and as the measured solution $x^i(t;p)$ in (30).

The computations are performed for the system (33) and (34), involving rescalling $\tilde{S}=S/N\left(0\right)$, $\tilde{E}=E/N\left(0\right)$, $\tilde{I}=I/N\left(0\right)$, $\tilde{R}=R/N\left(0\right)$ and to reduce the computational time, we use also the transformation $\tilde{t}=t/T$.

The additional challenging in this example is that we need to recover parameters with a very different order of magnitude. We apply both Algorithms 1 and 2 for more general case $I=L=4$ and $tol=$${10}^{-5}$.

Example 4.  (Inverse problem: noisy measurements at each $t_n$). We compute the solution of the inverse problem for noisy data (30), ${\sigma }^S=102500$, ${\sigma }^E=246$, ${\sigma }^I=205$, ${\sigma }^R=1025$, measured at each grid node $t_n$ ($\mathcal{M}=M$). The final time is $T=7$ months and $M=320$. The level of the noise of $(S^{\sigma },E^{\sigma },I^{\sigma },R^{\sigma })$ is illustrated on Figure 2.

We consider two choices of the initial guess ${\mathbf{p}}_0^1=(7.42\times {10}^{-11},1.12\times {10}^{-3},2.82,2.22)$ and ${\mathbf{p}}_0^2=(4\times {10}^{-10},2.75\times {10}^{-3},4.58,3.61)$ and examine the performance of Algorithms 1 and 2.

Vector $\tilde{\mathbf{p}}$ is chosen by giving one and the same deviation of each element of $\mathbf{p}$. The values of this deviation in percentages are given in Table 5. Errors of the recovered parameters $\mathbf{p}$ and solution $(S,E,I,R)$ for different regularization are presented in Table 5, as well.

We observe that the proposed numerical identification algorithm, realized by both algorithms, recovers the unknown parameters and solution $(S,E,I,R)$ with optimal accuracy. As can be expected, the regularization (28) is successful even for initial guess that differs significantly from their true values. The iteration process converges in four to five iterations.

In Figure 3 and Figure 4, we depict an exact and recovered solution of the system (33) for initial guess ${\mathbf{p}}_0^1$ and ${\mathbf{p}}_0^2$, computed by Algorithms 1 and 2, respectively. The better fitting of the recovered solution to the exact one is attained by Algorithm 2.

Example 5.  (Inverse problem: noisy measurements at $\mathcal{M}<M$ grid nodes). We repeat the experiments from Example 4 with only difference that now, we take measurements at $\tilde{M}=65$ grid nodes. The observations are generated from numerical solution of the direct problem at $\tilde{M}$ nodes, giving the same noise levels as in Example 4. Then, we generate simulated measurements interpolating data on the time mesh. The results are listed in Table 6. We establish that the values of the recovered parameters and corresponding errors are very closed to the ones in Table 5.

6. Conclusions

In this work, we constructed a numerical quasilinearization regularization approach for solving a parameter identification inverse problem for a fractional multi-order (${\alpha }_i$) nonlinear ODE system. The uniqueness and existence of the differential direct problem was discussed. We illustrated that the order of convergence of the numerical solution of the direct problem is $2-\underset{i}{max}\left\{{\alpha }_i\right\}$. In the case of noisy-free data, the order of convergence of the numerical solution of the inverse problem is the same as that of the the numerical solution of the direct problem, namely $2-\underset{i}{max}\left\{{\alpha }_i\right\}$. The precision of the recovered parameters is more significantly affected by the noise level in cases where derivatives have varying fractional orders over a large range, compared to situations with moderate and equal fractional order values. The convergence of the iteration process achieves a moderate number of iterations. Moreover, the precision of the determined parameters is optimal such that the recovered solution is sufficiently accurate. The recovering is successful even for $10\%$ noise added to the simulated data and initial guess, not closed to the true value vector, but the use of regularization technique is essential. Otherwise, the iteration process may not converge. The regularization substantially increases the range of convergence with respect to the initial guess and reduces the number of iterations.

In our future work, we plan to develop a numerical method for a parameter identification inverse problem for a fractional multi-order nonlinear ODE system, applying a generalized quasilinearization approach by incorporating the technique of lower and upper solutions. Moreover, we intend to extended the method for a time-fractional parabolic semilinear systems and to improve the temporal order of convergence. Additionally, we intend to develop a numerical method for recovering the fractional order of the derivative for a fractional multi-order nonlinear ODEs system.

Author Contributions

Conceptualization, M.N.K. and 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.; validation, M.N.K. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The original contributions presented in the study are included in the article.

Acknowledgments

The authors are very grateful to the anonymous reviewers whose valuable comments and suggestions improved the quality of the paper.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

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.

Figures and Tables

Figure 1. [Forumla omitted. See PDF.] at each iteration, [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] for [Forumla omitted. See PDF.] (left) and [Forumla omitted. See PDF.], [Forumla omitted. See PDF.] (right), Example 3.

Figure 2. The noise ([Forumla omitted. See PDF.]), added to the numerical solution of the direct problem in order to generate [Forumla omitted. See PDF.] in (30), Example 4.

Figure 3. Exact (solid line) and recovered (line with circles) solution, Algorithm 1, [Forumla omitted. See PDF.], initial guess [Forumla omitted. See PDF.], Example 4.

Figure 4. Exact (solid line) and recovered (line with circles) solution, Algorithm 2, [Forumla omitted. See PDF.], 10% deviation of [Forumla omitted. See PDF.], initial guess [Forumla omitted. See PDF.], Example 4.

Table 1

Errors and convergence rate of the solution $(S,I)$, direct problem, Example 1.

${\mathit{\alpha }}_1$	${\mathit{\alpha }}_2$		$\mathit{M}=20$	$\mathit{M}=40$	$\mathit{M}=80$	$\mathit{M}=160$	$\mathit{M}=320$	$\mathit{M}=640$	$\mathit{M}=1280$
0.5	0.5	${\mathcal{E}}_S^{\infty }$	2.552 × ${10}^{-3}$	1.018 × ${10}^{-3}$	3.902 × ${10}^{-4}$	1.454 × ${10}^{-4}$	5.323 × ${10}^{-5}$	1.928 × ${10}^{-5}$	6.940 × ${10}^{-6}$
		${\mathcal{CR}}_S^{\infty }$		1.326	1.383	1.424	1.450	1.465	1.474
		${\mathcal{E}}_I^{\infty }$	8.805 × ${10}^{-3}$	3.067 × ${10}^{-3}$	1.078 × ${10}^{-3}$	3.809 × ${10}^{-4}$	1.349 × ${10}^{-4}$	4.782 × ${10}^{-5}$	1.695 × ${10}^{-5}$
		${\mathcal{CR}}_I^{\infty }$		1.521	1.508	1.501	1.498	1.496	1.497
		${\mathcal{E}}_S^2$	2.203 × ${10}^{-3}$	7.655 × ${10}^{-4}$	2.682 × ${10}^{-4}$	9.462 × ${10}^{-5}$	3.349 × ${10}^{-5}$	1.187 × ${10}^{-5}$	4.209 × ${10}^{-6}$
		${\mathcal{CR}}_S^2$		1.525	1.513	1.503	1.498	1.496	1.496
		${\mathcal{E}}_I^2$	1.366 × ${10}^{-2}$	4.726 × ${10}^{-3}$	1.655 × ${10}^{-3}$	5.838 × ${10}^{-4}$	2.066 × ${10}^{-4}$	7.320 × ${10}^{-5}$	2.594 × ${10}^{-5}$
		${\mathcal{CR}}_I^2$		1.531	1.514	1.503	1.499	1.497	1.497
0.7	0.3	${\mathcal{E}}_S^{\infty }$	3.191 × ${10}^{-3}$	1.418 × ${10}^{-3}$	6.138 × ${10}^{-4}$	2.608 × ${10}^{-4}$	1.096 × ${10}^{-4}$	4.571 × ${10}^{-5}$	1.877 × ${10}^{-5}$
		${\mathcal{CR}}_S^{\infty }$		1.170	1.209	1.235	1.251	1.262	1.284
		${\mathcal{E}}_I^{\infty }$	1.057 × ${10}^{-2}$	4.060 × ${10}^{-3}$	1.581 × ${10}^{-3}$	6.213 × ${10}^{-4}$	2.457 × ${10}^{-4}$	9.847 × ${10}^{-5}$	3.252 × ${10}^{-5}$
		${\mathcal{CR}}_I^{\infty }$		1.381	1.361	1.348	1.339	1.319	1.599
		${\mathcal{E}}_S^2$	3.657 × ${10}^{-3}$	1.623 × ${10}^{-3}$	7.01 × ${10}^{-4}$	2.983 × ${10}^{-4}$	1.255 × ${10}^{-4}$	5.251 × ${10}^{-5}$	2.165 × ${10}^{-5}$
		${\mathcal{CR}}_S^2$		1.172	1.210	1.233	1.249	1.257	1.278
		${\mathcal{E}}_I^2$	1.594 × ${10}^{-2}$	5.925 × ${10}^{-3}$	2.252 × ${10}^{-3}$	8.693 × ${10}^{-4}$	3.390 × ${10}^{-4}$	1.325 × ${10}^{-4}$	4.206 × ${10}^{-5}$
		${\mathcal{CR}}_I^2$		1.428	1.396	1.373	1.358	1.355	1.656

Table 2

Errors of the recovered parameters $\mathbf{p}$ and solution $(S,I)$, Example 2.

	${\mathit{\alpha }}_1={\mathit{\alpha }}_2=0.5$		${\mathit{\alpha }}_1=0.4$, ${\mathit{\alpha }}_2=0.6$		${\mathit{\alpha }}_1=0.7$, ${\mathit{\alpha }}_2=0.3$	${\mathit{\alpha }}_1={\mathit{\alpha }}_2=1$
${\mathbf{p}}_{\mathbf{0}}$	(1, 2, 1)	(2, 3, 2)	(1, 2, 1)	(2, 3, 2)	(2, 3, 2)	(3, 2, 3)
ϵ	0.0005	0.001	0.001	0.005	0.0001	0
${\beta }_k$	0.800010	0.800025	0.800034	0.800262	0.800000	0.800000
${\mu }_k$	0.399992	0.499980	0.399972	0.399791	0.400000	0.400000
${\gamma }_k$	0.430014	0.430043	0.430047	0.430354	0.430000	0.430000
$E_{\beta }$	1.0402 × ${10}^{-5}$	2.4504 × ${10}^{-5}$	3.4993 × ${10}^{-5}$	2.6165 × ${10}^{-4}$	2.9138 × ${10}^{-7}$	1.1768 × ${10}^{-14}$
$E_{\mu }$	8.2974 × ${10}^{-6}$	1.9546 × ${10}^{-5}$	2.7937 × ${10}^{-5}$	2.0888 × ${10}^{-4}$	2.3914 × ${10}^{-7}$	5.5511 × ${10}^{-17}$
$E_{\gamma }$	1.4362 × ${10}^{-5}$	3.3832 × ${10}^{-5}$	4.7305 × ${10}^{-5}$	3.5370 × ${10}^{-4}$	4.1960 × ${10}^{-7}$	1.2990 × ${10}^{-14}$
${\mathcal{E}}_{\beta }$	1.3003 × ${10}^{-5}$	3.0630 × ${10}^{-5}$	4.3741 × ${10}^{-5}$	3.2706 × ${10}^{-4}$	3.6422 × ${10}^{-7}$	1.4710 × ${10}^{-14}$
${\mathcal{E}}_{\mu }$	2.0743 × ${10}^{-5}$	4.8865 × ${10}^{-5}$	6.9842 × ${10}^{-5}$	5.2222 × ${10}^{-4}$	5.9785 × ${10}^{-7}$	1.3878 × ${10}^{-16}$
${\mathcal{E}}_{\gamma }$	3.3400 × ${10}^{-5}$	7.8679 × ${10}^{-5}$	1.1001 × ${10}^{-4}$	8.2256 × ${10}^{-4}$	9.7582 × ${10}^{-7}$	3.0208 × ${10}^{-14}$
$E_S^{\infty }$	5.3987 × ${10}^{-5}$	5.3988 × ${10}^{-5}$	1.3685 × ${10}^{-4}$	1.3687 × ${10}^{-4}$	1.4339 × ${10}^{-4}$	2.5810 × ${10}^{-4}$
$E_I^{\infty }$	2.5134 × ${10}^{-4}$	2.5133 × ${10}^{-4}$	5.4781 × ${10}^{-4}$	5.4770 × ${10}^{-4}$	4.4899 × ${10}^{-4}$	1.6266 × ${10}^{-3}$
$E_S^2$	3.3493 × ${10}^{-5}$	3.3494 × ${10}^{-5}$	8.5864 × ${10}^{-5}$	8.5882 × ${10}^{-5}$	1.2549 × ${10}^{-4}$	1.2412 × ${10}^{-4}$
$E_I^2$	2.0661 × ${10}^{-4}$	2.0661 × ${10}^{-4}$	3.5098 × ${10}^{-4}$	3.5094 × ${10}^{-4}$	3.3904 × ${10}^{-4}$	1.1390 × ${10}^{-3}$
k	11	16	14	41	8	10

Table 3

Errors of the recovered parameters $\mathbf{p}$ and smooth solution $(S,I)$, Example 2.

	${\mathit{\alpha }}_1={\mathit{\alpha }}_2=0.5$		${\mathit{\alpha }}_1=0.4$, ${\mathit{\alpha }}_2=0.6$		${\mathit{\alpha }}_1=0.7$, ${\mathit{\alpha }}_2=0.3$
${\mathbf{p}}_{\mathbf{0}}$	(1, 2, 1)	(2, 3, 2)	(1, 2, 1)	(2, 3, 2)	(2, 3, 2)
ϵ	0.0001	0.001	0	0.005	0.005
${\beta }_k$	0.800000	0.800000	0.800000	0.800002	0.800000
${\mu }_k$	0.400000	0.400002	0.400000	0.399995	0.400000
${\gamma }_k$	0.430000	0.429998	0.430000	0.430006	0.429999
$E_{\beta }$	1.7116 × ${10}^{-8}$	8.0290 × ${10}^{-7}$	2.0005 × ${10}^{-12}$	2.3862 × ${10}^{-6}$	1.3503 × ${10}^{-7}$
$E_{\mu }$	3.5581 × ${10}^{-8}$	1.6727 × ${10}^{-6}$	7.3495 × ${10}^{-12}$	5.1874 × ${10}^{-6}$	2.2071 × ${10}^{-7}$
$E_{\gamma }$	4.1176 × ${10}^{-8}$	1.9328 × ${10}^{-6}$	7.8799 × ${10}^{-12}$	5.7160 × ${10}^{-6}$	3.2656 × ${10}^{-7}$
${\mathcal{E}}_{\beta }$	2.1394 × ${10}^{-8}$	1.0036 × ${10}^{-6}$	2.5006 × ${10}^{-12}$	2.9828 × ${10}^{-6}$	1.6879 × ${10}^{-7}$
${\mathcal{E}}_{\mu }$	8.8951 × ${10}^{-8}$	4.1817 × ${10}^{-6}$	1.8374 × ${10}^{-11}$	1.2968 × ${10}^{-5}$	5.5178 × ${10}^{-7}$
${\mathcal{E}}_{\gamma }$	9.5758 × ${10}^{-8}$	4.4948 × ${10}^{-6}$	1.8325 × ${10}^{-11}$	1.3293 × ${10}^{-5}$	7.5943 × ${10}^{-7}$
$E_S^{\infty }$	1.0147 × ${10}^{-4}$	1.0147 × ${10}^{-4}$	1.5576 × ${10}^{-4}$	1.5576 × ${10}^{-4}$	2.6470 × ${10}^{-4}$
$E_I^{\infty }$	2.2631 × ${10}^{-4}$	2.2631 × ${10}^{-4}$	2.8133 × ${10}^{-4}$	2.8133 × ${10}^{-4}$	8.3472 × ${10}^{-4}$
$E_S^2$	4.9615 × ${10}^{-5}$	4.9615 × ${10}^{-5}$	8.3308 × ${10}^{-5}$	8.3308 × ${10}^{-5}$	1.0833 × ${10}^{-4}$
$E_I^2$	1.4677 × ${10}^{-4}$	1.4677 × ${10}^{-4}$	1.8367 × ${10}^{-4}$	1.8367 × ${10}^{-4}$	4.8207 × ${10}^{-4}$
k	11	14	12	13	12

Table 4

Errors of the recovered parameters $\mathbf{p}$ and solution $(S,I)$, Example 3.

	${\mathit{\alpha }}_1={\mathit{\alpha }}_2=0.5$			${\mathit{\alpha }}_1=0.7$, ${\mathit{\alpha }}_2=0.3$
Noise	1%	5%	10%	1%	5%	10%
${\beta }_k$	0.7992	0.7954	0.7900	0.7985	0.7919	0.7834
${\mu }_k$	0.4008	0.4045	0.4099	0.4011	0.4060	0.4126
${\gamma }_k$	0.4281	0.4196	0.4081	0.4270	0.4144	0.3982
$E_{\beta }$	8.1197 × ${10}^{-4}$	4.5403 × ${10}^{-3}$	9.9972 × ${10}^{-3}$	1.5330 × ${10}^{-3}$	8.1016 × ${10}^{-3}$	1.6634 × ${10}^{-2}$
$E_{\mu }$	8.5923 × ${10}^{-4}$	4.5611 × ${10}^{-3}$	9.9945 × ${10}^{-3}$	1.1390 × ${10}^{-3}$	6.0346 × ${10}^{-3}$	1.2611 × ${10}^{-2}$
$E_{\gamma }$	1.9371 × ${10}^{-3}$	1.0318 × ${10}^{-2}$	2.1993 × ${10}^{-2}$	2.9961 × ${10}^{-3}$	1.5574 × ${10}^{-2}$	3.1774 × ${10}^{-2}$
${\mathcal{E}}_{\beta }$	1.0152 × ${10}^{-3}$	5.6730 × ${10}^{-3}$	1.2469 × ${10}^{-2}$	9.1622 × ${10}^{-3}$	1.0127 × ${10}^{-2}$	3.1528 × ${10}^{-2}$
${\mathcal{E}}_{\mu }$	2.0981 × ${10}^{-3}$	1.1453 × ${10}^{-2}$	2.4858 × ${10}^{-2}$	2.8476 × ${10}^{-3}$	1.5087 × ${10}^{-2}$	2.0793 × ${10}^{-2}$
${\mathcal{E}}_{\gamma }$	4.5050 × ${10}^{-3}$	2.3994 × ${10}^{-2}$	5.0916 × ${10}^{-2}$	6.9677 × ${10}^{-3}$	3.6219 × ${10}^{-2}$	7.3892 × ${10}^{-2}$
$E_S^{\infty }$	5.4023 × ${10}^{-5}$	5.4173 × ${10}^{-5}$	5.4314 × ${10}^{-5}$	1.4354 × ${10}^{-4}$	1.4424 × ${10}^{-4}$	1.4516 × ${10}^{-4}$
$E_I^{\infty }$	2.5187 × ${10}^{-4}$	2.5406 × ${10}^{-4}$	2.5714 × ${10}^{-4}$	4.4962 × ${10}^{-4}$	4.5218 × ${10}^{-4}$	4.5540 × ${10}^{-4}$
$E_S^2$	3.3502 × ${10}^{-5}$	3.3518 × ${10}^{-5}$	3.3507 × ${10}^{-5}$	1.2562 × ${10}^{-4}$	1.2618 × ${10}^{-4}$	1.2691 × ${10}^{-4}$
$E_I^2$	2.0686 × ${10}^{-4}$	2.0787 × ${10}^{-4}$	2.0923 × ${10}^{-4}$	3.3926 × ${10}^{-4}$	2.4010 × ${10}^{-4}$	3.4111 × ${10}^{-4}$
k	17	17	17	8	8	8

Table 5

Errors of the recovered parameters $\mathbf{p}$ and solution $(S,E,I,R)$ for different regularization, $ϵ=0.001$, Example 4.

${\mathbf{p}}_0$	${\mathbf{p}}_0^1$	${\mathbf{p}}_0^2$	${\mathbf{p}}_0^2$
Algorithm	Algorithm 1	Algorithm 2	Algorithm 2
deviation $\tilde{\mathbf{p}}$		5%	10%
${\beta }_k$	7.418 × ${10}^{-11}$	4.844 × ${10}^{-11}$	5.121 × ${10}^{-11}$
${\mu }_k$	5.270 × ${10}^{-4}$	5.271 × ${10}^{-4}$	5.273 × ${10}^{-4}$
${\lambda }_k$	2.820	2.120	2.241
${\gamma }_k$	2.226	1.673	1.769
$E_{\beta }$	2.848 × ${10}^{-11}$	2.743 × ${10}^{-12}$	5.512 × ${10}^{-12}$
$E_{\mu }$	1.968 × ${10}^{-6}$	2.182 × ${10}^{-6}$	2.319 × ${10}^{-6}$
$E_{\lambda }$	8.204 × ${10}^{-1}$	1.120 × ${10}^{-1}$	2.412 × ${10}^{-1}$
$E_{\gamma }$	6.477 × ${10}^{-1}$	9.479 × ${10}^{-2}$	1.904 × ${10}^{-1}$
${\mathcal{E}}_{\beta }$	6.234 × ${10}^{-1}$	6.003 × ${10}^{-2}$	1.206 × ${10}^{-1}$
${\mathcal{E}}_{\mu }$	3.750 × ${10}^{-3}$	4.156 × ${10}^{-3}$	4.416 × ${10}^{-3}$
${\mathcal{E}}_{\lambda }$	4.102 × ${10}^{-1}$	6.003 × ${10}^{-2}$	1.206 × ${10}^{-1}$
${\mathcal{E}}_{\gamma }$	4.102 × ${10}^{-1}$	6.003 × ${10}^{-2}$	1.206 × ${10}^{-1}$
${\mathcal{E}}_S^{\infty }$	4.435 × ${10}^{-5}$	4.724 × ${10}^{-5}$	4.876 × ${10}^{-5}$
${\mathcal{E}}_E^{\infty }$	6.307 × ${10}^{-2}$	1.576 × ${10}^{-2}$	3.077 × ${10}^{-2}$
${\mathcal{E}}_I^{\infty }$	6.937 × ${10}^{-2}$	1.628 × ${10}^{-2}$	3.179 × ${10}^{-2}$
${\mathcal{E}}_R^{\infty }$	9.305 × ${10}^{-2}$	1.696 × ${10}^{-2}$	3.312 × ${10}^{-2}$
${\mathcal{E}}_S^2$	2.823 × ${10}^{-5}$	2.901 × ${10}^{-5}$	3.029 × ${10}^{-5}$
${\mathcal{E}}_E^2$	6.041 × ${10}^{-2}$	3.445 × ${10}^{-2}$	6.643 × ${10}^{-2}$
${\mathcal{E}}_I^2$	6.431 × ${10}^{-2}$	3.459 × ${10}^{-2}$	6.670 × ${10}^{-2}$
${\mathcal{E}}_R^2$	3.821 × ${10}^{-2}$	9.985 × ${10}^{-3}$	1.927 × ${10}^{-2}$

Table 6

Errors of the recovered parameters $\mathbf{p}$ and solution $(S,E,I,R)$ for different regularization, $ϵ=0.001$, Example 5.

${\mathbf{p}}_0$	${\mathbf{p}}_0^1$	${\mathbf{p}}_0^2$	${\mathbf{p}}_0^2$
Algorithm	Algorithm 1	Algorithm 2	Algorithm 2
deviation $\tilde{\mathbf{p}}$		5%	10%
${\beta }_k$	4.752 × ${10}^{-11}$	4.844 × ${10}^{-11}$	5.121 × ${10}^{-11}$
${\mu }_k$	5.281 × ${10}^{-4}$	5.278 × ${10}^{-4}$	5.278 × ${10}^{-4}$
${\lambda }_k$	2.820	2.120	2.241
${\gamma }_k$	2.227	1.674	1.769
$E_{\beta }$	1.826 × ${10}^{-12}$	2.743 × ${10}^{-12}$	5.512 × ${10}^{-12}$
$E_{\mu }$	3.104 × ${10}^{-6}$	2.183 × ${10}^{-6}$	2.827 × ${10}^{-6}$
$E_{\lambda }$	8.205 × ${10}^{-1}$	1.120 × ${10}^{-1}$	2.413 × ${10}^{-1}$
$E_{\gamma }$	6.478 × ${10}^{-1}$	9.480 × ${10}^{-2}$	1.905 × ${10}^{-1}$
${\mathcal{E}}_{\beta }$	3.995 × ${10}^{-2}$	6.003 × ${10}^{-2}$	1.206 × ${10}^{-1}$
${\mathcal{E}}_{\mu }$	5.912 × ${10}^{-3}$	5.395 × ${10}^{-3}$	5.386 × ${10}^{-3}$
${\mathcal{E}}_{\lambda }$	4.102 × ${10}^{-1}$	6.003 × ${10}^{-2}$	1.206 × ${10}^{-1}$
${\mathcal{E}}_{\gamma }$	4.102 × ${10}^{-1}$	6.004 × ${10}^{-2}$	1.206 × ${10}^{-1}$

References

1. Angstmann, C.N.; Erickson, A.M.; Henry, B.I.; McGann, A.V.; Murray, J.M.; Nichols, J.A. A general framework for fractional order compartment models. SIAM Rev.; 2021; 63, pp. 375-392. [DOI: https://dx.doi.org/10.1137/21M1398549]

2. Ahmad, A.; Farman, M.; Ahmad, M.O.; Raza, N.; Abdullah, M. Dynamical behavior of SIR epidemic model with non-integer time fractional derivatives: A mathematical analysis. Int. J. Adv. Appl. Sci.; 2018; 5, pp. 123-129. [DOI: https://dx.doi.org/10.21833/ijaas.2018.01.016]

3. Arafa, A.A.M.; Rida, S.Z.; Khalil, M. Solutions of fractional order model of childhood diseases with constant vaccination strategy. Math. Sci. Lett.; 2012; 1, pp. 17-23. [DOI: https://dx.doi.org/10.12785/msl/010103]

4. Chakraverty, S.; Jena, R.M.; Jena, S.K. Time-fractional model of HIV-I infection of CD4+ T lymphocyte cells in uncertain environment. Time-Fractional Order Biological Systems with Uncertain Parameters; Synthesis Lectures on Mathematics & Statistics Springer: Cham, Switzerland, 2020.

5. Chakraverty, S.; Jena, R.M.; Jena, S.K. Time-Fractional Order Biological Systems with Uncertain Parameters; Springer: Cham, Switzerland, 2020; 144p.

6. Demirci, E.; Unal, A.; Özalp, N. A fractional order SEIR model with density dependent death rate. Hacet. J. Math. Stat.; 2011; 40, pp. 287-295.

7. González-Parra, G.; Arenas, A.J.; Chen-Charpentier, B.M. A fractional order epidemic model for the simulation of outbreaks of influenza A (H1N1). Math. Method Appl.; 2014; 37, pp. 2218-2226. [DOI: https://dx.doi.org/10.1002/mma.2968]

8. Goufo, E.F.D.; Maritz, R.; Munganga, J. Some properties of the Kermack-McKendrick epidemic model with fractional derivative and nonlinear incidence. Adv. Differ. Equ.; 2014; 2014, 278. [DOI: https://dx.doi.org/10.1186/1687-1847-2014-278]

9. Podlubny, I. Fractional Differential Equations; Academic Press: Cambridge, MA, USA, 1999.

10. Qureshi, S. Real life application of Caputo fractional derivative for measles epidemiological autonomous dynamical system. Chaos Solitons Fractals; 2020; 134, 109744. [DOI: https://dx.doi.org/10.1016/j.chaos.2020.109744]

11. Dokoumetzidis, A.; Magin, R.; Macheras, P. Fractional kinetics in multi-compartmental systems. J. Pharmacokinet. Pharmacodyn.; 2010; 37, pp. 507-524. [DOI: https://dx.doi.org/10.1007/s10928-010-9170-4]

12. Georgiev, S.G.; Vulkov, L.G. Parameter identification approach for a fractional dynamics model of honeybee population. Lect. Notes Comput. Sci.; 2022; 13127, pp. 40-48.

13. Yildiz, T.A. A fractional dynamical model for Honeybee colony population. Int. J. Biomath.; 2018; 11, 1850063.

14. Diethelm, K.; Uhlig, F. New approach to shooting methods for terminal value problems of fractional differential equations. J. Sci. Comput.; 2023; 97, 38. [DOI: https://dx.doi.org/10.1007/s10915-023-02361-9]

15. Ahmad, S.; Ullah, A.; Al-Mdallal, Q.M.; Khan, H.; Shah, K.; Khan, A. Fractional order mathematical modeling of COVID-19 transmission. Chaos Solitons Fractals; 2020; 139, 110256. [DOI: https://dx.doi.org/10.1016/j.chaos.2020.110256] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32905156]

16. Kozioł, K.; Stanisławski, R.; Bialic, G. Fractional-order SIR epidemic model for transmission prediction of COVID-19 disease. Appl. Sci.; 2020; 10, 8316. [DOI: https://dx.doi.org/10.3390/app10238316]

17. Tuan, N.H.; Mohammadi, H.; Rezapour, S. A mathematical model for COVID-19 transmission by using the Caputo fractional derivative. Chaos Solitons Fractals; 2020; 140, 110107. [DOI: https://dx.doi.org/10.1016/j.chaos.2020.110107] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33519107]

18. Diethelm, K. The Analysis of Fractional Differential Equations: An Application-Oriented Exposition Using Differential Operators of Caputo Type; Springer: Berlin/Heidelberg, Germany, 2004.

19. Diethelm, K.; Ford, N.J. Analysis of Fractional Differential Equations: Lecture Notes in Mathematics; Springer: Berlin/Heidelberg, Germany, 2002.

20. Li, C.; Tao, C. On the fractional Adams method. Comput. Math. Appl.; 2009; 58, pp. 1573-1588. [DOI: https://dx.doi.org/10.1016/j.camwa.2009.07.050]

21. Ascher, U.M.; Petzold, L.R. Computer Methods for Ordinary Differential Equations and Differential-Algebraic Equations; SIAM: Philadelphia, PA, USA, 1997.

22. Boichuk, A.A.; Chuiko, S.M. On approximate solutions of nonlinear boundary-value problems by the Newton-Kantorovich method. J. Math. Sci.; 2021; 258, pp. 594-617. [DOI: https://dx.doi.org/10.1007/s10958-021-05569-y]

23. Filipov, S.M.; Faragó, I. Implicit Euler time discretization and FDM with Newton method in nonlinear heat transfer modeling. Int. Sci. J. Math. Model.; 2018; 3, pp. 94-98.

24. Fan, L.; Yan, Y. A high order numerical method for solving nonlinear fractional differential equation with non-uniform meshes. Numerical Ethods and Applications (NMA 2018); Nikolov, G.; Kolkovska, N.; Georgiev, K. Lecture Notes in Computer Science Springer International Publishing: Cham, Switzerland, 2018; pp. 207-215.

25. Xu, M.; Niu, J.; Lin, Y. An efficient method for fractional nonlinear differential equations by quasi-Newton’s method and simplified reproducing kernel method. Math. Methods Appl. Sci.; 2017; 41, pp. 5-14. [DOI: https://dx.doi.org/10.1002/mma.4590]

26. Bellman, R.; Kalaba, R. Quasilinearization and Nonlinear Boundary-Value Problems; Elsevier Publishing Company: New York, NY, USA, 1965.

27. Ben-Romdhane, M.; Temimi, H.; Baccouch, M. An iterative finite difference method for approximating the two-branched solution of Bratu’s problem. Appl. Numer. Math.; 2019; 139, pp. 62-76. [DOI: https://dx.doi.org/10.1016/j.apnum.2019.01.003]

28. Feng, H.; Yue, X.; Wang, X.; Zhang, Z. Decoupling and quasi-linearization methods for boundary value problems in relative orbital mechanics. Nonlinear Dyn.; 2023; 111, pp. 199-215. [DOI: https://dx.doi.org/10.1007/s11071-022-07845-1]

29. Sinha, V.K.; Maroju, P. New development of variational iteration method using quasilinearization method for solving nonlinear problems. Mathematics; 2023; 11, 935. [DOI: https://dx.doi.org/10.3390/math11040935]

30. Lakshmikantham, V.; Vatsala, A. Generalized Quasilinearization for Nonlinear Problems; Mathematics and Its Applications Kluwer Academic Publishers: Dordrecht, The Netherlands, 1998; Volume 440.

31. Almuthaybiri, S.S.; Eloe, P.W.; Neugebauer, J.T. Quasilinearization and boundary value problems at resonance for Caputo fractional differential equations. Commun. Appl. Nonlinear Anal.; 2019; 26, pp. 80-100.

32. Vasundhara Devi, J.; McRae, F.A.; Drici, Z. Generalized quasilinearization for fractional differential equations. Comput. Math. Appl.; 2010; 59, pp. 1057-1062. [DOI: https://dx.doi.org/10.1016/j.camwa.2009.05.017]

33. Vasundhara Devi, J.; Suseela, C.H. Quasilinearization for fractional differential equations. Commun. Appl. Anal.; 2008; 12, pp. 407-418.

34. Hasanov, A.H.; Romanov, V.G. Introduction to Inverse Problems for Differential Equations; 1st ed. Springer: Cham, Switzerland, 2017; 261p. [DOI: https://dx.doi.org/10.1007/978-3-319-62797-7]

35. Lesnic, D. Inverse Problems with Applications in Science and Engineering; CRC Press: Abingdon, UK, 2021; 349.

36. Prilepko, A.I.; Orlovsky, D.G.; Vasin, I.A. Methods for Solving Inverse Problems in Mathematical Physics; Marcel Dekker: New York, NY, USA, 2000.

37. Samarskii, A.A.; Vabishchevich, P.N. Numerical Methods for Solving Inverse Problems in Mathematical Physics; de Gruyter: Berlin, Germany, 2007; 438p.

38. Kabanikhin, S.I. Inverse and Ill-Posed Problems; DeGruyer: Berlin, Germany, 2011.

39. Erman, S.; Demir, A.; Ozbilge, E. Solving inverse non-linear fractional differential equations by generalized Chelyshkov wavelets. Alex. Eng. J.; 2023; 66, pp. 947-956. [DOI: https://dx.doi.org/10.1016/j.aej.2022.10.063]

40. Salahshour, S.; Ahmadian, A.; Pansera, B.A.; Ferrara, M. Uncertain inverse problem for fractional dynamical systems using perturbed collage theorem. Commun. Nonlinear Sci. Numer. Simulat.; 2021; 94, 105553. [DOI: https://dx.doi.org/10.1016/j.cnsns.2020.105553]

41. Lu, Y.; Tang, Y.; Zhang, X.; Wang, S. Parameter identification of fractional order systems with nonzero initial conditions based on block pulse functions. Measurement; 2020; 158, 107684. [DOI: https://dx.doi.org/10.1016/j.measurement.2020.107684]

42. Kosari, M.; Teshnehlab, M. Non-linear fractional-order chaotic systems identification with approximated fractional-order derivative based on a hybrid particle swarm optimization-genetic algorithm method. J. AI Data Min.; 2018; 6, pp. 365-373.

43. Behinfaraz, R.; Badamchizadeh, M.; Ghiasi, A.R. An adaptive method to parameter identification and synchronization of fractional-order chaotic systems with parameter uncertainty. Appl. Math. Model.; 2016; 40, pp. 4468-4479. [DOI: https://dx.doi.org/10.1016/j.apm.2015.11.033]

44. Liu, F.; Burrage, K. Novel techniques in parameter estimation for fractional dynamical models arising from biological systems. Comput. Math. Appl.; 2011; 62, pp. 822-833. [DOI: https://dx.doi.org/10.1016/j.camwa.2011.03.002]

45. Georgiev, S. Mathematical identification analysis of a fractional-order delayed model for tuberculosis. Fractal Fract.; 2023; 7, 538. [DOI: https://dx.doi.org/10.3390/fractalfract7070538]

46. Georgiev, S.; Vulkov, L. Numerical coefficient reconstruction of time-depending integral-and fractional-order SIR models for economic analysis of Covid-19. Mathematics; 2022; 10, 4247. [DOI: https://dx.doi.org/10.3390/math10224247]

47. Abdullaev, U.G. Quasilinearization and inverse problems of nonlinear dynamics. J. Optim. Appl.; 1995; 85, pp. 509-523. [DOI: https://dx.doi.org/10.1007/BF02193053]

48. Abdulla, U.G.; Poteau, R. Identification of parameters in systems—Biology. Math. Biosc.; 2018; 305, pp. 133-145. [DOI: https://dx.doi.org/10.1016/j.mbs.2018.09.004]

49. Kilbas, A.A.; Marzan, S.A. Nonlinear differential equations with the Caputo fractional derivative in the space of continuously differentiable functions. Diff. Equ.; 2005; 41, pp. 84-89. [DOI: https://dx.doi.org/10.1007/s10625-005-0137-y]

50. Lin, W. Global existence theory and chaos control of fractional differential equations. J. Math. Anal. Appl.; 2007; 332, pp. 709-726. [DOI: https://dx.doi.org/10.1016/j.jmaa.2006.10.040]

51. Koleva, M.N.; Vulkov, L.G. Two-grid quasilinearization approach to ODEs with applications to model problems in physics and mechanics. Comput. Phys. Commun.; 2010; 181, pp. 663-670. [DOI: https://dx.doi.org/10.1016/j.cpc.2009.11.015]

52. Tomorich, R. Sansitivity Analysis of Dynamical Systems; Academic Press: New York, NY, USA, 1963.

53. Beckenbach, E.F.; Bellman, R.E. Inequalities; Springer: Berlin, Germany, New York, NY, USA, 1962.

54. Rockafellar, R.T. Convex Analysis; Princeton University Press: Princeton, NJ, USA, 1969.

55. Stynes, M.; O’Riordan, E.; Gracia, J.L. Error analysis of a finite difference method on graded meshes for a time-fractional diffusion equation. J. Numer. Anal.; 2017; 55, pp. 1057-1079. [DOI: https://dx.doi.org/10.1137/16M1082329]

56. Zhang, Y.; Sun, Z.; Liao, H. Finite difference methods for the time fractional diffusion equationo non-uniform meshes. J. Comput. Phys.; 2014; 265, pp. 195-210. [DOI: https://dx.doi.org/10.1016/j.jcp.2014.02.008]

57. Tikhonov, A.; Arsenin, V. Solutions of Ill-Posed Problems; John Wiley & Sons: New York, NY, USA, 1977.

58. Zhou, Y.; Stynes, M. Optimal convergence rates in time-fractional discretisations: The L1, ${\mathrm{L}1}^{¯}$ and Alikhanov schemes. East Asian J. Appl. Math.; 2022; 12, pp. 503-520.