On oscillatory second order impulsive neutral difference equations

1.   Introduction

The inverse problem for differential equations is part of the fascinating branches of mathematics. It is developed in connection with diverse branches of applied science ^[1]. The type of boundary condition determines the identity of the problem. If the boundary is given for an initial time the problem is recognized as an initial value problem (IVP); if the boundary is described for the final, it is called a terminal value problem (TVP). If the boundary is described at both times, it is a Sturm-Liouville problem ^[2].

For ordinary differential equations (ODEs), the terminal value problem can be transformed into an initial value problem, and they are generally well-posed problems. Surprisingly, this is not true for fractional differential equations (FDEs). Not only can a TVP for FDEs not be transformed into an IVP, but also such problems are not well-posed in general ^[3,4]. It has been recently noticed that a TVP for FDEs is well-posed in a small neighborhood of the boundary ^[5,6,7,8,9].

In recent decades, literature on fractional operators has been extensively increased for describing non-local dynamical processes ^[10]. The non-local property of such processes is usually captured by an integral with a memory kernel. Recently, in the notable paper ^[11], high dimensional problems are reduced to an integral equation with a memory kernel. This is a green light for replacing ordinary derivatives with fractional derivatives in complex processes. Thus, fractional differential equations can describe abundant models well.

There is no unified definition for a fractional derivative. However, if we restrict some properties (such as being singular or non-singular, local or non-local) we can obtain some completely separate classifications. Among them, the Caputo derivative has received more interest and attention in the literature, mainly for well-interpreted boundary conditions and more similarity to the ordinary derivative (the Caputo derivative of the constant function is zero) ^[12].

Similarly, fractional integrals have diverse definitions. The route toward the generalized fractional operator will be clear if we monitor some selected classes of fractional integrals. The Riemann-Liouville fractional integral is a generalization of the integer order integral operator

and the Hadamard fractional integral is a generalization of the $n$th order integral operator

where $n\in\mathbb{N}$ ^[13]. Katugampola ^[14,15] introduced his fractional derivative by replacing $t_i$ with $t_i^\rho,$ $\rho > 0$. Fu et al. have replaced exponential functions to introduce exponential fractional derivatives ^[16]. It is tempting to replace $t_i$ with a general function $g(t_i),$ (where $g$ is a weighted function) to study generalized fractional derivatives. Such unified generalization of fractional operators with respect to another function has been studied in ^[17,18].

It is a common aspect of pure mathematics to see an application of a seemingly not applicable concept later. Analogously, for generalized fractional derivatives, we see some noticeable applications that have been published just very recently. For example, the subdiffusion equation with a generalized fractional derivative has been effectively invoked to describe subdiffusion in a medium with an evolving structure over time ^[19]. The application of generalized fractional operators for the Fokker-Planck equation has been noticed in ^[16,20]. We will show another importance of such generalization in inverse problems.

Considering the generality of the fractional derivative with respect to another derivative and its application, it is highly motivated to study fractional differential equations with respect to another derivative.

A $\nu$-dimensional terminal value problem for a system of generalized fractional differential equations (GFDEs) is described by

and

where ${^{G}_{a}D}^{\bar{\alpha}, \mathbf{g}}_{t}$ is a generalized fractional derivative with respect to the vector function $\mathbf{g} = [g_1, \ldots, g_\nu]^T$ ($g_i :[a, b]\rightarrow \mathbb{R}$ are strictly increasing functions with continuous derivatives $g_i'$ on $(a, b)$), $\bar{\alpha} = [\alpha_1, \ldots, \alpha_\nu]^T$ ($\alpha_i \in(0, 1)$), $\mathbf{y}:[a, b]\rightarrow \mathbb{R}^n$ is an unknown vector function, $\mathbf{f}:[a, b]\times \mathbb{R}^n \rightarrow \mathbb{R}^n$ is a given vector function describing an evolutionary processes, and $\mathbf{y}_b\in \mathbb{R}^n$ is a terminal value.

In a real dynamical process, we do not expect each separate component to have the same memory. Therefore, the study of FDEs with vector order is more important than that for fixed order. We note that the vector order introduced in ^[21] is an element-wise operator and defined by

This paper contributes the following achievements:

$^*$ Study of the TVP for systems of high dimensional nonlinear generalized FDEs,

$^*$ Computing generalizing fractional operators by classical fractional operators,

$^*$ Introducing suitable weighted space that relates previous studies to generalized fractional derivative,

$^*$ Introducing and analyzing a comprehensive collocation method that covers all generalized operators,

$^*$ Obtaining a lower bound for well-posedness of TVPs with various weight functions,

$^*$ Applying generalized derivatives for estimating the past population of diabetes cases.

To the best of our knowledge, these are investigated for the first time. The important achievement of this paper is that it gives us more analyzed options in modeling dynamical processes, to choose a better weight function. If we need an inverse problem with a longer interval, the result of this paper introduces a generalized derivative that guarantees the well-posedness in a study model.

A generalized fractional derivative can be transformed into a classical fractional derivative ^[13,22]. This fact was noticed and employed for solving generalized fractional differential equations, especially in ^[22]. In this respect, we introduced a weighted space with respect to another function to carry out our analysis from generalized fractional-order derivatives toward classical fractional derivatives.

Numerical methods in parallel to theoretical analysis for solving fractional differential equations have developed rapidly. The Jacobi spectral methods for solving related functional integral and fractional differential equations have been extensively studied in ^{[23,24,25,26]}. Applying this method for terminal value problems can be found in ^[27]. The collocation methods on piecewise polynomials spaces for fractional differential equations were studied in ^[6,9]. The spectral methods provide fast convergence methods that transform a related problem into a high-dimensional algebraic equation, while the collocation method provides a high order method with more options for control of the error. As far as we know, the collocation method for terminal value problems with generalized fractional differential equations is not yet studied and not analyzed. An important aim of this paper also is to propose and analyze such methods for these problems.

This paper is organized as follows. In Section 2, we review the generalized fractional operators (integral and derivative operators). In Section 3, we transform TVPs for systems of FDEs into weakly singular Fredholm-Volterra integral equations with vanishing delay. In Section 4, we obtain a lower bound for the well-posedness of such problems. In Section 5, we propose a numerical method, and in Section 6 we provide an error analysis for the proposed method. In Section 7, after validating the proposed method, we compare the effects of various choices of weight function on modeling with a TVP.

2.   Preliminaries

The generalized fractional operator is defined with respect to the weight function $g:[a, b]\rightarrow \mathbb{R}$ with the following properties:

$C_1:$ $g\in C^1[a, b];$

$C_2:$ $g$ is a strictly increasing function;

$C_3:$ $g^{-1}:g([a, b])\rightarrow [a, b]$ exists and is continuous (by $C_2,$ $g([a, b]) = [g(a), g(b)]$).

Definition 2.1. ^[18] Let $f\in C[a, b]$ ($a, b\in \mathbb{R}$). The generalized fractional integral ${^{G}_{a}I}^{\alpha, g}_{t}$ of order $\alpha$ ($\alpha\in \mathbb{C},$ $Re(\alpha) > 0$) with the weight function $g$ is defined by

Remark 1. Particular choices of $g$ result famous generalizations and classical definitions of fractional operators:

$\bf(I).$ If $g(x) = x,$ then ${^{G}_{a}I}^{\alpha, g}_{t}$ is the classical Riemann-Liouville fractional integral,

$\bf(II).$ If $g(x) = \ln(x)$ on $[1, b],$ ($b > 1$), then ${^{G}_{a}I}^{\alpha, g}_{t}$ is the Hadamard fractional integral, in this case $g^{-1}(x) = e^x,$

$\bf(III).$ If $g(x) = \frac{e^{\lambda x}}{\lambda}$ ($\lambda > 0$) on $[a, b],$ then ${^{G}_{a}I}^{\alpha, g}_{t}$ is the exponential fractional integral, in this case $g^{-1}(x) = \frac{\ln(\lambda x)}{\lambda},$ ^[16],

$\bf(IV).$ If $g(x) = x^\rho,$ then ${^{G}_{a}I}^{\alpha, g}_{t}$ is the Katugampola fractional integral, in this case $g^{-1}(x) = x^{\frac{1}{\rho}}$ ^[14,15].

Definition 2.2. Let $X[a, b]$ be a space of real valued functions on $[a, b]$. The weighted space with respect to the weight function $g$ is defined by

with the weighted norm

where $\|.\|$ is a norm of the space $X[g(a), g(b)].$

For example, $L_g[a, b]$ and $C_g[a, b]$ are weighted spaces of the Lebesgue integrable functions and the continuous functions with Lebesgue norm and supremum norm, respectively.

Remark 2. Suppose $k\in C[a, b].$ Since $g^{-1}$ is continuous, $kog^{-1}$ is continuous. Thus, $k\in C_g[a, b]$ and $C[a, b]\subset C_g[a, b].$ Similarly, we can prove that $C_g[a, b]\subset C[a, b].$ For supremum norm in this space, we have

Thus, $C_g[a, b] = C[a, b].$ However, this type of equivalency can not happen for all spaces like $L[a, b].$ In this case, by substituting $x = g(y)$ and $dx = g'(y)dy,$ we have

Throughout the paper we generally use the supremum norm unless we mention it. The following theorem describes the generalized fractional integral by the Riemann-Liouville fractional integral.

Theorem 2.1. Suppose $\alpha\in \mathbb{C},$ $Re(\alpha) > 0,$ and $g:[a, b]\rightarrow \mathbb{R}$ has the properties $C_1$–$C_3$. Then, for $h\in C[a, b],$

where ${^{RL}_{g(a)}I}^{\alpha}_{t}$ stands for the Riemann-Liouville fractional integral, and the notation "$o$" stands for the combination operator.

Proof. Setting $k(t) = h(g(t)),$ we get

Here, we integrated with substitution $x = g(\tau),$ ($dx = g'(\tau)d\tau$). Thus, we obtain

This completes the proof.

Remark 3. Theorem 2.1 shows that ${^{RL}_{g(a)}I^{\alpha}_{t}} k(g^{-1}(t))$ is well-defined if and only if ${^{G}_{a}I}^{\alpha, g}_{t}k(t)$ is well-defined. Let $X[g(a), g(b)]$ be the space that the operator ${^{RL}_{g(a)}I^{\alpha}_{t}}$ is well-defined. Then, the operator ${^{G}_{a}I}^{\alpha, g}_{t}k(t)$ is well-defined on the weighted space $X_g[a, b].$

Definition 2.3. ^[18] Let $\alpha\in \mathbb{C},$ $Re(\alpha) > 0,$ and $n = [Re(\alpha)]+1.$ The generalized fractional derivatives on $[a, b]$ ($0\leq a < b$) are defined by

for $\alpha\notin \mathbb{N}$ and for a Borel function $f:[a, b]\rightarrow \mathbb{R}$ in which the integral in Eq (2.4) is well-defined. For $\alpha\in \mathbb{N},$ it is the classical integer order definition, i.e., ${^{K}_{a}D}^{\alpha, \rho}_{t}f(t) = \frac{d^{n-1}}{dt^{n-1}}$.

Theorem 2.2. Let $\alpha\in \mathbb{C},$ and $0 < Re(\alpha) < 1.$ Then,

Proof. Setting $h(t) = f(g(t)),$ we get

This completes the proof.

Computations of generalized derivatives and integrals of functions $(g(t)-g(a))^{v},$ for $0 < Re(\alpha) < 1$ are a direct consequence of Theorems 2.1 and 2.2. Thus, we have

and

Now, we can use Theorem 2.2 to find a suitable space for defining a generalized fractional derivative. Regularly, the spaces $AC^{n}[a, b]$ and $C^n[a, b]$ are the best choices for classical fractional derivatives with order $\alpha$ such that $n-1 < Re(\alpha) < n,$ (Theorems 2.1 and 2.2 of ^[28], also see ^[13]). We recall that

where $AC[a, b]$ is the set of absolute continuous functions. Thus, we can state the following lemma.

Lemma 2.1. Let $f\in AC^n_g[a, b].$ Then, ${^{G}_{a}D}^{\alpha, g}_{t}f(t)$ exists almost everywhere on [a, b]. Moreover, if $f\in C^n_g[a, b],$ then ${^{G}_{a}D}^{\alpha, g}_{t}f\in C[a, b].$

Proof. The proof is a straightforward result of Eq (2.5), and Theorems 2.1 and 2.2 of ^[28].

Remark 4. Equations (2.2) and (2.5) can be rewritten as

and

Lemma 2.2. Let $\alpha\in(0, 1],$ $\rho\in [0, 1),$ and $f\in AC[a, b].$Then,

Proof. It follows from Eq (2.9) that

Applying Eq (2.10), we get

Taking into account that $(g(t))o(g^{-1}(t)) = t,$ we obtain

Finally, from properties of classical fractional operators, we obtain

which completes the proof.

Remark 5. All the mentioned results can be generalized to higher-dimensional spaces by the vector order notation. For vector order $\bar{\alpha} = [\alpha_1, \ldots, \alpha_\nu]^{T}$ ($Re(\alpha_i)\in (0, 1),$ $i = 1, \ldots, \nu$) and vector weights $\mathbf{g} = [g_1, \ldots, g_\nu]^{T}$ we have

where $\mathbf{f} = [f_1, \ldots, f_\nu]^T$ belongs to the $\nu$-dimensional space $(AC[a, b])^{\nu}$. It is a straightforward generalization of a one-dimensional case since all operators operate element-wise ^[21].

As an application of Eqs (2.7), (2.8)and (2.11), we can solve a linear initial value problem for FDEs:

with the initial value $\mathbf{y}(0) = \mathbf{y}_0$ and imposing $\alpha_i = \alpha,$ $g_i = g$ for $i = 1, \cdots, \nu$. Applying ${^{G}_{a}I_t}^{\alpha, g}$ to both sides of Eq (2.14) and using Eq (2.11), we have

Here, it is worth to mentioning that, for general $\bar{\alpha}$ and $\mathbf{g},$ the commutative property

does not hold for vector order fractional integrals as well as for vector order fractional derivatives.

The iterative method

with the initial value

is utilized to obtain the exact solution $\mathbf{y} = \lim_{n\rightarrow \infty} \mathbf{y}_n$. Recursively, we obtain

and similarly

This pattern suggests

where the dot ("$.$'') operator stands for element-wise multiplication.

3.   Transforming inverse problem into integral equations

Applying the fractional integral to both sides of the system (1.2), it follows from Eq (2.2) that

or equivalently,

Computationally, it is important to note that

Putting $t = b$ in Eq (3.2) and using the terminal condition, we obtain

Substituting $\mathbf{y}(a)$ from Eq (3.3) into Eq (3.2), we obtain

Equation (3.4) can be represented by a weakly singular Fredholm-Volterra integral equation with vanishing delay $\mathbf{g}(t)$

or equivalently,

Remark 6. If $g(t) = pt,$ the corresponding system has proportional delay ^[29]. Also, with $g(t) = t^{\frac{3}{2}}$ and the non-local boundary condition

the system (3.6) is a Lighthill system. This system is used for describing the temperature distribution of the surface of a projectile moving through a laminar layer ^[30].

4.   Local existence results

4.1.   Continuity of the operator ${^{G}_{a}I_t}^{\bar{\alpha}, \mathbf{g}}\mathbf{f}(t, \mathbf{y}(t))$ with respect to $\mathbf{y}(t)$

We suppose the vector function $\mathbf{f} = [f_1, \ldots, f_\nu]^T$ satisfies the following conditions.

$\bf(H1)$ The vector function $\mathbf{f}$ is continuous with respect to its variables on $[0, T]\times \mathbb{R}^\nu,$

$\bf(H2)$ Functions $f_i,$ $i = 1, \ldots, \nu,$ are Lipschitz functions with Lipschitz constants $L_i,$ i.e.,

for $\mathbf{y}, \mathbf{w}\in \mathbb{R}^\nu,$ and we set

It is a straightforward conclusion of Theorem 2.1 and Remark 2 that the operator

transforms a continuous vector function into a continuous vector function, and thus it is well-defined.

We use the norm

to state continuity of the operator ${^{G}_{a}I_t}^{\bar{\alpha}, \mathbf{g}}\mathbf{f}(t, .).$

Theorem 4.1. Let assumptions (H1) and (H2) hold for the vector function $\mathbf{f}.$ Then, the operators ${^{G}_{a}I_t}^{\bar{\alpha}, \mathbf{g}}\mathbf{f}(t, .)$ and ${^{G}_{a}I_b}^{\bar{\alpha}, \mathbf{g}}\mathbf{f}(b, .)$ are continuous, i.e.,

and

where the vector functions $\mathbf{u}$ and $\mathbf{v}$ belong to $(C[a, b])^{\nu}$ and

Proof. The supremum norm of the $i$th component of the generalized fractional integral operator can be estimated by

By assumptions $C_1$–$C_3,$ $g_i(x)$ and $(g_i(t)-g_i(x))^{1-\alpha_i}$ for $x\in[0, t]$ are non-negative. Thus,

Consequently, the inequalities (4.1) and (4.2) follow from Eq (4.4).

4.2.   Existence

By Theorem 4.1, the operator $\mathcal{T}:(C[a, b])^{\nu}\rightarrow (C[a, b])^{\nu}$ described by

is well-defined.

Theorem 4.2. Let assumptions (H1) and (H2) hold for $\mathbf{f}: [0, T]\times \mathbb{R}^\nu\rightarrow \mathbb{R}^\nu.$Let

Then, the system (3.6) has a unique solution on $(C[a, b])^{\nu}.$

Proof. Systems (3.5) and (3.6) can be rewritten in the operator form

By using Theorem 4.1 for the vector functions $\mathbf{u}$ and $\mathbf{v}\in(C[a, b])^{\nu},$ we have

From condition (4.6), $M < 1,$ and hence

Thus, $\mathcal{T}$ is a contracting operator, and by the Banach fixed point theorem, it has a unique solution on $(C[a, b])^{\nu}.$ This completes the proof.

5.   Collocation method on piecewise polynomial spaces

Let us adopt notations of ^[29]. Let $I_h = \{t_n: \ a = t_0 < t_1 < \ldots < t_N = b\}$ be a partition of $[a, b],$ $h_n = t_{n+1}-t_n,$ $\sigma_n = (t_n, t_{n+1}]$ ($n = 0, \ldots, N-1$) and $h = \max_{n = 0, \ldots, N-1}h_n.$ We recall that the grading mesh points of the form $t_n: = (b-a)(n/N)^r$ are an appropriate choice for integral equations with weakly singular kernels^[29]. The piecewise polynomial space can be defined by

where $\Pi_{m-1}$ is the space of polynomials of degree less than $m.$ An approximate solution of the system (1.2) has the dense representation

where $L_i$ are Lagrange interpolation formulas, and $U_{n, i}: = \mathbf{u}_N(t_{n, i})$ and collocation points are defined by $t_{n, i}: = t_n+c_ih_n$ for given collocation parameters $0\leq c_1 < \ldots < c_m\leq 1.$ The operator $P_N$ that projects $\Big(C[0, T]\Big)^\nu$ into $\Big(\mathbb{S}_{m-1}^{-1}(I_h)\Big)^\nu$ is described by

Thus, we can explicitly and uniquely define $P_N$ by

The quadrature approximation of the operator ${^{G}_{a}I_t}^{\bar{\alpha}, \mathbf{g}}\mathbf{f}(t, .)$ is computed by

Let $t: = t_n+vh_n \in [t_n, t_{n+1}].$ Then,

where the weights of the quadrature are determined by

for $j = 1, \ldots, m$ and $n = 0, \ldots, N-1.$ Similarly, a quadrature method for ${^{G}_{a}I_b}^{\bar{\alpha}, \mathbf{g}}\mathbf{f}(t, .)$ can be computed by

or, equivalently,

Setting $t = t_{o, k}$ for $o = 0, \ldots, N-1$ and $k = 1, \ldots, m,$ and using the introduced quadrature method, the unknown $U_{o, k}$ can be obtained by solving the system of nonlinear equations described by

After solving this system, we can obtain the dense approximate solution by Eq (5.1). The nonlinear system (5.6) is solved by the recursive formula

where $n_t$ is the constant number of iterations (can be chosen by a user), or can be adopted by

with a given tolerance $\mbox{Tol}.$

In this research, we also need to introduce a numerical method for solving related initial value problems. Let us consider the system (1.2) with a given initial condition

A numerical approach by recursively solving the $\nu m$ dimensional algebraic systems

for $o = 0, \ldots, N-1$ is utilized to obtain the corresponding dense solution. However, for the IVP in each interval, we need to solve only a nonlinear equation of dimension $\nu\times m,$ while for the TVP we need to solve a nonlinear equation of dimension $\nu\times m\times N.$

5.1.   Existence of the approximate solution

Does the nonlinear system (5.6) have a solution? Is it unique? What is the convergence of the proposed iterative method (5.7)? This subsection answers these questions.

Theorem 5.1. Let assumptions (H1) and (H2) hold for the vector function $\mathbf{f}: [0, T]\times \mathbb{R}^\nu\rightarrow \mathbb{R}^\nu.$ Let

for a given $N.$ Then, the approximate solution of the system (5.6) exists and is unique, and the iterative method (5.7) converges to the solution of the system (5.6) on collocation points.

Proof. Since $\mathbf{u}_N\in \Big(\mathbb{S}_{m-1}^{-1}(I_h)\Big)^\nu$ we have $P_N\mathbf{u}_N = \mathbf{u}_N,$ and the dense approximate solution (5.1) satisfies

where

Let $\mathbf{u}, \mathbf{v}\in \Big(\mathbb{S}_{m-1}^{-1}(I_h)\Big)^\nu.$ The $i$th component of the $Q_N\mathbf{u}(t)-Q_N\mathbf{v}(t)$ satisfies

Therefore,

where $M$ is defined as in Eq (4.3). Similarly,

The operator $P_N\mathcal{T}_N$ is a contractor since

by the hypothesis of the theorem. It follows from the Banach fixed point theorem that the operator

has a unique solution, and the recursive formula

converges to the fixed point of $P_N\mathcal{T}_N$ with any initial solution

We recall that $\Big(\mathbb{S}_{m-1}^{-1}(I_h)\Big)^\nu$ is a finite dimensional space and thus a Banach space.

6.   Analysis of the method

In this section, we show that the error of the discretized collocation method is proportional to the error of the quadrature method for the exact solution:

Theorem 6.1. Let assumptions (H1) and (H2) hold for $\mathbf{f}: [0, T]\times \mathbb{R}^\nu\rightarrow \mathbb{R}^\nu.$ Let $\Lambda$ defined by (5.10) satisfy

Then, the approximate solution described by (5.1) and (5.6) exists, and

Proof. Since $\mathbf{u}_N\in \Big(\mathbb{S}_{m-1}^{-1}(I_h)\Big)^\nu,$ it follows that $P_N\mathbf{u}_N = \mathbf{u}_N.$ Hence, the dense approximate solution (5.1) satisfies

where

By using Eq (6.1),

Setting

and subtracting Eq (6.5) from Eq (5.10), we obtain

The first difference of Eq (6.6) satisfies

Let us consider the the $i$th component of the $Q_N\mathbf{u}_N(t)-Q_N\mathbf{y}(t),$ and we have

It follows from Eq (6.8) that

Similarly, we can prove that

Substituting Eqs (6.9) and (6.10) into Eq (6.7), we obtain

Taking the norm from Eq (6.6) and using Eq (6.11), we obtain

This leads to an interesting bound for the error,

The following lemma supports Theorem 6.1 by a bound for $\|P_N\|$.

Lemma 6.1. Let $P_N:\left(\mathcal{C}[0, T]\right)^{\nu}\rightarrow \left(S_{m-1}^{-1}(I_h) \right)^{\nu}$ be a projection defined by Eq (5.2). Let $\|.\|$ be the induced norm of projection. Then, there exist $N_0\in\mathbb{N}$ and a constant $C$ such that

for all $N > N_0.$

Proof. We prove this statement for the 1-D case, and the proof of the $\nu$-dimensional case is similar. Let $\|y\| < 1.$

Setting $t = t_n+vh_n$ and $\lvert y(z)\rvert = \max \lvert y(t_{n, j})\rvert,$ we obtain

where

is independent of $y$ and $t.$

Remark 7. We note that for $m = 1,$ we have $\|P_N\|\leq 1.$ Thus, $\Lambda$ in Eqs (6.2) and (5.10) will be equal to $M$ defined in (4.6). This means that the conditions of (6.1), (5.1) and (4.2) coincide.

7.   Numerical experiments

To start numerical implementation, we compute the weight coefficients for the given parameters $m$ and $c$. Let $m = 1$ and $c_1 = \theta,$ $L_1(t) = 1.$ The corresponding methods are known as the $\theta$ method. The $\theta$ method approximates the solution with a constant function in each interval. In uniform mesh, it is similar to Haar scale functions. The weight of the proposed approximate quadrature can be computed by

7.1.   Validation example

Consider a system of GFDEs described by Eq (1.2) with the vector order $\bar{\alpha} = [0.5, 0.5]^T,$ the weight function $\mathbf{g} = [g, g]^T,$ the source function

and the terminal value

Before giving the exact solution, we notice that $\mathbf{f}$ satisfies the condition (H1) and (H2). According to Theorem 4.2, for an appropriately chosen value of $b,$ the solution falls into $(C[a, b])^2.$ The exact solution of this system is

This is in agreement with the assertion of Theorem 4.2.

We notice that we could not choose $a = 0$ for case II, since the function $\ln$ is not defined at zero. For a better comparison, we choose the initial terminal $a = 1,$ and we find $b$ such that

Equation (7.1) guarantees the condition (4.6) of Theorem 4.1. Thus, the corresponding system has a unique solution with an arbitrary value of $\mathbf{y}_b.$ Table 1 shows lower bounds of $b$ for various choices of weight functions $g$ obtained by Eq (7.1).

Remark 8. Table 1 highlights case II (Hadamard fractional operator) as well as case IV (Katugampola type) as candidates for having bigger bounds of well-posedness. The parameters of cases III and IV change $b$ with different rates. The most rapid change occurs in the Katugampola type. We recall that the terminal value problem may not be well-posed ^[4], at all.

Remark 8 guides us toward available options for modeling inverse problems in applied science. Actually, we need well-posedness for larger intervals in general.

In Figure 1, we illustrate the exact and the approximate solutions for various weight functions $g,$ with the upper bound $b = 2,$ the collocation parameter $c = 0.1,$ the graded mesh exponent $r = 2$ and the number of steps $N = 20$. The maximum errors and estimated orders of convergence for both components of the solution are compared in Tables 2 and 3, respectively. The maximum error is computed by

It is known that, for $\mathbf{y}\in (C[a, b])^2,$ $R(\mathbf{y}) = O(N^{-m})+O(N^{r\beta})$ (see for example ^[9]). In this example, $m = 1$ and $r = 2;$ thus, according to Theorem 6.1,

Tables 2 and 3 show that the order of converge tends to $1.$ Thus, the reported results are in complete agreement with the theoretical results.

7.2.   Past population history of diabetes cases

Diabetes is a silent epidemic of high blood sugar levels, mainly due to decreasing activities and an inappropriate diet. An epidemiology model using FADEs for simulating the population of diabetes cases is proposed in ^[31]. This model applies linear two-dimensional FDEs

to describe the population of diabetes cases. Here,

is a $2\times 2$ matrix, and $B$ is a 2-dimensional column vector. By some physical interpretation of the model ^[31], we impose $0\leq a_{11} < 5,$ $0\leq a_{12}\leq 1$ and $0\leq a_{21} < 2.$ The interpretation of $a_{22}$ is the rate of death in ^[31]. However, the population of the studied zone is not a constant parameter, and indeed it is increasing. Thus, we let $a_{22}$ be positive, and we impose $-1\leq a_{22} < 8.$

Let us use the system (7.2) as a black box for modeling the populations of diabetes cases with model parameters $A$ and $B.$ We use available data, which was reported in ^[32], for the population of diabetes cases from $1990$–$2017$ to find the model parameters. Then, we use the obtained model for predicting the past dynamics by a terminal value problem.

We set $\bar{\alpha} = [0.9, 0.9]^T$ to involve history in our model. Our interest is to study the effect of $g$ on the results of this model. Let the output of the black box be $o(t) = y_2(t),$ indicating the population of diabetes cases for an input value of $t.$

We use $o$ for simulating the dynamics of the prevalence number of diabetes. Let $D(t)$ show available data for the prevalence number of diabetes. Then, $D(t)/D_M$ is normalized data, where $D_M = \max |D(t)|.$

Setting $A = [0, 0;0, 0],$ we obtain

This means, by tending $\alpha\rightarrow 0,$ the matrix $B$ can be regarded as an approximation for the slope of the curve $\mathbf{y}(t).$ Thus, a reasonable physical interpretation suggests choosing $B$ as the prevalence rate of diabetes with some gain related to $g$.

According to ^[32], a prevalence rate is a number between $4500$ and $5500$ per $100k$ population. Thus, we use $B = \eta[0, 0.045]^T$ to reduce the parameters in the related optimization problem, where the parameter $\eta\in \mathbb{R}$ is a model parameter.

Now, the inverse problem of modeling dictates to determining $A$ and $\eta$ such that minimizes the following objective function:

The output $o(t)$ will be calculated from the block box as an initial value problem on $[1990,2017]$. The modeled IVP can be solved by Eq (2.19) or a numerical method described by (5.9). We applied the latter one.

Let's apply the inverse problem to the obtained model. We emphasize that classical fractional-order TVPs (i.e., TVPs with Caputo fractional derivative) for high dimensional systems are not well-posed on the large intervals. A cure for such a modeling problem is the ultimate aim of this paper. To this end, we use appropriate generalized fractional derivatives with a larger well-posedness bound.

Let $\mathbf{g}(x) = [x^{\rho t}, x^{\rho t}]^T.$ The case $\rho = 1$ corresponds to the Caputo fractional derivative. In Figure 2(a) we depict the result of the fitted model for $\rho = 1.$ We apply (5.9) with $N = 200,$ $c = 1$ and graded mesh of exponent $r = 2.$ The fitting error (7.3) in this case is $E = 0.24998.$ Since $\Lambda = 113.2473,$ the conditions (6.2) and (4.6) do not hold. Thus, Theorems 4.2, 5.1 and 6.1 do not guarantee the corresponding numerical and exact solutions exist. Also, there is no guarantee for convergence of the approximate solution. Figure 2(b) shows the approximate solution by the prescribed method in Eqs (5.1) and (5.7) ($N = 200,$ $c = 1$ and $r = 2$) for the corresponding TVP. It is not a surprise that the numerical approximation for the TVP does not converge to the corresponding solution of the IVP. Thus, this model cannot be used for predicting past dynamics.

Mathematically, the constraints for elements of $A$ impose that the Lipschitz constant of $\mathbf{f}$ be less than $8.$ As a result, we can improve the domain of well-posedness by decreasing $\rho$ in $\mathbf{g}(x) = [x^{\rho t}, x^{\rho t}]^T$. On the other hand, this restriction may increase the optimization error, i.e., $E.$ Using the same procedure as we reported for the case $\rho = 1,$ we report the result for $\rho = 0.2, 0.4, 0.8$ in Figures 3–5, respectively (the case $\rho = 0.6$ is similar to case $\rho = 0.8$ and it does not add more information).

For the case $\rho = 0.8,$ the fitting error $E = 0.1530$ is decreased (Figure 3(a)–(c)). The $\lambda$ corresponding to interval $I_1: = [1990,2017]$ is decreased to $13.1412.$ It is not in the zone of the predicted bound for well-posedness and convergence of the given TVP. However, it converges, and the solutions of the TVP and IVP coincide, as shown in Figure 3(b). The solution for prediction of past time on $I_2 = [1940,2017]$ with the same TVP is illustrated in Figure 3(c). Obviously, increasing the interval length increases $\lambda.$ For interval $I_2$ we get $\lambda = 33.8248.$ There is no guarantee of well-posedness, and the results cannot be reliable (a similar pattern repeats for the case $\lambda = 0.6$).

For the case $\rho = 0.4,$ the fitting error $E = 0.37360$ is increased (Figure 4(a)–(c)). However, the values $\lambda$ corresponding to $I_1$ and $I_2$ are decreased, and both are near the well-posedness and convergence zone. We can rely on these results.

The last case is $\rho = 0.2:$ see Figure 5(a)–(c). The optimization result will cross some boundaries of our restricted conditions on $A$. The values of $\lambda$ for both intervals are less than $1.$ Well-posedness of the corresponding TVP and convergence of the proposed numerical method are ensured. However, the error decreases to $E = 0.9076.$ It behaves not better than a model with linear regression.

Conclusively, there is a contrast between the fitting error and the well-posedness condition. A compromise can be made by considering both conditions for a given model. In this case, we propose the model with $\rho = 0.4,$ for estimating past dynamics.

8.   Discussion and conclusions

An extensive application of FDEs is involved in many branches of science, such as biology, epidemiology, medicine and chemistry. Always, these systems are higher dimensional and generally involve some terms of nonlinear Lotka Volterra dynamics with initial conditions. Therefore, as IVPs they are well-posed. However, the inverse problem is not well-posed for higher dimensional equations. It is useless and dangerous to use fractional derivatives for modeling such processes, especially when we have boundary conditions, without considering their well-posedness. Fortunately, this paper's results for modeling with general fractional derivatives open a new window to overcome this problem. However, we encounter another problem: increasing fitting error. For this aim, a compromising solution between well-posedness and less fitting error is considered.

Conflict of interest

The authors declare no conflict of interest.