The computational domain consists of the aortic root, aorta, left coronary artery (with branches), and right coronary artery (with branches). The diameters, lengths, and topology of vessels can be extracted from CCTA scans. A simplified version of arterial network is presented in Figure 1. We simulated stenosis as a separate segment with decreased diameter.

**Figure 1.** A simplified network of major coronary arteries. Segment 6 represents 66% stenosis. The model parameters for each segment are presented in Table A1 in Appendix A. We impose cardiac output function (Figure 2) on the inlet of segment 1. On the terminal ends of segments 4, 7, and 8, we impose hydraulic resistance and outflow pressure (6). Boundary condition on the terminal end of the aorta (segment 2) involves a 2-element Windkessel model (7) .

One-dimensional vessels are connected to each other in junction points to create an arterial structure. The conditions of mass conservation and total pressure continuity are imposed at the junction points:

$$\sum\_{i} Q\_{i} = 0,\tag{4}$$

$$\frac{\rho u\_i^2}{2} + P\_i = \frac{\rho u\_j^2}{2} + P\_{j\prime} \quad i \neq j. \tag{5}$$

Equation (4) represents an algebraic sum of influxes and effluxes, where *i* is the index of a vessel connected to a junction. *ui*, *uj* and *Pi*, and *Pj* in (5) are the velocities and blood pressures of vessels with indices *i* and *j* near the junction point.

**Figure 2.** Cardiac output function.

On the inlet of the aorta (segment 1 on Figure 1), we set cardiac the output as a periodic time function *Q*(*t*) (Figure 2). The shape of the cardiac output was proposed in [23]. It can be adjusted according to the patient's heart rate (HR), stroke volume (SV), or other data (peak-to-mean flow ratio, QT-interval).

On the outlet of coronary arteries (segments 4,7,8 on Figure 1), we impose a pressure drop condition:

$$P\_k - P\_{out} = R\_k Q\_{k'} \tag{6}$$

where *k* is the index of a segment, *Pk* is the blood pressure at the boundary point, *Qk* is the blood flux at the boundary point, and *Rk* is the hydraulic resistance, *Pout* is the outflow pressure. Outflow pressure can be described as a value of blood pressure at which the microcirculation between arteries and veins stops. It ranges between 20 and 60 mm Hg [5] and can be adjusted according to the patient's data. Resistances *Rk* are distributed according to empiric Murray's law through an algorithm described in [24]. Resistances increase during the systolic phase to simulate contractions of myocardium tissue that hinder coronary blood flow [19].

The boundary condition on the outlet of the aorta (segment 2 on Figure 1) differs from boundary condition on the terminal coronary arteries (6) since the former one represents the whole systemic circle as well as microcirculation. In order to describe the behavior of the microcirculation vessels, the two-element Windkessel model [25] is extensively used:

$$Q\_a(t) = \frac{P\_a(t) - P\_{out}}{R\_a} + C\frac{dP\_a}{dt} \,\prime \tag{7}$$

where *Qa* and *Pa* are the blood flow and pressure in the aorta, and *Ra* is the hydraulic resistance of the systemic circle and microcirculation. In (7), compliance *C* is introduced. which represents the ability of blood vessels to distend and store blood volume. Larger values of *C* correspond to greater vessel elasticity. Compliance *C* can be adjusted according to patient's systolic and diastolic blood pressures, and resistance *Ra* can be derived from the systemic vascular resistance—the ratio between mean blood pressure and cardiac output [24]. The elastic index *c* from (3), on the other hand, increases with the rigidity of the vessels and, thus, has a different physiological interpretation.

We solve the hyperbolic system (1) and (2) inside each vessel numerically with the help of an explicit grid-characteristic method [26], which is monotone and first-order accurate. Compatibility conditions imposed on junctions with Equations (4) and (5) and boundary points with conditions (6) and (7) form the system of nonlinear equations which is solved with the Newton method. Compatibility conditions are discretized implicitly with a first-order approximation. Discretizations and convergence studies are presented in [27].

The described model can be used to calculate FFR at any point of the coronary arteries. We calculate FFR as a ratio between mean pressure in the coronary artery distal to stenosis (*Ph dist*) and mean aortic pressure (*Ph aortic*) during vasodilation of the coronary vessels (hyperemia) [2]:

$$FFR = \frac{\overline{\mathcal{P}}\_{dist}^{\text{fr}}}{\overline{\mathcal{P}}\_{aortic}^{\text{fr}}}.\tag{8}$$

Hyperemia is simulated by decreasing terminal resistances *Rk* by 70% [28]. FFR values below 0.8 are considered to be significant. This means that the coronary vessel has an abnormal narrowing (stenosis) that should be surgically treated.

By adjusting *C* and *Pout* in (6) and (7), we can reproduce patient's systolic and diastolic blood pressures. However, the range of possible blood profiles is quite limited. To improve the mathematical model, additional elements can be introduced into the arterial structure and the Windkessel model [8]. Instead, in this paper, we propose to use the fractional time derivative in Equation (7).

#### *2.2. Fractional-Order Boundary Conditions*

We impose boundary condition on the terminal end of the aorta using the fractionalorder Windkessel model, which can be written as follows:

$$Q\_{\mathfrak{a}}(t) = \frac{P\_{\mathfrak{a}}(t) - P\_{out}}{R\_{\mathfrak{a}}} + \mathcal{C}\_{\mathfrak{a}} D\_t^{\mathfrak{a}} P\_{\mathfrak{a}}(t), \tag{9}$$

where *Dα t* is a fractional differentiation operator; *α* is a fractional differentiation order, which is assumed to be between 0 and 2 in this work; and *Cα* is a pseudo-compliance (pseudo-capacitance). Fractional differentiation order *α* determines the relative degree of interaction between the capacitance of the microvasculature vessels, elastic compliance of the vessels, and the dissipation forces inside them. This, in turn, defines the physiological meaning of *Cα*.

Windkessel models with fractional derivatives were extensively studied before [13,29]. However, combining a one-dimensional hemodynamic model with the fractional-order Windkessel boundary condition is a new approach that allows us to represent a greater variety of storage and dissipation effects that can be represented with a single additional parameter *α*.

A large number of different definitions have been proposed for the fractional differentiation operator *Dα t* [30]. We use the Caputo fractional derivative in this work:

$$D\_t^a P(t) = \frac{1}{\Gamma(\lceil \pi \rceil - a)} \int\_0^t \frac{P^{(\lceil \pi \rceil)}(t')}{(t - t')^{1 + a - \lceil \pi \rceil}} dt',\tag{10}$$

where Γ is a gamma function, *α* is a ceiling of *α* (smallest integer greater than *α*), and *P*(*α*)(*t* ) is a derivative with an integer order *α*. There are many other definitions of the fractional derivative: the Atangana–Baleanu fractional integral [31], Riemann–Liouville fractional derivative [32], Riesz derivative [33], etc. The choice of the Caputo fractional derivative is due to its simplicity in representation (relative to other fractional derivatives) and availability of well-studied numerical methods with approximation estimates for various problems. Another useful representation derived in [34] can be obtained using integration by parts in (10):

$$D\_t^a P(t) = \frac{1}{\Gamma(-a)} \int\_0^t \frac{P(t')}{(t - t')^{1 + a}} dt'.\tag{11}$$

This representation is valid for 0 < *α* < 2. We use it to approximate the integral in (11) with a trapezoidal rule. For the interval [0, *t*] with a grid {*t n* = *nτ* : *n* = 0, 1, 2, .., *N*}, assuming a constant time-step *τ* and *P*(0) = 0, *P* (0) = 0, numerical approximation of *Dα τPN* can be expressed as

$$\_rD\_\tau^\kappa P\_N = \frac{\tau^{-\kappa}}{\Gamma(2-\alpha)} \sum\_{n=0}^N a\_{n,N} P\_{N-n\prime} \tag{12}$$

where coefficients *an*,*N* are

$$a\_{n,N} = \begin{cases} 1, & n = 0, \\ (n+1)^{1-a} - 2n^{1-a} + (n-1)^{1-a}, & 0 < n < N, \\ (1-a)N^{-a} - N^{1-a} + (N-1)^{1-a}, & n = N. \end{cases}$$

The error and stability analysis of this numerical approximation was described in detail in [30,35]. It was demonstrated that the order of approximation error is *O*(*τ*2−*α*).

We use an explicit approximation scheme for (9) with the discretization of fractional derivative described above. This allows us to determine the blood pressure *PN* at the terminal end of the aorta from the values of pressure and flux at the previous time steps. Then, we calculate the outflow *QN* at the terminal end of the aorta using compatibility conditions (1) and (2) and wall-state Equation (3).

#### *2.3. Model Personalization*

One of the most important problems in patient-specific blood flow modeling is to identify the parameters of the model. Diameters, lengths, and overall arterial structure can be extracted from CCTA images with the help of segmentation and skeletonization algorithms [36]. Parameter *c* in (3) represents the pulse wave velocity and can be estimated from the patient's age, blood pressure, heart rate, and stroke volume with the help of machine learning methods [20]. Terminal resistances *Rk*, *Ra* in (6), (7) and (9) are calculated from systemic vascular resistance (the ratio between mean pressure and cardiac output) and the diameters of the terminal arteries [24].

Outflow pressure *Pout* and compliance *C* in (7) are estimated to reproduce measured systolic and diastolic blood pressure. A number of algorithms for *Pout* and *C* estimation are presented in [25]. The following procedure to estimate these parameters is used in this work. We calculate initial value of *C* as a ratio between stroke volume and pulse pressure (*PP* = *Psys* − *Pdia*) and the initial value of *Pout* as 50% of diastolic pressure. After this, *C* and *Pout* are iteratively adjusted until the measured diastolic and systolic blood pressures match:

$$P\_{out}^{i+1} = P\_{out}^i \frac{P\_{sys}^{trur} + P\_{dia}^{trur}}{P\_{sys}^i + P\_{dia}^i}, \\ \mathcal{C}^{i+1} = \mathcal{C}^i \frac{P\_{sys}^{trur} - P\_{dia}^{trur}}{P\_{sys}^i - P\_{dia}^i}. \tag{13}$$

Adjustment of *Cα* in (9) is performed using the same procedure, but initial estimation is usually less precise since the dimension and interpretation of *Cα* changes with *α*.

Unfortunately, relying solely on systolic and diastolic pressures may produce incorrect diagnostic outcomes. For example, the hemodynamic significance of stenosis may vary for the same values of systolic and diastolic pressures. In order to obtain more accurate estimates, additional available information about the pressure profile, such as mean pressure, must be taken into account.

We propose to estimate an additional parameter, fractional derivative order *α*, based on the value of the mean pressure *Pmean*. To do this, we first iteratively adjust *α* to match the mean pressure and then adjust *Cα* and *Pout* for each *α* to match the systolic and diastolic blood pressures. After this, we calculate *Pmean* and compare it with the measured mean pressure. If the calculated mean pressure is higher than the measured one, the order *α* should be decreased, and vice-versa. As we will see from the results, the relationship between *α* and *Pmean* is very close to linear. Therefore, in most situations, it is sufficient to perform two preliminary calculations for *α* = 1 and *α* = 1.5 or *α* = 0.5 to estimate *α*, which provides the appropriate value for the mean pressure.

The procedure of parameter identification described above is summarized on Figure 3.

**Figure 3.** Parameter identification procedure.

#### *2.4. Patient Data*

We tested the parameter identification procedure (Figure 3) on a publicly available dataset [5]. This dataset contains data from ten patients, including the geometry of arterial networks and the location stenoses in various coronary vessels. We kept the numeration of patients from [5] but we excluded Patient 9 from our study since the measurement of mean pressure is unavailable. As a result, our study included nine patients (Table 1 with 13 stenoses).

**Table 1.** Characteristics of the patient dataset (mean ± standard deviation). Details are presented in [5]. *θ* = *Pmean*−*Pdia Psys*−*Pdia* is a measure of blood profile thickness

Figure 4 presents examples of two patient-specific structures for Patient 1 and Patient 4. Table 2 presents patient metadata as well as measured FFR values. These two patients were chosen based on the value *θ*, which is the measure of the blood profile thickness:

$$\theta = \frac{P\_{\text{mean}} - P\_{\text{dia}}}{P\_{\text{sys}} - P\_{\text{dia}}}.\tag{14}$$

The average value of *θ* across all nine patients is *θave* = 0.45. Patient 1 has the lowest value *θ*1 = 0.36, and Patient 4 has the closest to the average value *θ*4 = 0.44. As a result, Patient 1 has the thinnest blood pressure profile and Patient 4 has the most typical profile.

**Figure 4.** Patient-specific network of coronary arteries: (**a**) patient 1 with 70% stenosis (segment 5) and (**b**) patient 4 with prolonged 50–60% stenosis (segment 5). Parameters of each segment are presented in Tables A2 and A3.

**Table 2.** Characteristics of Patient 1 and Patient 4. Stroke volume (SV) was estimated from patient age and BMI values presented in [5].

#### **3. Results**

#### *3.1. Blood Pressure and FFR Sensitivity to Order α*

The proposed model was applied to calculate the blood pressure profiles for various fractional differentiation orders *α* in a simplified network of coronary arteries (Figure 1). We also calculated the FFR for 66% of the stenoses for various *α*.

First, we calculated aortic blood pressure for *α* = 1.0 and adjusted the model parameters to acquire the physiological systolic (125 mm Hg) and diastolic (75 mm Hg) blood pressures. Then, we performed calculations for other values of *α* with the same set of model parameters. The value of compliance *C* remains constant for various *α*. This is technically incorrect since the physiological meaning and dimensional formula for *C* depend on *α*, but it helps us to explore changes in pressure profiles with the change of *α* (Figure 5a).

Second, we adjusted *C* and *Pout* for each order *α* to obtain the same values of systolic and diastolic blood pressures (Figure 5b). As *α* decreases, the pressure peak shifts to the left, and the pressure profile becomes thinner. This results in a drop in mean pressure. Conversely, as *α* increases, the profile peak shifts to the right, the profile becomes thicker, and the mean pressure grows.

**Figure 5.** Aortic blood pressure for various orders *α*. (**a**) All parameters, except for *α*, are fixed. (**b**) Model parameters were adjusted to yield the same values of the systolic and diastolic blood pressure.

Systolic and diastolic blood pressures are one of the most commonly available patientspecific parameters. All blood profiles on Figure 5b have the same systolic and diastolic blood pressures, but the mean pressures are different for each *α*. If the mean pressure is available for a given patient, we can use it to select an appropriate *α* and calculate FFR. Figure 6 demonstrates how the calculated mean pressure and FFR depend on *α*, assuming that the systolic and diastolic blood pressures are the same (125/75 mm Hg). The relationship between mean pressure and *α* is very close to linear within the considered interval (from 0.25 to 1.5). Therefore, this simplifies the process of *α* identification from patient's mean pressure: if we calculate the mean pressures for any two values of *α*, we can derive an appropriate *α* for a given mean pressure value using linear interpolation.

At the same time, FFR drops with increasing *α* (Figure 6b). This relation resembles exponential decay. For 0 < *α* < 1, FFR drops rapidly from 0.95 (*α* = 0.25) to 0.78 (*α* = 1.0). The threshold between the significant and insignificant lesions is 0.8, so the choice of *α* affects the diagnostic outcome. For 1 < *α* < 2, FFR is almost constant.

**Figure 6.** Mean pressure and FFR for a simplified network of coronary arteries. (**a**) Mean pressure and *α*. (**b**) FFR and *α*. The horizontal red line corresponds to FFR = 0.8—threshold value between significant (FFR < 0.8) and insignificant (FFR > 0.8) stenoses .

#### *3.2. Patient-Specific Calculations*

In this section, we describe applying a parameter estimation algorithm (Figure 3) to estimate FFR for nine patients (Table 1). We start with two examples: Patient 1 and Patient 4.

Patient 1 had stenosis in the left anterior descending artery (LAD) with a corresponding measured FFR value of 0.89 and a mean aortic pressure of 111 mm Hg. Blood flow calculations for *α* = 1 (utilizing boundary condition (7)) yielded *FFRα*=1 = 0.83 and *Pα*=1 *mean* = 119 mm Hg. Then, we applied the proposed model with fractional-order Windkessel boundary conditions (9) with the parameter identification procedure to fit the given systolic, diastolic, and mean pressures. As a result, the following estimates were obtained: *Pα*=0.56 *mean* = 111 mm Hg was achieved for *α* = 0.56, and the resulting FFR was *FFRα*=0.56 = 0.87. The error in FFR estimation was significantly reduced after applying optimization of *α* for mean pressure.

Patient 4 had stenosis in the LAD with a corresponding measured FFR value of 0.82 and a mean aortic pressure of 94 mm Hg. Blood flow calculations for *α* = 1 (utilizing boundary condition (7)) yielded *FFRα*=1 = 0.82 and *Pα*=1 *mean* = 94 mm Hg. No further optimization was required in this case since the calculated mean pressure matched the measured mean pressure with good accuracy. The calculated FFR value also matched the measured one.

FFR estimations for other patients are presented in Table 3. The original approach to estimate FFR involves a boundary condition (7) without fractional derivative. We ignored the measured mean pressure and adjusted our model to achieve the measured systolic and diastolic blood pressures. The fractional-order approach involves adjusting the fractional derivative order to obtain measured mean pressure. The RMSE for the original approach was 0.05, and the RMSE for the fractional order approach was 0.04. The RMSE was mainly defined by large errors in the FFR estimations of patients 6, 8, and 10. We assumed that stenosis degree and length were not identified properly for these patients. The FFR estimation was improved for patients with a "thin" blood profile (*θ* < 0.4), including Patient 1, Patient 5, and Patient 7. Patients 2, 3, 4, and 6 had am optimal fractional order *αopt* close to 1.0 (or equal to 1.0), so FFR estimations for both approaches were similar. The FFR estimations for patient 6 were less precise for the fractional order approach, but the difference was very small. Patients 8 and 10 had an optimal fractional order *αopt* > 1.0, and the FFR estimation was similar for both approaches. This is due to the fact that FFR is almost constant for *α* > 1.0 (Figure 6).

**Table 3.** FFR estimations for the patient dataset. Patient data: *Pmean* is the measured mean pressure, mm Hg; *θ* is a measure of the pressure profile thickness (14); Loc. is a location of stenosis; *FFR* is the invasively measured FFR. The original approach for FFR estimation (order *α* = 1.0): *FFRα*=1 is the calculated FFR value with a boundary condition (7); *Pest mean* is the calculated mean pressure with order *α* = 1. The fractional derivative approach (order *α* = *αopt*) involves adjusting order *α* so that the calculated mean pressure matches the measured one: *FFRα*=*αopt* is the calculated FFR value with the boundary condition (9), and fractional order *α* = *αopt*; *αopt* is the optimal fractional order. Patient 9 was excluded due to the absence of a mean pressure measurement.

Results of the FFR estimation showed that the fractional-order approach provided benefits for a certain group of patients (patients 1, 5, and 7) with a thin blood pressure profile. These patients may have a common cardiovascular pathology. According to [14], low fractional orders can be used to simulate the blood profiles of patients with hypertension. Patients 1, 5, and 7 have a BMI > 30 kg/m2 which is a good predictor of hypertension. However, their blood pressure levels are not the highest in the dataset. We need a larger sample of patients to make further conclusions.

#### **4. Discussion**

We proposed coupling the well-established one-dimensional hemodynamic model of coronary blood flow with a fractional-order boundary condition, as well as a procedure for estimating its parameters. The fractional derivative can be a useful tool for more accurate modeling of the pressure profile. The actual blood pressure profiles depend on many factors, such as age, height, weight, medical history, and artery elasticity. The most commonly available characteristics of blood pressure profiles are systolic and diastolic blood pressures. Unfortunately, relying solely on these two values alone may produce incorrect diagnostic outcomes. For example, the hemodynamic significance of stenosis may vary for the same systolic and diastolic blood pressures. This fact has motivated researchers to look for new tools to model blood pressure profiles.

We used the fractional derivative order to match calculated mean blood pressure with the measured one. Adjusting mean pressure can be performed in other ways: introducing a larger Windkessel model, expanding the arterial network, etc. The fractional-derivative Windkessel model is a good compromise: we introduced a very small amount of new parameters (1–2) and gained the ability to simulate a whole spectrum of dissipative and storage mechanisms with the help of fractional order *α*. Our approach does require additional information on the blood pressure profile. These data can be obtained with simple noninvasive procedures. Unfortunately, in many cases, these data are unavailable, and all the pressure profile information is reduced to systolic and diastolic pressure values. This was the case for Patient 9 who was excluded from our study.

The proposed approach has a number of shortcomings that need to be resolved in the future. First, in some cases, the only data available regarding a patient's pressure profile were systolic and diastolic blood pressures. Identifying the fractional derivative order *α* in this case would require a completely different approach that can be based on other data, such as patient medical history. Second, calculating fractional derivatives calls for significant computational resources. This negates one of the main advantages of one-dimensional blood flow models—computational efficiency. Instead of the basic approach presented in this work, new efficient numerical approximations can be used. Third, the proposed approach should be tested on a larger number of patients with a wide range of FFR values. The FFR is nearly constant for fractional differentiation orders *α* > 1, so for some patients, its adjustment will be useless.

Further research will focus on integrating more efficient approaches to identify model parameters and to approximate a fractional derivative. Other areas of work include collaborating with clinicians to find effective methods for pulse wave assessment and further validation on a larger amount of patients. The proposed approach has great potential to provide an alternative means to simulating arterial stiffness and pulse waves.

**Author Contributions:** Conceptualization, T.G. and R.Y.; methodology, T.G.; software, T.G.; validation, T.G.; formal analysis, R.Y.; investigation, T.G.; data curation, T.G.; writing—original draft preparation, T.G.; writing—review and editing, R.Y.; visualization, T.G. and R.Y.; supervision, T.G. All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by the Russian Science Foundation under project no. 22-71-10087.

**Data Availability Statement:** The datasets used in this work are described in [5] and are freely available at https://doi.org/10.6084/m9.figshare.8047742.v2 (accessed on 2 April 2023).

**Acknowledgments:** We thank Alexander Lapin (Sechenov University) for useful discussion and literature recommendations.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**

The following abbreviations are used in this manuscript:

#### **Appendix A**

**Table A1.** Parameters of the vessels for simplified structure (Figure 1).

**Table A2.** Parameters of the vessels for Patient 1 (Figure 4).

**Table A3.** Parameters of the vessels for Patient 4 (Figure 4).

#### **References**

**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.

## *Article* **An Application of the Distributed-Order Time- and Space-Fractional Diffusion-Wave Equation for Studying Anomalous Transport in Comb Structures**

**Lin Liu 1,\*, Sen Zhang 1, Siyu Chen 1, Fawang Liu 2,\*, Libo Feng 2, Ian Turner 2, Liancun Zheng 1 and Jing Zhu 1**

**Abstract:** A comb structure consists of a one-dimensional backbone with lateral branches. These structures have widespread application in medicine and biology. Such a structure promotes an anomalous diffusion process along the backbone (*x*-direction), along with classical diffusion along the branches (*y*-direction). In this work, we propose a distributed-order time- and space-fractional diffusion-wave equation to model a comb structure in the more general setting. The distributed-order time- and space-fractional diffusion-wave equation is firstly formulated to study the anomalous diffusion in the comb model subject to an irregular convex domain with the motivation that the time-fractional derivative considers the memory characteristic and the space one with the variable diffusion coefficient possesses the nonlocal characteristic. The finite element method is applied to obtain the numerical solution. The stability and convergence of the numerical discretization scheme are derived and analyzed. Two numerical examples of relevance to the comb model are given to verify the correctness of the numerical method. Moreover, the influence of the involved parameters on the three-dimensional and axial projection drawing particle distribution subject to an elliptical domain are analyzed, and the physical meanings are interpreted in detail.

**Keywords:** distributed-order fractional derivative; anomalous diffusion; comb model; constitutive relationship

**PACS:** 26A33; 65M12; 65M60; 35R11; 74Q15

### **1. Introduction**

A comb model is used to study anomalous diffusion in a medium of a specific structure [1]. Systematic research on this class of models is of great theoretical significance and application to comb structures such as dendritic spines [2] that arise in medicine and biology. An example of an experimental setup to probe the dynamics of actin polymerization is given in Figure 1a. The image at the top gives the optical micrograph of the microfluidics structure, and the image at the bottom is of the microfluidic micrographs fluorescently labeled, polymerized actin filaments [3]. Figure 1b presents the electron tomogram of a spiny dendrite [4]. From the two practical problems, it is easily seen that the particle transport is not random but in the form of a comb, and this specific structure is named a comb model [5]. Comb models are a powerful tool for studying many other diffusion phenomena, such as the diffusion of cancer cells [6], the fractal glioma development under RF electric field treatment [7] and the diffusion of percolation clusters [8]. For this study, we assume for all the practical problems that the transport of particles in a comb form can be simplified to the structure exhibited in Figure 2. As the figure shows, the comb model contains a straight backbone on the *x*-axis with the lateral branches attached perpendicularly to the backbone, which plays

**Citation:** Liu, L.; Zhang, S.; Chen, S.; Liu, F.; Feng, L.; Turner, I.; Zheng, L.; Zhu, J. An Application of the Distributed-Order Time- and Space-Fractional Diffusion-Wave Equation for Studying Anomalous Transport in Comb Structures. *Fractal Fract.* **2023**, *7*, 239. https://doi.org/ 10.3390/fractalfract7030239

Academic Editors: Stanislaw Migorski and Riccardo Caponetto

Received: 30 December 2022 Revised: 21 January 2023 Accepted: 10 February 2023 Published: 7 March 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

the role of traps [8]. Through the special structure of the comb model, the diffusion process along the *x*-direction only happens on the *x*-axis and the transport between any different fingers must take place through the backbone. As is well known, one of the most striking characteristics for the classical normal diffusion is the linear growth with time of the second moment of the particle positions. For the special structure of the comb model, one important pattern of diffusion can be deduced, whereby the subdiffusion with exponent 1/2 arises subject to the classical Fick's constitution model with a linear form [9]

$$\vec{J}(\mathbf{x}, y, t) = \left( -k\_x \delta(y) \frac{\partial u(\mathbf{x}, y, t)}{\partial \mathbf{x}}, -k\_y \frac{\partial u(\mathbf{x}, y, t)}{\partial y} \right), \tag{1}$$

where *J*(*x*, *y*, *t*) refers to the diffusion flux vector, *u*(*x*, *y*, *t*) denotes the distribution function at the special positions (*x*, *y*) and time *t*, *kxδ*(*y*) indicates the diffusion coefficient along the *x*-direction while *ky* is the diffusion coefficient along the *y*-direction, *δ*(*y*) refers to the Dirac delta function which reflects the special structure of the comb model.

**Figure 1.** (**a**) Optical micrograph of a segment of the microfluidics comb-like structures (on top). On the bottom: microfluidic micrographs of fluorescently labeled, polymerized actin filaments in a comb-like structure [4]. (**b**) One of the examples of a physical environment suitable for the comb model. Electron tomogram of a spiny dendrite. Image taken from Internet (http://www.cacr.caltech. edu/projects/ldviz/results/levelsets/, accessed on 1 July 2013).

**Figure 2.** The schematic of a comb model.

Due to the geometrical structure and the non-uniformity of the medium transmission, the classical Fick's law in conventional diffusion with the paradox of an infinite transport velocity [10,11] is no longer applicable. In order to study the transmission mechanism of the concentration field for the anomalous diffusion in the comb model, three modifications for the Fick's model are presented. The first is the introduction of the relaxation parameter *ξ* and the second considers the time-fractional derivatives with the motivation that the relaxation parameter makes the transport process attach a finite transport velocity, while the time-fractional derivative takes the memory characteristic of the transport process into account. Furthermore, as discussed in [12], due to the special structure of the comb model, the highly inhomogeneous characteristic happens along the *x*-axis, and this characteristic can be reflected by the space-fractional derivative [13]. Thus, as the third modification, the second space-integer derivative is modified to a space-fractional derivative with variable coefficients, which considers a left and right nonlocal characteristic. Based on the time-fractional Cattaneo model [14] and the two-dimensional space-fractional constitutive model [15], the following time- and space-fractional Cattaneo constitution relationship with variable diffusion coefficients is formulated

$$\begin{split} \vec{f}(\mathbf{x}, \boldsymbol{y}, t) + \xi \frac{\partial \vec{f}(\mathbf{x}, \boldsymbol{y}, t)}{\partial t} \\ = \xi\_0^{RL} D\_t^{1-a} \left( -\delta(\boldsymbol{y}) \left( d^+(\mathbf{x}, \boldsymbol{y})\_\perp D\_x^\gamma \boldsymbol{u}(\mathbf{x}, \boldsymbol{y}, t) - d^-(\mathbf{x}, \boldsymbol{y})\_\mathcal{X} D\_\mathcal{K}^\gamma \boldsymbol{u}(\mathbf{x}, \boldsymbol{y}, t) \right), -e(\mathbf{x}, \boldsymbol{y}) \frac{\partial \boldsymbol{u}(\mathbf{x}, \boldsymbol{y}, t)}{\partial \boldsymbol{y}} \right), \end{split} \tag{2}$$

where *δ*(*y*)*d*+(*x*, *y*) and *δ*(*y*)*d*−(*x*, *y*) represent the left and right variable diffusion coefficient along the *x*-direction, respectively; *e*(*x*, *y*) refers to the variable diffusion coefficient along the *y*-direction; *RL* 0 *D*1−*α t* refers to time-fractional derivative of order 1 − *α* (0 < *α* < 1) with the Riemann–Liouville definition; *LDγ x* and *xDγ R*, respectively, denote the left and right Riemann–Liouville space-fractional derivatives of order *γ* (0 < *γ* < 1). The definitions are, respectively, given by

$$\begin{split} \, \_0^{RL}D\_t^{1-\alpha}u(x,y,t) &= \frac{1}{\Gamma(\alpha)} \frac{\partial}{\partial t} \int\_0^t \frac{1}{(t-\tau)^{1-\alpha}} u(x,y,\tau) d\tau, \\ \, \_LD\_x^{\gamma}u(x,y,t) &= \frac{1}{\Gamma(1-\gamma)} \frac{\partial}{\partial x} \int\_L^x \frac{1}{(x-s)^{\overline{\gamma}}} u(s,y,t) ds, \\ \, \_xD\_R^{\gamma}u(x,y,t) &= \frac{-1}{\Gamma(1-\gamma)} \frac{\partial}{\partial x} \int\_x^R \frac{1}{(s-x)^{\overline{\gamma}}} u(s,y,t) ds, \end{split}$$

where the symbol Γ(·) refers to the Euler gamma function, *L* and *R* refer to the left and right boundaries along the *x*-direction.

By combining the constitutive relation (2) with the following mass conservation equation

$$\frac{\partial u(\mathbf{x}, \mathbf{y}, t)}{\partial t} + \nabla \cdot \vec{J}(\mathbf{x}, \mathbf{y}, t) = 0,\tag{3}$$

we obtain the time- and space-fractional Cattaneo governing equation

$$\begin{aligned} &\xi \frac{\partial^{1+a}u(\mathbf{x},y,t)}{\partial t^{1+a}} + \frac{\partial^a u(\mathbf{x},y,t)}{\partial t^a} \\ &= \frac{\partial}{\partial x} \left[ \delta(y) \left( d^+(\mathbf{x},y)\_L D\_X^\gamma u(\mathbf{x},y,t) - d^-(\mathbf{x},y)\_X D\_X^\gamma u(\mathbf{x},y,t) \right) \right] + \frac{\partial}{\partial y} \left[ \epsilon(\mathbf{x},y) \frac{\partial u(\mathbf{x},y,t)}{\partial y} \right], \end{aligned} \tag{4}$$

where the symbols *∂αu*(*x*,*y*,*t*) *∂tα dα* and *∂*1+*αu*(*x*,*y*,*t*) *∂t*1+*α* denote the Caputo fractional derivative operators of order *α* and 1 + *α*, respectively, in which the definitions are given as

$$\frac{\partial^{\alpha}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial t^{\alpha}} = \frac{1}{\Gamma(1-\alpha)} \int\_{0}^{t} \frac{1}{(t-\tau)^{\alpha}} \frac{\partial u(\boldsymbol{x},\boldsymbol{y},\tau)}{\partial \tau} d\tau,$$

$$\frac{\partial^{1+\alpha}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial t^{1+\alpha}} = \frac{1}{\Gamma(1-\alpha)} \int\_{0}^{t} \frac{1}{(t-\tau)^{\alpha}} \frac{\partial^{2}u(\boldsymbol{x},\boldsymbol{y},\tau)}{\partial \tau^{2}} d\tau.$$

As a generalization of the integer derivative, the fractional operator considers the memory and nonlocal characteristics and has important applications in a variety of fields. For the fractional governing Equation (4), a limitation is that it is suitable for describing the probability density distribution of a very narrow class of diffusion processes because it is characterized by a unique time-fractional exponent [16]. Motivated by this idea, by integrating the fractional-order derivatives with respect to the order of differentiation, a distribution-order operator was proposed by Caputo [17]. The governing equation with the distributed-order operator exhibits memory and nonlocal effects over various timefractional scales and becomes a powerful tool to describe transport phenomena in complex heterogeneous media. However, as far as we are aware, the distributed-order time and space diffusion-wave equation has not been considered to study the anomalous diffusion in the comb model.

Motivated by the above discussions, as an original contribution to the literature, we discuss and analyze the following distributed-order time- and space-fractional diffusionwave equation to study the anomalous diffusion in the comb model

*ξ* - 2 1 *ϕ*1(*β*) *∂βu*(*x*, *y*, *t*) *∂tβ dβ* + - 1 0 *ϕ*0(*α*) *∂αu*(*x*, *y*, *t*) *∂tα dα* (5) = *∂ ∂x δ*(*y*) *d*+(*x*, *y*)*LDγ x u*(*x*, *y*, *t*) − *d*−(*x*, *y*)*xDγ Ru*(*x*, *y*, *t*) + *∂ ∂y e*(*x*, *y*) *∂u*(*x*, *y*, *t*) *∂y* + *f*(*x*, *y*, *t*),

where *ϕ*1(*β*) and *ϕ*0(*α*) denote weight functions and *β* ∈ (1, 2), *ϕ*1(*β*) ≥ 0, *ϕ*1(*β*) ≡ 0, 0 < 2 1 *ϕ*1(*β*)*dβ* < ∞, *α* ∈ (0, 1), *ϕ*0(*α*) ≥ 0, *ϕ*0(*α*) ≡ 0, 0 < 1 0 *ϕ*0(*α*)*dα* < ∞ and *f*(*x*, *y*, *t*) is a source term.

As Ref. [6] indicated, the diffusion in the comb model, which is described by the Fick's model, is an example of a subdiffusive one-dimensional medium where a continuous-time random walk takes place along the backbone while the diffusion along the *y* direction has a traps effect. The classical Fick's model possesses the local characteristic, and the fractional derivative is proposed, considering the memory and nonlocal characteristics. For a further development, the distributed-order time-fractional derivative is proposed by integrating the fractional-order derivatives. In conclusion, Equation (5) is a development to describe the continuous-time random walk, considering various memory and nonlocal characteristics.

The anomalous diffusion in the comb model has been applied to many different fields in medicine and biology. In some applications of the comb model treating with the numerical method, the infinite regions are approximated by the rectangular domains with large sides [18]. However, by extending the computational modeling to irregular domains, we can broaden the potential applicability of the comb model. Based on these discussions, in this paper, the initial and boundary conditions are given by

$$u(\mathbf{x}, y, 0) = \phi\_0(\mathbf{x}, y), \ u\_t(\mathbf{x}, y, 0) = \phi\_1(\mathbf{x}, y), \ (\mathbf{x}, y) \in \overline{\Omega},\tag{6}$$

$$
\mu(x, y, t) = 0, \ (x, y, t) \in \partial\Omega \times (0, T], \tag{7}
$$

respectively, where Ω is an irregular convex domain.

The new governing Equation (5) is a generalization and development of the classical model to study the anomalous diffusion in the comb model. By choosing *ϕ*0(*α*)=*δ*(*α* − *α*0), *ϕ*1(*β*)=*δ*(*β* − *α*0 − 1), *d*+(*x*, *y*)=*d*−(*x*, *y*) = *kx*, *e*(*x*, *y*) = *ky*, *γ* = 1 and *f*(*x*, *y*, *t*) = 0, Equation (5) reduces to the time-fractional Cattaneo governing equation [14]

$$
\xi \frac{\partial^{a\_0+1}\mu(\mathbf{x},y,t)}{\partial t^{a\_0+1}} + \frac{\partial^{a\_0}\mu(\mathbf{x},y,t)}{\partial t^{a\_0}} = k\_\mathbf{x}\delta(y)\frac{\partial^2\mu(\mathbf{x},y,t)}{\partial \mathbf{x}^2} + k\_\mathbf{y}\frac{\partial^2\mu(\mathbf{x},y,t)}{\partial y^2}.
$$

For the choice *ϕ*0(*α*)=*δ*(*α* − *α*0), *ξ* = 0, *d*+(*x*, *y*)= *kx*, *d*−(*x*, *y*) = 0, *e*(*x*, *y*) = *ky* and *f*(*x*, *y*, *t*) = 0, Equation (5) reduces to the time- and space-fractional governing equation [19]

$$\frac{\partial^{\alpha\_0} u(\mathbf{x}, y, t)}{\partial t^{\alpha\_0}} = k\_x \delta(y) \frac{\partial^{\gamma + 1} u(\mathbf{x}, y, t)}{\partial \mathbf{x}^{\gamma + 1}} + k\_y \frac{\partial^2 u(\mathbf{x}, y, t)}{\partial y^2}.$$

Finally, by choosing *ξ* = 0, *ϕ*0(*α*)=*δ*(*α* − 1), *d*+(*x*, *y*)=*d*−(*x*, *y*) = *kx*, *e*(*x*, *y*) = *ky*, *γ* = 1 and *f*(*x*, *y*, *t*) = 0, we obtain the anomalous diffusion based on the classical Fick's model [2]

$$\frac{\partial \mu(\mathbf{x}, y, t)}{\partial t} = k\_x \delta(y) \frac{\partial^2 \mu(\mathbf{x}, y, t)}{\partial x^2} + k\_y \frac{\partial^2 \mu(\mathbf{x}, y, t)}{\partial y^2}.$$

The study of distributed-order time- and space-fractional diffusion-wave equations is of great significance to further improve and predict the anomalous diffusion phenomena in the comb-like structures.

#### **2. The Structure of the Paper**

In this paper, the two-dimensional irregular convex domain is defined as Ω = {(*x*, *y*)|*xL*(*y*) ≤ *x* ≤ *xR*(*y*), *yD*(*x*) ≤ *y* ≤ *yU*(*x*)}, where *xL*(*y*), *xR*(*y*), *yD*(*x*), *yU*(*x*) are the left, right, lower and upper boundaries of Ω, respectively. Denote *x*min = min (*x*,*y*)∈Ω *xL*(*y*),

*x*max = max (*x*,*y*)∈Ω *xR*(*y*), *y*min = min (*x*,*y*)∈Ω *yD*(*x*), *y*max = max (*x*,*y*)∈Ω *yU*(*x*). Then, the inner product is defined as

$$u(\mu, v)\_{\Omega} = \int\_{y\_{\min}}^{y\_{\max}} \int\_{x\_{\mathbb{L}}(y)}^{x\_{\mathbb{R}}(y)} u(\mathbf{x}, y) v(\mathbf{x}, y) d\mathbf{x} dy = \int\_{x\_{\min}}^{x\_{\max}} \int\_{y\_{\mathbb{D}}(x)}^{y\_{\mathbb{D}}(x)} u(\mathbf{x}, y) v(\mathbf{x}, y) d\mathbf{x} dy,$$

and the *L*2-norm is given as ||*u*||*L*2(Ω) = (*u*, *u*)Ω 1/2.

The fractional derivative space in one dimension was firstly established by Roop and Ervin [20] and then developed further by Bu et al. [21,22], Yang et al. [23], Hao et al. [24] and Wang et al. [25,26]. Due to the special form of the governing equation with a space-fractional derivative of order *γ* in *x* and the integer derivative in *y*, some new definitions and lemmas of the fractional derivative spaces are defined.

**Definition 1.** *The definitions for the left (right) fractional derivative space with semi-norm and norm are, respectively, given as*

$$\begin{split} \left| \left. \left| u \right| \right|\_{\boldsymbol{f}\_{L}^{\boldsymbol{1}}(\Omega)} = \left( \left|| \left.

We denote *J γ*,1 *L*,0 (Ω), *J γ*,1 *R*,0(Ω), *Hγ* 0 (Ω), *H*1 0 (Ω) and *J γ*,1 *S*,0 (Ω) as the closure of *C*∞ 0 (Ω) with respect to || · ||*J γ*,1 *L* (Ω) , || · ||*J γ*,1 *R* (Ω) , || · ||*Hγ*(Ω), || · ||*H*1(Ω) and || · ||*J γ*,1 *S* (Ω) , respectively, where *C*∞ 0 (Ω) is the space of smooth functions with compact support in Ω. Based on the above definitions, some useful and important lemmas are introduced.

**Lemma 1.** *For u*(*x*, *y*), *v*(*x*, *y*) ∈ *J γ*,1 *L*,0 (Ω) B *J γ*,1 *R*,0(Ω) (0 < *γ* < 1)*, we have*

$$\left( \,\_{L}D\_{\mathbf{x}}^{\gamma}u(\mathbf{x},\mathbf{y}), \frac{\partial v(\mathbf{x},\mathbf{y})}{\partial \mathbf{x}} \right) = -\left( \,\_{L}D\_{\mathbf{x}}^{(\gamma+1)/2}u(\mathbf{x},\mathbf{y}) \,, \,\_{\mathbf{x}}D\_{\mathbf{R}}^{(\gamma+1)/2}v(\mathbf{x},\mathbf{y}) \right),\tag{8}$$

$$\left( \,\_x D\_R^{\gamma} u(\mathbf{x}, y), \frac{\partial v(\mathbf{x}, y)}{\partial \mathbf{x}} \right) = - \left( \,\_x D\_R^{(\gamma + 1)/2} u(\mathbf{x}, y) \,\_t D\_x^{(\gamma + 1)/2} v(\mathbf{x}, y) \right). \tag{9}$$

**Proof.** See [27].

**Lemma 2.** *For u* ∈ *Hγ* 0 (Ω) *and* 0 < *γ* < (*γ* + 1)/2*, then*

$$||u||\_{L^{2}(\Omega)} \leq \mathcal{C}\_{2}||\_{L}D\_{x}^{\gamma}u||\_{L^{2}(\Omega)} \leq \mathcal{C}\_{1}||\_{L}D\_{x}^{(\gamma+1)/2}u||\_{L^{2}(\Omega)^{\prime}}\tag{10}$$

$$||u||\_{L^{2}(\Omega)} \leq \mathbb{C}\_{4} ||\_{\mathfrak{x}} D\_{\mathfrak{R}}^{\gamma} u||\_{L^{2}(\Omega)} \leq \mathbb{C}\_{3} ||\_{\mathfrak{x}} D\_{\mathfrak{R}}^{(\gamma+1)/2} u||\_{L^{2}(\Omega)}.\tag{11}$$

*For u* ∈ *H*1 0 (Ω)*, we have*

$$||u||\_{L^{2}(\Omega)} \leq \mathcal{C}\_{5} \left||\left.\frac{\partial u}{\partial y}\right||\_{L^{2}(\Omega)}\right|\,,\tag{12}$$

*where C*1*, C*2*, C*3*, C*4 *and C*5 *are positive constants independent of u.*

**Proof.** The proof of this Lemma follows that given in [28].

**Lemma 3.** *For u*(*x*, *y*) ∈ *J γ*,1 *L*,0 (Ω) ∩ *J γ*,1 *R*,0(Ω)*, we have*

$$\left( {}\_{L}D\_{\mathbf{x}}^{(\gamma+1)/2} \boldsymbol{u}\_{\times \mathbf{x}} D\_{\mathbf{R}}^{(\gamma+1)/2} \boldsymbol{u} \right) = \cos \left( \pi (\gamma+1)/2 \right) \left| \left| {}\_{L}D\_{\mathbf{x}}^{(\gamma+1)/2} \boldsymbol{u} \right| \right|\_{L^{2}(\mathbb{R}^{2})}^{2} \tag{13}$$

$$
\left(\frac{\partial \mu}{\partial y}, \frac{\partial \mu}{\partial y}\right) = \left\| \left| \frac{\partial \mu}{\partial y} \right| \right\|\_{L^2(\mathbb{R}^2)}^2. \tag{14}
$$

**Proof.** Similar to the derivation process in [27], we can obtain the results immediately.

**Lemma 4.** *If u* > 0*, γ* = *n* − 1/2*, n* ∈ N*, then J γ*,1 *L*,0 (Ω)*, J γ*,1 *R*,0(Ω)*, J γ*,1 *S*,0 (Ω)*, Hγ* 0 (Ω) *and H*1 0 (Ω) *are equivalent with equivalent norms and semi-norms.*

**Proof.** The proof of this lemma can be found in [28].

#### **3. Derivation of the Finite Element Scheme for the Comb Model**

*3.1. Finite Element Fully Variational Formulation*

In the following section, for the sake of simplicity, denote *d*1(*x*, *y*) = *δ*(*y*)*d*+(*x*, *y*), *d*2(*x*, *y*) = *δ*(*y*)*d*−(*x*, *y*). Due to the irregular shape of the solution domain, the traditional rectangular mesh cannot be used. The finite element method is applied to obtain the solution of the governing Equation (5), subject to the initial conditions (6) and irregular boundary conditions (7).

Firstly, in the governing Equation (5), the distributed-order time-fractional derivatives are discretized using the mid-point quadrature rule [29]

$$\int\_{1}^{2} \varphi\_{1}(\beta) \frac{\partial^{\beta}}{\partial t^{\beta}} \left( \frac{\partial^{2} u}{\partial x^{2}} \right) d\beta = \sum\_{i=0}^{M\_{1}} \omega\_{i}^{(1)} \frac{\partial^{\beta\_{i}}}{\partial t^{\beta\_{i}}} \left( \frac{\partial^{2} u}{\partial x^{2}} \right) + O(h\_{\beta}^{2}),$$

$$\int\_{0}^{1} \varphi\_{0}(\alpha) \frac{\partial^{\alpha+1}}{\partial t^{\alpha+1}} d\alpha = \sum\_{i=0}^{M\_{2}} \omega\_{i}^{(2)} \frac{\partial^{\alpha\_{i}+1} u}{\partial t^{\alpha\_{i}+1}} + O(h\_{\alpha}^{2}),$$

where *ω*(1) *i* = *hβϕ*1(*αi*), *ω*(2) *i* = *hαϕ*0(*αi*), *hα* = 1 *M*2+1 and *hβ* = 1 *M*1+1 are fractional parameter steps and *αi* = *ihα*+(*i*+1)*hα* 2 for *i* = 0, 1, . . . , *M*2, *βi* = *ihβ*+(*i*+1)*hβ* 2 for *i* = 0, 1, . . . , *M*1.

Let *τ* = *T*/*N* be the time step and *tk* = *kτ* (*k* = 0, 1, ... , *N*) where *N* is a positive integer. Denote *uk*−1/2 = *uk*+*uk*−1 2 for *k* = 1, ... , *N*. For *u*(*x*, *y*, *t*) ∈ *C*(Ω × [0, *T*]), denote *uk* = *uk*(·) = *u*(·, *tk*). We introduce ∇*tuk*−1/2 = *uk*−*uk*−1 *τ* . At *t* = *tk*−1/2, the L2-scheme [30] to approximate the fractional derivative of order *βi* (1 < *βi* < 2) with the Caputo definition is given as

$$\frac{\partial^{[\boldsymbol{\beta}\_{l}]}u^{k-1/2}}{\partial t^{\boldsymbol{\beta}\_{l}}} = \frac{\tau^{1-\boldsymbol{\beta}\_{l}}}{\Gamma(3-\boldsymbol{\beta}\_{l})} \left[ a\_{0}^{\{\boldsymbol{\beta}\_{l}\}} \nabla\_{t}u^{k-1/2} - \sum\_{j=1}^{k-1} (a\_{k-1-j}^{\{\boldsymbol{\beta}\_{l}\}} - a\_{k-j}^{\{\boldsymbol{\beta}\_{l}\}}) \nabla\_{t}u^{j-1/2} - a\_{k-1}^{\{\boldsymbol{\beta}\_{l}\}} u\_{t}^{0} \right] + R\_{0}^{k\_{l}\delta\_{l}} \tag{15}$$

$$= \nabla\_{t}^{\{\boldsymbol{\beta}\_{l}\}}u^{k-1/2} + R\_{0}^{k\_{l}\delta\_{l}}.\tag{15}$$

where *a* (*βi*) *j* = (*j* + 1)2−*βi* − *j* 2−*βi* , *j* = 0, 1, 2, ..., *k* − 1, |*Rk*,*βi* 0 | ≤ *C*(*τ*3−*βi*).

**Lemma 5.** *For the above a* (*βi*) *j , define vector S* = [*S*1, *S*2, *S*3, ..., *SN*] *T and constant P, it holds that*

$$\frac{2\tau^{1-\beta\_l}}{\Gamma(3-\beta\_l)}\sum\_{k=1}^N \left[a\_0^{(\beta\_l)}S\_k - \sum\_{j=1}^{k-1} (a\_{k-j-1}^{(\beta\_l)} - a\_{k-j}^{(\beta\_l)})S\_j - a\_{k-1}^{(\beta\_l)}P\right]S\_k \geq \frac{T^{1-\beta\_l}}{\Gamma(2-\beta\_l)}\sum\_{k=1}^N S\_k^2 - \frac{T^{2-\beta\_l}}{\tau\Gamma(3-\beta\_l)}P^2.$$

**Proof.** See [30].

The time Caputo fractional derivative of order *αi* (0 < *αi* < 1) is discretized by using the L1-scheme [30], and the scheme at *t* = *tk*−1/2 is given as

$$\frac{\partial^{a\_i} u^{k-1/2}}{\partial t^{a\_i}} = \frac{\tau^{1-a\_i}}{2\Gamma(2-a\_i)} \sum\_{j=1}^k b\_{k-j}^{(a\_i)} \nabla\_t u^{j-1/2} + \frac{\tau^{1-a\_i}}{2\Gamma(2-a\_i)} \sum\_{j=1}^{k-1} b\_{k-j-1}^{(a\_i)} \nabla\_t u^{j-1/2} + R\_1^{k, a\_i}$$

$$= \nabla\_t^{(a\_i)} u^{k-\frac{1}{2}} + R\_1^{k, a\_i} \tag{16}$$

where *b* (*αi*) *j* = (*j* + 1)1−*αi* − *j* 1−*αi* , *j* = 0, 1, 2, ..., *k* − 1, |*Rk*,*αi* 1 | ≤ *C*(*τ*2−*αi*).

**Lemma 6.** *For the above b* (*αi*) *j , choosing any positive integer M and vector* [*v*1, *v*2, *v*3, ..., *vM*] ∈ R*M, we have*

$$\sum\_{k=1}^{M} \sum\_{j=1}^{k} b\_{k-j}^{(\alpha\_i)} (\upsilon\_{j'} \, \upsilon\_k) + \sum\_{k=1}^{M} \sum\_{j=1}^{k-1} b\_{k-j-1}^{(\alpha\_i)} (\upsilon\_{j'} \, \upsilon\_k) \ge 0.1$$

**Proof.** See [31].

At time *t* = *tk*−1/2, the semi-discrete scheme for the governing Equation (5) is given as follows:

$$\begin{split} &\frac{\partial}{\partial\sum\_{i=0}^{M\_1}\omega\_i^{(1)}\nabla\_t^{(\mathcal{G}\_i)}u^{k-1/2} + \frac{1}{2}\sum\_{i=0}^{M\_2}\omega\_i^{(2)}\nabla\_t^{(a\_i)}u^{k-1/2} \\ &=\frac{\partial}{\partial x}\left[d\_1(\mathbf{x},\mathbf{y})\_LD\_\mathbf{x}^\gamma u^{k-1/2} - d\_2(\mathbf{x},\mathbf{y})\_xD\_R^\gamma u^{k-1/2}\right] + \epsilon(\mathbf{x},\mathbf{y})\frac{\partial^2 u^{k-1/2}}{\partial y^2} + f^{k-1/2}.\end{split} \tag{17}$$

Define *V* = *Hγ* 0 (Ω) ∩ *H*1 0 (Ω) to be the numerical solution space. In this work, we choose to use triangular elements to mesh Ω. Because the domain is irregularly shaped, we refer to this throughout this paper as an unstructured mesh. Denote {Γ*h*} as a family of unstructured triangulations of domain Ω, where *h* is the maximum diameter of any triangle in Γ*h*. Then, we obtain the conforming finite element subspace *Vh* ∈ *V* as

$$V\_{\hbar} = \{ v\_{\hbar} | v\_{\hbar} \in \mathbb{C}(\overline{\Omega}) \cap V\_{\prime} v\_{\hbar} |\_{K} \text{ is linear for all } K \in \Gamma\_{\hbar} \text{ and } v\_{\hbar} |\_{\partial \Omega} = 0 \}.$$

Assume that *uk h* is the approximation of *u*(*x*, *y*, *t*) at time *t* = *tk*. We can derive the fully discrete formulation of (5)–(7): find *uk h* ∈ *Vh*, for any *k* = 0, 1, . . . , *N*, such that

$$\begin{split} &\mathbb{E}\sum\_{i=0}^{M\_{1}}\frac{\omega\_{i}^{(1)}\mathbbm{1}^{1-\beta\_{i}}}{\Gamma(3-\beta\_{i})}\left[a\_{0}^{\theta\_{i}}(\nabla\_{i}u\_{h}^{k-1/2},v\_{h})-\sum\_{j=1}^{k-1}(a\_{k-1-j}^{(\delta\_{i})}-a\_{k-j}^{(\delta\_{i})})(\nabla\_{i}u\_{h}^{j-1/2},v\_{h})-a\_{k-1}^{(\delta\_{i})}\left((u\_{h}^{0})\_{t},v\_{h}\right)\right] \\ &+\frac{1}{2}\sum\_{i=0}^{M\_{2}}\frac{\omega\_{i}^{(2)}\mathbbm{1}^{1-\kappa\_{i}}}{\Gamma(2-\kappa\_{i})}\sum\_{j=1}^{k}b\_{k-j}^{(a\_{j})}(\nabla\_{i}u\_{h}^{j-1/2},v\_{h})+\frac{1}{2}\sum\_{i=0}^{M\_{2}}\frac{\omega\_{i}^{(2)}\mathbbm{1}^{1-\kappa\_{i}}}{\Gamma(2-\kappa\_{i})}\sum\_{j=1}^{k-1}b\_{k-j-1}^{(a\_{i})}(\nabla\_{i}u\_{h}^{j-1/2},v\_{h}) \\ &=-\mathcal{B}(u\_{h}^{k-1/2},v\_{h})+(f^{k-1/2},v\_{h}), \end{split} \tag{18}$$

with the initial conditions and boundary conditions given by

$$
u\_h^0 =
u\_{0h\prime} \left.
u\_h^k \right|\_{\left(\partial \Omega\right)} = 0,\tag{19}$$

where *u*0*h* ∈ *Vh* is a reasonable approximation for *u*0. The expression for *B*(*u*, *v*) is given as

$$B(\boldsymbol{u}, \boldsymbol{v}) = \left(d\_1(\mathbf{x}, \boldsymbol{y})\_L D\_\mathbf{x}^\gamma \boldsymbol{u}, \frac{\partial \boldsymbol{v}}{\partial \mathbf{x}}\right) - \left(d\_2(\mathbf{x}, \boldsymbol{y})\_x D\_R^\gamma \boldsymbol{u}, \frac{\partial \boldsymbol{v}}{\partial \mathbf{x}}\right) + \left(\boldsymbol{e}(\mathbf{x}, \boldsymbol{y}) \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{y}'}, \frac{\partial \boldsymbol{v}}{\partial \boldsymbol{y}}\right). \tag{20}$$

#### *3.2. Implementation of Finite Element Method with an Unstructured Mesh*

In this section, we provide details of the implementation of the finite element method with an unstructured mesh. Firstly, we use the software Gmsh [32] to partition the convex domain Ω with unstructured triangular elements. For every triangular element *ep*, define *Ne* as total number of the triangles and *NP* is the number of elements. Using piecewise linear polynomials on every triangular element *ep*, for each time step, we can write *uk h* in the form *uk h* = *NP* ∑ *n*=1 *uk nϕn*(*x*, *y*), where *ϕn*(*x*, *y*) is the basis function and *uk n* is the unknown to be solved for. Denote *ϕn*(*xm*, *ym*) = *δnm* (*n*, *m* = 1, 2, ... , *Np*), where *δnm* refers to the Kronecker delta function, *vh* = *ϕm*(*x*, *y*), the mass matrix *M* = (*ϕn*, *ϕm*)*Np*×*Np* , the stiffness matrix *A* = *B*(*ϕn*, *ϕm*)*Np*×*Np* , *Fk* = (*Fk* 1 , *Fk* 2 ,..., *Fk Np* )*T* where *Fk m* = *f k*+*f k*−1 2 , *ϕm* ,

*ω*0 *m* = (*u*0 *h*)*t* ,*ϕm* , *Uk* = # *uk* 1, *uk* 2,..., *uk Np* \$*T* , *W*0 = # *ω*0 1, *ω*0 2,..., *ω*0 *Np* \$*T* , then we can rewrite (18) in matrix forms as follows

$$\begin{split} & \left[ \left( 2\xi \sum\_{i=0}^{M\_1} \omega\_i u\_0^{(\beta\_i)} + 2 \sum\_{i=0}^{M\_2} r\_i b\_0^{(a\_i)} \right) M + \tau A \right] \mathbb{I} l^k \\ &= \left[ \left( 2\xi \sum\_{i=0}^{M\_1} \omega\_i u\_0^{(\beta\_i)} + 2 \sum\_{i=0}^{M\_2} r\_i b\_0^{(a\_i)} \right) M + \tau A \right] \mathbb{I} l^{k-1} + 2\xi \tau \sum\_{i=0}^{M\_1} \omega\_i u\_{k-1}^{(\beta\_i)} \mathbb{W}^0 + 2\tau \mathbb{I} r^k \\ &+ 2\xi \sum\_{i=0}^{M\_1} \omega\_i \sum\_{j=1}^{k-1} \left( a\_{k-1-j}^{(\beta\_i)} - a\_{k-j}^{(\beta\_i)} \right) \left[ M \mathbb{I} l^j - M \mathbb{I} l^{j-1} \right] - 2 \sum\_{i=0}^{M\_2} r\_i \sum\_{j=1}^{k-1} \left( b\_{k-j}^{(a\_i)} + b\_{k-j-1}^{(a\_i)} \right) \left[ M \mathbb{I} l^j - M \mathbb{I} l^{-1} \right], \end{split}$$

where *ωi* = *ω*(1) *i τ*1−*βi* Γ(3−*βi*) , *ri* = *ω*(2) *i τ*1−*αi* 2Γ(2−*αi*) .

The critical point to obtain the solution is to approximate the first and the second terms of *B*(*ϕn*, *ϕm*). By applying the Gauss quadrature [23], we obtain

 *d*1(*x*, *y*)*LDγ x ϕn*, *∂ϕm ∂x* = ∑ *E*∈Γ*h* - *E d*1(*x*, *y*)*LDγ x ϕn ∂ϕm ∂x dxdy* ≈ ∑ *E*∈Γ*h* ∑ (*xci*,*yci*)∈*GE λid*1(*x*, *y*)*LDγ x ϕn*| (*xci*,*yci*) *∂ϕm ∂x* (*xci*,*yci*) , *d*2(*x*, *y*)*RDγ x ϕn*, *∂ϕm ∂x* = ∑ *E*∈Γ*h* - *E d*2(*x*, *y*)*xDγ Rϕn ∂ϕm ∂x dxdy* ≈ ∑ *E*∈Γ*h* ∑ (*xci*,*yci*)∈*GE κid*2(*x*, *y*)*xDγ Rϕn*| (*xci*,*yci*) *∂ϕm ∂x* (*xci*,*yci*) ,

where *GE* is the set of Gauss points in a certain element *E* and *λi*, *κi* are the weight coefficients corresponding to the Gauss points (*xci*, *yci*). In this article, we used four Gauss points in each triangle.

**Remark 1.** *I would be precise on the number of Gauss points used and the accuracy of the approximation used.*

In this paper, we use four Gauss points in every triangle. The detailed computation process can be summarized in Algorithm 1.

**Algorithm 1** Calculate (*d*1(*x*, *y*)*LDγ x ϕn*, *∂ϕm ∂x* ) and (*d*2(*x*, *y*)*RDγ x ϕn*, *∂ϕm ∂x* ) using finite element method on an unstructured mesh

#### **4. Stability and Convergence**

In this section, we analyze the stability and convergence of the discrete scheme on the irregular convex domain. In the following section, for the sake of simplicity, we consider the constant diffusion coefficient case *d*1(*x*, *y*) = *d*2(*x*, *y*) = *e*(*x*, *y*) = 1. Denote (·, ·)=(·, ·)Ω, || · ||0 = || · ||*L*2(Ω). Prior to presenting the numerical analysis, we firstly give the definitions of the semi-norm |·|(*γ*,1) and norm || · ||(*γ*,1)

$$\left|u\right|\_{\left(\gamma,1\right)} = \left(\left|\left|\_{L}D\_{\mathbf{x}}^{\left(\gamma+1\right)/2}u\right|\right|\_{0}^{2} + \left|\left|\left|\frac{\partial u}{\partial y}\right|\right|\_{0}^{2}\right)^{1/2}\right.\tag{21}$$

$$||u||\_{\left(\gamma,1\right)} = \left(||u||\_0^2 + |u|\_{\left(\gamma,1\right)}^2\right)^{1/2}.\tag{22}$$

In what follows, the constant *C* may be different in various sections.

#### *4.1. Stability*

**Lemma 7.** *For any u* ∈ *V, the semi-norm* |*u*|(*γ*,1) *and norm* ||*u*||(*γ*,1) *are equivalent and there exists positive constants C*1 *and C*2 *independent of u, such that*

$$C\_1 ||u||\_{\left(\gamma,1\right)} \le |u|\_{\left(\gamma,1\right)} \le ||u||\_{\left(\gamma,1\right)} \le C\_2 |u|\_{H^1\left(\Omega\right)}.$$

**Proof.** From Lemma 2, we immediately have

$$||u||\_0 \le C||\_L D\_x^{(\gamma+1)/2} u||\_0 \le C|u|\_{(\gamma,1)}.$$

By applying Lemma 2 and the Definitions (21) and (22), we have

$$||u||\_{\left(\gamma,1\right)} = \mathbb{C}\_{2} (||u||\_{0}^{2} + |u|\_{\left(\gamma,1\right)}^{2})^{1/2} \leq \mathbb{C}\_{2} |u|\_{\left(\gamma,1\right)} \leq \mathbb{C} |u|\_{\tilde{l}\_{L}^{\gamma,1}\left(\Omega\right)} \leq \mathbb{C}\_{2} |u|\_{H^{1}\left(\Omega\right)}.$$

By using Lemma 4 and applying the definitions of the norm and semi-norm, we have

$$||u||\_{\left(\gamma,1\right)} = \left(||u||\_0^2 + |u|\_{\left(\gamma,1\right)}^2\right)^{1/2} \ge \left(|u|\_{\left(\gamma,1\right)}^2\right)^{1/2} = |u|\_{\left(\gamma,1\right)}.\tag{23}$$

The proof is completed.

**Lemma 8.** *For any u*, *v* ∈ *Hγ* 0 (Ω) ∩ *H*1 0 (Ω)*, there exists constants C*1 *and C*2 *such that the function B*(*u*, *v*) *satisfies* |*B*(*u*, *v*)| ≤ *C*1||*u*||(*γ*,1)||*v*||(*γ*,1) *and B*(*u*, *u*) ≥ *C*2||*u*||2 (*γ*,1) *.*

**Proof.** Firstly, by using Lemma 1, Lemma 4 and Lemma 7, the Definitions (21) and (22) and applying the Cauchy–Schwartz inequality, namely (*u*, *v*) ≤ ||*u*||0||*v*||0, we have

|*B*(*u*, *v*)| ≤ (*LD*(*γ*+1)/2 *x u*,*x D*(*γ*+1)/2 *R v*) + (*xD*(*γ*+1)/2 *R u*,*L D*(*γ*+1)/2 *x v*) + *∂u ∂y* , *∂v ∂y* ≤ *LD*(*γ*+1)/2 *x u* 0 *xD*(*γ*+1)/2 *R v* 0 + *xD*(*γ*+1)/2 *R u* 0 *LD*(*γ*+1)/2 *x v* 0 + *∂u ∂y* 0 *∂v ∂y* 0 ≤*C* |*u*|(*γ*,1)|*v*|(*γ*,1) + |*u*|(*γ*,1)|*v*|(*γ*,1) + |*u*|(*γ*,1)|*v*|(*γ*,1) ≤ *C*||*u*||(*γ*,1)||*v*||(*γ*,1), *B*(*u*, *u*) ≥ (*LD*(*γ*+1)/2 *x u*, *xD*(*γ*+1)/2 *R u*) + (*xD*(*γ*+1)/2 *R u*, *LD*(*γ*+1)/2 *x u*) + *∂u ∂y* , *∂u ∂y* ≥ *C* (*LD*(*γ*+1)/2 *x u*, *xD*(*γ*+1)/2 *R u*) + *∂u ∂y* , *∂u ∂y* ≥ *C*|*u*| 2 *H*1(Ω) ≥ *C*||*u*||2 (*γ*,1).

**Theorem 1.** *(Stability) The fully discrete scheme (18) is unconditionally stable and it holds that*

$$||u\_{\mathbf{h}}^{N}||\_{\left(\gamma,1\right)} \leq \mathbb{C} \left[ \max\_{1 \leq k \leq N} ||f^{k-1/2}||\_{0}^{2} + ||\phi\_{1}||\_{1}^{2} + ||\phi\_{0}||\_{\left(\gamma,1\right)} \right].$$

**Proof.** Denote *vh* = ∇*tuk*−1/2 *h* . By multiplying 2*τ* with each item and summing *k* from 1 to *N*, the discrete scheme (18) changes as

$$2\pi\xi\sum\_{i=0}^{M\_1}\omega\_i^{(1)}\sum\_{k=1}^N \left(\nabla\_t^{(\beta)}u\_h^{k-1/2}, \nabla\_t u\_h^{k-1/2}\right) + 2\pi\sum\_{k=1}^N \sum\_{i=0}^{M\_2} \omega\_i^{(2)} \left(\nabla\_t^{(a)}u\_h^{k-1/2}, \nabla\_t u\_h^{k-1/2}\right)$$

$$-2\pi\sum\_{k=1}^N B\left(u\_h^{k-1/2}, \nabla\_t u\_h^{k-1/2}\right) - 2\pi\sum\_{k=1}^N \left(f^{k-1/2}, \nabla\_t u\_h^{k-1/2}\right) = 0. \tag{24}$$

For the first term with the initial condition (*u*0 *h*)*t* = *φ*1, by using Lemma 5, we have

$$2\pi \mathfrak{F} \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \sum\_{k=1}^N \left(\nabla\_t^{(\mathfrak{f}\_i)} u\_h^{k-1/2} \, , \nabla\_t u\_h^{k-1/2} \right) \ge 2\pi \mathfrak{F} a\_1 \sum\_{k=1}^N \left| \left| \nabla\_t u\_h^{k-1/2} \right| \right|\_0^2 - \left| \mathfrak{F} a\_2 \right| \left| \phi\_1 \right| \Big|\_{0^\sigma}^2 \tag{25}$$

where *a*1 = *M*1 ∑ *i*=0 *ω*(1) *i T*1−*βi* 2Γ(2−*βi*) , *a*2 = *M*1 ∑ *i*=0 *ω*(1) *i T*2−*βi* Γ(3−*βi*) .

By applying Lemma 6, we have

$$2\pi \sum\_{k=1}^{N} \sum\_{i=0}^{M\_2} \omega\_i^{(2)} \left(\nabla\_t^{(a\_i)} u\_h^{k-1/2}, \nabla\_t u\_h^{k-1/2}\right) \ge 0. \tag{26}$$

Define a symmetrical and continuous function

$$B\_0(\boldsymbol{u}, \boldsymbol{v}) = \left| \left( {\_1D\_x^{(\gamma+1)/2} \boldsymbol{u}\_{\times x} D\_R^{(\gamma+1)/2} \boldsymbol{v} } \right) \right| + \left| \left( {\_1D\_R^{(\gamma+1)/2} \boldsymbol{u}\_{\times L} D\_x^{(\gamma+1)/2} \boldsymbol{v} } \right) \right| + \left| \left( \frac{\partial \boldsymbol{u}}{\partial \boldsymbol{y}'} \frac{\partial \boldsymbol{v}}{\partial \boldsymbol{y}} \right) \right|. \tag{27}$$

then we have *B*(*u*, *v*) ≤ *CB*0(*u*, *v*). For the newly defined function *B*0 *uk*−1/2 *h* , ∇*tuk*−1/2 *h* , we have *B*0 *uk*−1/2 *h* , ∇*tuk*−1/2 *h* = 1 2*τ B*0 *uk h*, *uk h* − *B*0 *uk*−1 *h* , *uk*−1 *h* [21,22]. Perform the summation of *k* from 1 to *N*, and we derive the following inequality

$$2\pi \sum\_{k=1}^{N} B\left(u\_h^{k-1/2}, \nabla\_t u\_h^{k-1/2}\right) \le 2\pi \mathbb{C} \sum\_{k=1}^{N} B\_0\left(u\_h^{k-1/2}, \nabla\_t u\_h^{k-1/2}\right)$$

$$= 2\pi \mathbb{C} \sum\_{k=1}^{N} \frac{1}{2\pi} \left[B\left(u\_h^k, u\_h^k\right) - B\left(u\_h^{k-1}, u\_h^{k-1}\right)\right] = \mathbb{C}\left(B\left(u\_h^N, u\_h^N\right) - B\left(u\_h^0, u\_h^0\right)\right). \tag{28}$$

By using the important inequality 2*ab* ≤ *a*2 2*ε* + 2*εb*2, the fifth item changes as

$$\begin{split} \sum\_{k=1}^{N} 2\tau \langle f^{k-1/2}, \nabla t u\_h^{k-1/2} \rangle &\leq \sum\_{k=1}^{N} \tau \left[ \frac{||f^{k-1/2}||\_0^2}{2a\_1 \tilde{\xi}} + 2a\_1 \tilde{\xi} ||\nabla t u\_h^{k-1/2}||\_0^2 \right] \\ &\leq \frac{T}{2a\_1 \tilde{\xi}} \max\_{1 \leq k \leq N} ||f^{k-1/2}||\_0^2 + 2a\_1 \tilde{\xi} \sum\_{k=1}^{N} \tau ||\nabla\_t u\_h^{k-1/2}||\_0^2. \end{split} \tag{29}$$

Then, by using the inequality of (25)–(29), Equation (24) changes as

$$B(\boldsymbol{u}\_h^N, \boldsymbol{u}\_h^N) \le \mathbb{C} \left[ \frac{T}{2a\_1 \mathfrak{F}} \max\_{1 \le k \le N} ||f^{k-1/2}||\_0^2 + a\_2 \mathfrak{F} ||\phi\_1||\_0^2 + B(\boldsymbol{u}\_{h'}^0, \boldsymbol{u}\_h^0) \right].$$

Applying Lemma 8, we have

$$\|\|u\_{\hbar}^{N}||\_{\left(\gamma,1\right)}^{2} \leq C \left[ \max\_{1 \leq k \leq N} ||f^{k-1/2}||\_{0}^{2} + ||\phi\_{1}||\_{0}^{2} + ||\phi\_{0}||\_{\left(\gamma,1\right)} \right].$$

Therefore, the scheme is unconditionally stable.

#### *4.2. Convergence*

Prior to providing the convergence of the discrete scheme, we first give an approximation property. Define the interpolation operator *Ih*: *Hs*+1(Ω) → *Vh*, for any *u* ∈ *Hμ*(Ω), 1 < *μ* ≤ *s* + 1, there exists a constant *C* depending only on Ω such that [28]

$$||\mu - I\_h \mu||\_{H^1(\Omega)} \le Ch^{\mu - 1} ||\mu||\_{H^\mu(\Omega)}.\tag{30}$$

For any *u* ∈ *V* and *vh* ∈ *Vh*, we define a projection operator *Ph*: *V* → *Vh* possessing the following property

$$B(P\_{\hbar}u\_{\prime}v\_{\hbar}) = B(u\_{\prime}v\_{\hbar}).\tag{31}$$

Then, we have the following lemma.

**Lemma 9.** *If u* ∈ *Hμ*(Ω) ∩ *V,* 1 < *μ* ≤ *s* + 1*, there exists a constant C independent of h and u such that*

$$||P\_{\hbar}u - u||\_{\left(\gamma, 1\right)} \leq C\hbar^{\mu - 1}||u||\_{\mu}.\tag{32}$$

**Proof.** Because

$$||P\_h\mu - \mu||\_{\left(\gamma, 1\right)}^2 \le \mathbb{C}B\left(P\_h\mu - \mu, P\_h\mu - \mu\right) \le \mathbb{C}B\left(P\_h\mu - \mu, I\_h\mu - \mu\right),$$

and

$$B(P\_\hbar \mu - \mu, I\_\hbar \mu - \mu) \le C||P\_\hbar \mu - \mu||\_{\left(\gamma, 1\right)}||I\_\hbar \mu - \mu||\_{\left(\gamma, 1\right)}.$$

Using the approximation properties, we have

$$||P\_\hbar \mu - \mu||\_{\left(\gamma, 1\right)} \le C||I\_\hbar \mu - \mu||\_{\left(\gamma, 1\right)} \le C||I\_\hbar \mu - \mu||\_{H^1\left(\Omega\right)} \le Ch^{\mu - 1}||\mu||\_{\left(\mu, 1\right)}$$

**Theorem 2.** *(Convergence) Assume that uN* = *u*(*x*, *y*, *tN*) *is the exact solution with u, ∂αi u ∂t αi* , *∂βi u ∂t βi* ∈ *L*∞(*Hμ*(Ω); 0, *T*)*,* 1 < *μ* ≤ *s* + 1*, then the numerical solution uN h satisfies*

$$\begin{split} \left|| \left| u\_{h}^{\mu} - u(t\_{n}) \right| \right|\_{\left( \gamma, 1 \right)}^{2} &\leq \mathcal{C} \pi^{2\min\{ |\mathcal{S} - \beta\_{l}|, |\mathcal{Q} - a\_{i}| \}} \\ &+ \mathcal{C} h^{2\left( \mu - 1 \right)} \left[ \left| \left| u^{N} \right| \right|\_{\mu}^{2} + ||\phi||\_{\mu}^{2} + ||\phi||\_{\mu}^{2} + \max\_{1 \leq k \leq N} \left| \left| \left| \frac{\partial^{\ell\_{l}} u^{k - 1/2}}{\partial t^{\ell\_{l}}} \right| \right| \right|\_{\mu}^{2} + \max\_{1 \leq k \leq N} \left| \left| \frac{\partial^{\ell\_{l}} u^{k - 1/2}}{\partial t^{\ell\_{l}}} \right| \right|\_{\mu}^{2} \right]. \end{split}$$

**Proof.** Let *en* = *un h* − *u*(*tn*), then the newly defined *en* satisfies

$$\begin{split} \xi \sum\_{i=0}^{M\_1} \omega\_i^{(1)} (\nabla\_i^{(\delta)} \boldsymbol{\varepsilon}^{k-1/2}, \boldsymbol{v}\_h) + \xi \sum\_{i=0}^{M\_1} \omega\_i^{(1)} (R\_0^{\delta, k}, \boldsymbol{v}\_h) + \sum\_{i=0}^{M\_2} \frac{\omega\_i^{(2)} \tau^{1-a\_i}}{2\Gamma(2-a\_i)} \sum\_{j=1}^k b\_{k-j}^{(a\_i)} (\nabla\_i \boldsymbol{\varepsilon}^{j-1/2}, \boldsymbol{v}\_h) \\ + \sum\_{i=0}^{M\_2} \frac{\omega\_i^{(2)} \tau^{1-a\_i}}{2\Gamma(2-a\_i)} \sum\_{j=1}^{k-1} b\_{k-j-1}^{(a\_i)} (\nabla\_i \boldsymbol{\varepsilon}^{j-1/2}, \boldsymbol{v}\_h) + \sum\_{i=0}^{M\_3} \omega\_i^{(2)} (R\_1^{a\_i k}, \boldsymbol{v}\_h) + \mathcal{B}(\boldsymbol{\varepsilon}^{k-1/2}, \boldsymbol{v}\_h) = 0. \end{split} \tag{33}$$

Define *en* = *ρn* + *θn*, where *ρn* = *Phu*(*tn*) − *u*(*tn*) and *θn* = *un h* − *Phu*(*tn*). Then, for any *vh*, by using the Definition (31), we have

$$B(\rho^{k-1/2}, v\_h) = B\left(P\_h u(t\_{k-1/2}) - u(t\_{k-1/2}), v\_h\right) = 0.$$

In addition, choosing the interpolations as the initial values of *u*0 *h* at time *t*0, i.e., *u*0 *h* = *Ihφ*0, we obtain

$$\mathcal{B}(\theta^0, \theta^0) \le \mathcal{C} \mathcal{B}\_0(\theta^0, \theta^0) \le \mathbb{C} \left[ ||u(t\_0) - u^0\_h||\_{\left(\gamma, 1\right)}^2 + ||P\_h u(t\_0) - u(t\_0)||\_{\left(\gamma, 1\right)}^2 \right] \le \mathcal{C} h^{2\left(\mu - 1\right)} ||\phi\_0||\_{\mu}^2.$$

Similarly, choosing (*u*0 *h*)*t* = *Ihφ*1, the norm ||(*θ*0 *h*)*t* ||2 0 satisfies the following relationship

$$\|\| (\theta\_h^0)\_t \|\|\_0^2 \le \mathbb{C} \| |(u(t\_0) - u\_h^0)\_t| \|\_{\left(\gamma, 1\right)}^2 + | |(P\_h u(t\_0) - u(t\_0))\_t| \|\_{\left(\gamma, 1\right)}^2 \le C h^{2(\mu - 1)} ||\phi\_1||\_{\mu}^2.$$

Defining *vh* = ∇*tθk*−1/2 and summing *k* from 1 to *N*, Equation (33) can be rewritten as

$$\begin{split} & \quad \xi \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \sum\_{k=1}^N \left( \nabla\_i^{(\delta\_1)} \theta^{k-1/2}, \nabla\_i \theta^{k-1/2} \right) + \sum\_{k=1}^N \sum\_{i=0}^{M\_2} \omega\_i^{(2)} \left( \nabla\_i^{(a\_i)} \theta^{k-1/2}, \nabla\_i \theta^{k-1/2} \right) + \sum\_{k=1}^N \mathcal{B} \left( \theta^{k-1/2}, \nabla\_i \theta^{k-1/2} \right) \\ & = -\xi \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \sum\_{k=1}^N \left( \nabla\_i^{(\delta\_1)} \rho^{k-1/2}, \nabla\_i \theta^{k-1/2} \right) - \sum\_{k=1}^N \sum\_{i=0}^{M\_2} \omega\_i^{(2)} \left( \nabla\_i^{(a\_i)} \rho^{j-1/2}, \nabla\_i \theta^{k-1/2} \right) \\ & - \sum\_{k=1}^N \xi \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \left( \mathcal{R}\_0^{\delta,k}, \nabla\_i \theta^{k-1/2} \right) - \sum\_{k=1}^N \sum\_{i=0}^{M\_2} \omega\_i^{(2)} \left( \mathcal{R}\_1^{\delta,k}, \nabla\_i \theta^{k-1/2} \right). \end{split}$$

Note that ||·||0 ≤ ||·||(*γ*,1). By applying Lemma 9, the norm ||∇(*βi*) *t ρk*−1/2||2 0 can be estimated as

$$\begin{split} \left\lVert \left| \nabla\_{t}^{\langle \beta\_{i} \rangle} \rho^{k-1/2} \right| \right\rVert\_{0}^{2} &= \left\lVert \left| \nabla\_{t}^{\langle \beta\_{i} \rangle} \rho^{k-1/2} - \frac{\partial^{\beta\_{i}} \rho^{k-1/2}}{\partial t^{\beta\_{i}}} + \frac{\partial^{\beta\_{i}} \rho^{k-1/2}}{\partial t^{\beta\_{i}}} \right| \right\rVert\_{0}^{2} \\ &\leq \left\lVert \left| \nabla\_{t}^{\langle \beta\_{i} \rangle} \rho^{k-1/2} - \frac{\partial^{\beta\_{i}} \rho^{k-1/2}}{\partial t^{\beta\_{i}}} \right| \right\rVert\_{0}^{2} + \left\lVert \left| \frac{\partial^{\beta\_{i}} \rho^{k-1/2}}{\partial t^{\beta\_{i}}} \right| \right\rVert\_{0}^{2} \\ &= \left\lVert \left| R\_{0}^{\beta\_{i}k} \right| \right|\_{0}^{2} + \left\lVert \left| \frac{\partial^{\beta\_{i}}}{\partial t^{\beta\_{i}}} \right| P\_{h} u(t\_{k-1/2}) - u(t\_{k-1/2}) \right\rVert \right\rVert\_{0}^{2} \\ &\leq C \tau^{2(3-\beta\_{i})} + C h^{2(\mu-1)} \left\lVert \frac{\partial^{\beta\_{i}} u^{k-1/2}}{\partial t^{\beta\_{i}}} \right\rVert\_{\mu}^{2} . \end{split}$$

Similarly, we derive

$$\left| \left| \nabla\_t^{(a\_i)} \rho^{k-1/2} \right| \right|\_0^2 \le C r^{2(2-a\_i)} + C h^{2(\mu-1)} \left| \left| \frac{\partial^{x\_i} u^{k-1/2}}{\partial t^{a\_i}} \right| \right|\_{\mu}^2.$$

By applying the inequality in Lemma 9, the norms ||*ρN*||2 (*γ*,1) and ||*ρ*0||2 (*γ*,1) satisfy

$$\begin{aligned} ||\rho^N||\_{\left(\gamma,1\right)}^2 &= ||P\_\hbar u(t\_N) - u(t\_N)||\_{\left(\gamma,1\right)}^2 \leq \mathsf{C}h^{2\left(\mu-1\right)}||u^N||\_{\mu}^2, \\ ||\rho^0||\_{\left(\gamma,1\right)}^2 &= ||P\_\hbar u(t\_0) - u(t\_0)||\_{\left(\gamma,1\right)}^2 \leq \mathsf{C}h^{2\left(\mu-1\right)}||u^0||\_{\mu}^2. \end{aligned}$$

By using Lemma 5 and the initial condition (*θ*0 *h*)*t* = *φ*1, the following inequality holds

$$\left| \sqrt{\sum}\_{i=0}^{M\_1} \omega\_i^{(1)} \sum\_{k=1}^N \left( \nabla\_t^{(\mathfrak{f}\_i)} \theta^{k-1/2}, \nabla\_t \theta^{k-1/2} \right) \right| \ge a\_1 \overline{\xi} \sum\_{k=1}^N ||\nabla\_t \theta^{k-1/2}||\_0^2 - \frac{\overline{\xi} a\_2}{2\pi} ||(\theta\_h^0)\_t||\_0^2. \tag{34}$$

Applying Lemma 6, we have the following inequality

$$\sum\_{k=1}^{N} \sum\_{l=0}^{M\_2} \frac{\omega\_l^{(2)} \tau^{1-a\_l}}{2\Gamma(2-a\_l)} \sum\_{j=1}^{k} b\_{k-j}^{(a\_l)} (\nabla\_l \theta^{(-1/2)}, \nabla\_l \theta^{k-1/2}) + \sum\_{k=1}^{N} \sum\_{l=0}^{M\_2} \frac{\omega\_l^{(2)} \tau^{1-a\_l}}{2\Gamma(2-a\_l)} \sum\_{j=1}^{k-1} b\_{k-j-1}^{(a\_j)} (\nabla\_l \theta^{(-1/2)}, \nabla\_l \theta^{k-1/2}) \ge 0. \tag{35}$$

Applying the mid-point formula, we have the result

$$\sum\_{k=1}^{N} \mathcal{B} \left( \theta^{k-1/2}, \nabla\_i \theta^{k-1/2} \right) \ge \mathbb{C} \sum\_{k=1}^{N} \mathcal{B}\_0 \left( \theta^{k-1/2}, \nabla\_i \theta^{k-1/2} \right) \\ = \mathbb{C} \frac{1}{2\pi} \left[ \mathcal{B}\_0 \left( \theta^N, \theta^N \right) - \mathcal{B}\_0 \left( \theta^0, \theta^0 \right) \right]. \tag{36}$$

By using the important inequality −*ab* ≤ *εa*2 + *b*2/(4*ε*), we have

$$\begin{split} & -\mathfrak{F} \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \sum\_{k=1}^N \left( \nabla\_t^{(\beta\_i)} \rho^{k-1/2}, \nabla\_t \theta^{k-1/2} \right) \\ \leq & \xi \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \sum\_{k=1}^N \chi\_1 \left[ \mathsf{C} \tau^{2(3-\beta\_i)} + \mathsf{C} h^{2(\mu-1)} \left| \left| \frac{\partial^{\theta\_i} u^{k-1/2}}{\partial t^{\theta\_i}} \right| \right|\_{\mu}^2 \right] + \xi \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \sum\_{k=1}^N \frac{1}{4\chi\_1} ||\nabla\_t \theta^{k-1/2}||\_0^2 \\ \leq & \xi \mathsf{C} \left[ \tau^{2(3-\beta\_i)} + h^{2(\mu-1)} \max\_{1 \leq k \le N} \left| \left| \frac{\partial^{\theta\_i} u^{k-1/2}}{\partial t^{\theta\_i}} \right| \right|\_{\mu}^2 \right] + \xi \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \sum\_{k=1}^N \frac{1}{4\chi\_1} ||\nabla\_t \theta^{k-1/2}||\_0^2. \end{split} \tag{37}$$

$$\begin{split} & -\sum\_{k=1}^{N} \sum\_{i=0}^{M\_2} \omega\_i^{(2)} \left( \nabla\_t^{(a\_i)} \rho^{j-1/2}, \nabla\_t \theta^{k-1/2} \right) \\ & \leq \sum\_{i=0}^{M\_2} \omega\_i^{(2)} \sum\_{k=1}^{N} \chi\_2 \left( \mathcal{C} \tau^{2(2-a\_i)} + \text{Ch}^{2(\mu-1)} \left| \left| \frac{\partial^{a\_i} u^{k-1/2}}{\partial t^{a\_i}} \right| \right|\_{\mu}^2 \right) + \sum\_{k=1}^{N} \sum\_{i=0}^{M\_2} \frac{\omega\_i^{(2)}}{4\chi\_2} ||\nabla\_t \theta^{k-1/2}||\_0^2 \\ & \leq \text{NC} \left( \tau^{2(2-a\_i)} + h^{2(\mu-1)} \max\_{1 \leq k \leq N} ||\frac{\partial^{a\_i} u^{k-1/2}}{\partial t^{a\_i}}||\_{\mu}^2 \right) + \sum\_{k=1}^{N} \sum\_{i=0}^{M\_2} \frac{\omega\_i^{(2)}}{4\chi\_2} ||\nabla\_t \theta^{k-1/2}||\_0^2 \end{split} \tag{38}$$

where *χ*1 = *a*3 *a*1 , *χ*2 = *a*4 *ξa*1 , *a*3 = *M*1 ∑ *i*=0 *ω*(1) *i* and *a*4 = *M*2 ∑ *i*=0 *ω*(2) *i* . Similarly,

$$1 - \sum\_{k=1}^{N} \tilde{\xi} \sum\_{i=0}^{M\_1} \omega\_i^{(1)} (R\_0^{\otimes i, k}, \nabla\_t \theta^{k - 1/2}) \le N \mathbb{C} \tau^{2(3 - \beta\_i)} + \tilde{\xi} \sum\_{k=1}^{N} \sum\_{i=0}^{M\_1} \omega\_i^{(1)} \frac{1}{4\chi\_3} ||\nabla\_t \theta^{k - 1/2}||\_{0^+}^2 \tag{39}$$

$$-\sum\_{k=1}^{N} \sum\_{i=0}^{M\_2} (R\_1^{a\_i;k}, \nabla\_t \theta^{k-1/2}) \le N \mathbb{C} \tau^{2(2-a\_i)} + \sum\_{i=0}^{M\_2} \sum\_{k=1}^{N} \omega\_i^{(2)} \frac{1}{4\chi\_4} ||\nabla\_t \theta^{k-1/2}||\_{0^\prime}^2 \tag{40}$$

where *χ*3 = *a*3 *a*1 and *χ*4 = *a*4 *ξa*1 .

By using the inequalities (34)–(40), we obtain

$$\begin{split} &\frac{1}{2\tau}[B(\theta^{N},\theta^{N})-B(\theta^{0},\theta^{0})] \leq NC\Big{[}\tau^{2(3-\beta\_{i})}+h^{2(\mu-1)}\max\_{1\leq k\leq N}\left|\left|\frac{\partial^{\theta\_{i}}\mu^{k-1/2}}{\partial t^{\beta\_{i}}}\right|\Big|\_{\mu}^{2}\right] \\ &+NC\Big{(}\tau^{2(2-a\_{i})}+h^{2(\mu-1)}\max\_{1\leq k\leq N}\left|\left|\frac{\partial^{a\_{i}}\mu^{k-1/2}}{\partial t^{a\_{i}}}\right|\Big|\_{\mu}^{2}\right)+\frac{\tilde{\xi}^{a}a\_{2}}{2\tau}Ch^{2(\mu-1)}||\phi\_{1}||\_{\mu}^{2}. \end{split}$$

By utilizing Lemma 8, the above equation changes as

$$\mathbb{P}\left(|\theta^{N}|\right)^{2}\_{\left(\gamma,1\right)} \leq \mathsf{C}\tau^{2\min\left\{|\left|\boldsymbol{\beta}-\boldsymbol{\beta}\_{l}\right|,\left|\boldsymbol{2}-\boldsymbol{a}\_{l}\right|\right\}} + \mathsf{C}\mathcal{H}^{2(\mu-1)}\left[\max\_{1\leq k\leq N} \left|\left|\frac{\partial^{\boldsymbol{\beta}\_{l}}\boldsymbol{u}^{k-1/2}}{\partial\boldsymbol{t}^{\beta\_{l}}}\right|\right|\_{\mu}^{2} + \max\_{1\leq k\leq N} \left|\left|\frac{\partial^{\boldsymbol{\mu}\_{l}}\boldsymbol{u}^{k-1/2}}{\partial\boldsymbol{t}^{\mu\_{l}}}\right|\right|\_{\mu}^{2} + ||\boldsymbol{\phi}\_{l}||\_{\mu}^{2} + ||\boldsymbol{\phi}\_{l}||\_{\mu}^{2}\right].$$

The simplified form is given as

$$\begin{split} ||\boldsymbol{u}\_{h}^{N} - \boldsymbol{u}(t\_{N})||\_{\left(\gamma,1\right)}^{2} &\leq ||\boldsymbol{\rho}^{N}||\_{\left(\gamma,1\right)}^{2} + ||\boldsymbol{\theta}^{N}||\_{\left(\gamma,1\right)}^{2} \leq \mathsf{C}\tau^{2\min\left\{|\boldsymbol{\beta}-\boldsymbol{\beta}\_{i}|,|2-a\_{i}|\right\}} \\ &+ \mathsf{C}\mathsf{H}^{2\left(\boldsymbol{\mu}-1\right)} \left[ ||\boldsymbol{\mu}^{N}||\_{\boldsymbol{\mu}}^{2} + ||\boldsymbol{\phi}\_{0}||\_{\boldsymbol{\mu}}^{2} + ||\boldsymbol{\phi}\_{1}||\_{\boldsymbol{\mu}}^{2} + \max\_{1\leq k\leq N} ||\frac{\partial^{\boldsymbol{\theta}\_{i}}\mathsf{u}^{k-1/2}}{\partial\boldsymbol{\theta}^{\boldsymbol{\theta}\_{i}}}||\_{\boldsymbol{\mu}}^{2} + \max\_{1\leq k\leq N} ||\frac{\partial^{\boldsymbol{\theta}\_{i}}\mathsf{u}^{k-1/2}}{\partial\boldsymbol{\theta}^{\boldsymbol{\theta}\_{i}}}||\_{\boldsymbol{\mu}}^{2} \right]. \end{split}$$

The proof is completed.

**Remark 2.** *By using the triangular linear basis function, i.e., s* = 1*, it can be concluded from Theorem 2 that the error satisfies*

$$||u\_h^n - u(t\_n)||\_{\left(\gamma, 1\right)} \le C(r^{\min\{3-\beta i, 2-a\_i\}} + h).$$

#### **5. Numerical Examples**

In this section, we present two numerical examples: one is in a rectangular domain with the main purpose to demonstrate the effectiveness of our theoretical analysis, and the other is in an elliptical domain for analyzing the effects of different parameters on the particle distributions. In the mid-point quadrature rule [29], we choose *M*1 = 9, *M*2 = 9.

**Example 1.** *Firstly, we consider the following two-dimensional distributed-order time- and spacefractional diffusion-wave equation on a rectangular domain*

$$\begin{split} &\mathbb{E}\int\_{1}^{2}\varrho\_{1}(\boldsymbol{\beta})\frac{\partial^{\boldsymbol{\beta}}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial t^{\boldsymbol{\beta}}}d\boldsymbol{\beta}+\int\_{0}^{1}\varrho\_{0}(\boldsymbol{a})\frac{\partial^{\boldsymbol{a}}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial t^{\boldsymbol{a}}}d\boldsymbol{a} \\ &=\frac{\partial}{\partial\boldsymbol{x}}\Big[d\_{1}(\boldsymbol{x},\boldsymbol{y})\frac{\partial^{\boldsymbol{\gamma}}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial\boldsymbol{x}^{\boldsymbol{\gamma}}}-d\_{2}(\boldsymbol{x},\boldsymbol{y})\frac{\partial^{\boldsymbol{\gamma}}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial(-\boldsymbol{x})^{\boldsymbol{\gamma}}}\Big]+\frac{\partial}{\partial\boldsymbol{y}}\Big[\boldsymbol{\varepsilon}(\boldsymbol{x},\boldsymbol{y})\frac{\partial u(\boldsymbol{x},\boldsymbol{y},t)}{\partial\boldsymbol{y}}\Big]+f(\boldsymbol{x},\boldsymbol{y},t), \end{split}$$

subject to

$$\begin{aligned} u(\mathbf{x}, y, 0) &= \mathbf{x}^2 (1 - \mathbf{x})^2 y^2 (1 - y)^2, \; u\_t(\mathbf{x}, y, 0) = \mathbf{0}, \; (\mathbf{x}, y) \in \overline{\Omega}, \\ u(\mathbf{x}, y, \mathbf{t}) &= \mathbf{0}, \; (\mathbf{x}, y, \mathbf{t}) \in \partial \Omega \times [0, T], \end{aligned}$$

where Ω = (0, 1) × (0, 1). The exact solution of this problem is given by *u*(*x*, *y*, *t*) = (*t* 2 + 1)*x*2(1 − *x*)2*y*2(1 − *y*)2.

In Table 1, we take the special case with *ξ* = 1, *γ* = 0.8, *d*1(*x*, *y*) = *d*2(*x*, *y*) = *e*(*x*, *y*) = *x*2 + *y*2 to compute the *Hγ* error, *L*2 error and convergence order of *h* with *τ* = 1 1000 at *t* = 1 with different weight coefficients *ϕ*1(*β*) = *wi*(*β*) and *ϕ*0(*α*) = *ri*(*α*), *i* = 1, 2, 3, which are given by *w*1(*β*) = 0.5*δ*(*β* − 1.5) + 0.5*δ*(*β* − 1.8), *r*1(*α*) = 0.5*δ*(*α* − 0.5) + 0.5*δ*(*α* − 0.8), *w*2(*β*) = 2−*β* 2 , *r*2(*α*) = 1−*α* 2 , *w*3(*β*) = *β*2/2, *r*3(*α*) = *α*2/2. By examining the spatial convergence orders shown in Table 1, we notice that the expected convergence orders proved in Theorem 2 are obtained. With a different choice of the weight coefficient, the numerical solutions are in agreement with the theoretical analysis which indicates the validity of the proposed method.

**Example 2.** *In this example, we consider the following two-dimensional distributed-order timeand space-fractional diffusion-wave equation on an elliptical domain* Ω = {(*x*, *y*)| *x*2 *R*2 *a* + *y*2 *R*2 *b* < 1}

$$\begin{split} &\mathbb{E}\int\_{1}^{2}\varphi\_{1}(\boldsymbol{\beta})\frac{\partial^{\boldsymbol{\beta}}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial t^{\beta}}d\boldsymbol{\beta}+\int\_{0}^{1}\varphi\_{0}(\boldsymbol{a})\frac{\partial^{\boldsymbol{\alpha}}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial t^{\boldsymbol{\alpha}}}d\boldsymbol{a} \\ &=\frac{\partial}{\partial\boldsymbol{x}}\Big[d\_{1}(\boldsymbol{x},\boldsymbol{y})\frac{\partial^{\boldsymbol{\gamma}}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial\boldsymbol{x}^{\boldsymbol{\gamma}}}-d\_{2}(\boldsymbol{x},\boldsymbol{y})\frac{\partial^{\boldsymbol{\gamma}}u(\boldsymbol{x},\boldsymbol{y},t)}{\partial(-\boldsymbol{x})^{\boldsymbol{\gamma}}}\Big]+\frac{\partial}{\partial\boldsymbol{y}}\Big[\boldsymbol{\varepsilon}(\boldsymbol{x},\boldsymbol{y})\frac{\partial u(\boldsymbol{x},\boldsymbol{y},t)}{\partial\boldsymbol{y}}\Big]+f(\boldsymbol{x},\boldsymbol{y},t)\_{\boldsymbol{\gamma}}.\end{split}$$

**Table 1.** The *Hγ* error, *L*2 error and convergence order of *h* with *τ* = 1 1000 at *t* = 1 for the case *ξ* = 1, *γ* = 0.8, *d*1(*x*, *y*) = *d*2(*x*, *y*) = *e*(*x*, *y*) = *x*2 + *y*2 with different *ϕ*1(*β*) and *ϕ*0(*α*).

*subject to*

$$u(x,y,0) = \frac{1}{100} \left(\frac{x^2}{R\_a^2} + \frac{y^2}{R\_b^2} - 1\right)^2, \ u\_t(x,y,0) = 0, \quad (x,y) \in \overline{\Omega}\_r$$

$$u(x,y,t) = 0, \quad (x,y,t) \in \partial\Omega \times [0,T],$$

*where T* = 1. *The exact solution of this problem is given by u*(*x*, *y*, *t*) = *t* 2+1 100 *x*2 *R*2 *a* + *y*2 *R*2 *b* − 1 2 *.*

In the following discussions, for the sake of simplicity, all the numerical results listed in the tables and figures are evaluated at *Ra* = 0.5, *Rb* = 1. In Table 2, the *Hγ* error, *L*2 error and convergence order of *h* with *τ* = 1 1000 at *t* = 1 with different *ϕ*1(*β*) = *wi*(*β*), *ϕ*0(*α*) = *ri*(*α*), *i* = 1, 2 are presented, where *w*1(*β*) = *δ*(*β* − 1.8), *r*1(*α*) = *δ*(*α* − 0.8), *w*2(*β*) = 2−*β* 2 , *r*2(*α*) = 1−*α* 2 . The linear triangular elements are applied for this numerical example to verify the theoretical analysis. As we can see, the *Hγ* spatial convergence order is close to 1 while the *L*2 spatial convergence order is close to 2, which coincide with the theoretical analysis in Theorem 2. Through the above analysis, we see that our finite element algorithm also works well for an elliptical domain.

The solution behaviors with the effects of the different involved parameters, such as the relaxation parameter and the weight coefficient, are highlighted by graphical illustrations and analyzed in detail. We choose *ξ* = 1, *f*(*x*, *y*, *t*) = 0, *d*1(*x*, *y*) = 1−*x* 2 *δ*(*y*), *d*2(*x*, *y*) = 1+*x* 2 *δ*(*y*), *e*(*x*, *y*) = *x*2 + *y*2 + 1, *t* = 1 to observe the behaviors of the temporal evolution of the particle distribution, with the effect of the weight coefficients as shown in Figure 3. We use the exponential function with the form *δ*(*x*) ≈ 1 2 √*πσ e*− *x*2 4*σ* to approximate the Dirac delta function in the numerical simulation. Similar to [16], the weight coefficients are chosen as the power-law form with *ϕ*1(*β*) = *nβn*−1, *ϕ*0(*α*) = *nαn*−1 where *n* = 1/2, 1, 3. We observe that the impacts of the weight coefficients are significant on the solution behaviors. For *n* = 1/2, the weight coefficient is monotonically decreasing with the increase of *β* while monotonically increasing with the increase of *α*, and at this stage, the distribution presents as a diffusion form. For *n* = 1, the weight coefficient is constant which means that the weight for every fractional parameter is equal and the wave characteristic appears. With the increase of *n*, for *n* = 3, the weight coefficient is monotonically increasing with the increase of *β* and the decrease of *α*. As shown in Figure 3, the wave characteristic of the distribution becomes stronger.

**Table 2.** The *Hγ* error, *L*2 error and convergence order of *h* with *τ* = 1 1000 at *t* = 1 for the case *ξ* = 1, *γ* = 0.8, *d*1(*x*, *y*) = *d*2(*x*, *y*) = *e*(*x*, *y*) = *x*2 + *y*2 with different *ϕ*1(*β*) and *ϕ*0(*α*).

**Figure 3.** The three-dimensional and axial projection drawing particle distribution when *n* = 1/2, *n* = 1 and *n* = 3.

Figure 4 presents the influence of parameter *ξ* on the particle distributions for *t* = 1, *f*(*x*, *y*, *t*) = 0, *d*1(*x*, *y*) = 1−*x* 2 *δ*(*y*), *d*2(*x*, *y*) = 1+*x* 2 *δ*(*y*), *e*(*x*, *y*) = *x*2 + *y*2 + 1, *ϕ*1(*β*) = 1 and *ϕ*0(*α*) = 1. As the relaxation parameter increases from *ξ* = 0 to *ξ* = 1, the central region of the particle distribution begins to cave inward, which indicates that the distributions have a wave characteristic. The larger the relaxation parameter is, the stronger the wave characteristic will be. The reason is that the relaxation parameter is added on the distributedorder time-fractional derivative of order (1, 2), which possesses the wave characteristics. With an increase in the relaxation parameter, the fractional derivative of order (1, 2) with the wave characteristic plays a greater role in the particles' transport.

**Figure 4.** The three-dimensional and axial projection drawing particle distribution when *ξ* = 0, *ξ* = 1/2 and *ξ* = 1.

#### **6. Conclusions**

In this paper, we presented an original distributed-order time- and space-fractional diffusion-wave equation to analyze the anomalous diffusion in comb structures. The solution of the governing equation was obtained using the finite element method for the case where the coefficients are taken as constant. Two examples were given: one was in a rectangular domain and the other one was in an elliptical domain. In the two examples, the *Hγ* error, *L*2 error and convergence order of *h* with *τ* = 1 1000 at *t* = 1 subject to different weight coefficients showed that the results demonstrated the effectiveness of the numerical method. For the elliptical domain, the influence of the involved parameters, such as the relaxation parameter and the weight coefficient on the particle distribution, were analyzed, and the physical meaning of the diffusion-wave characteristics was discussed in detail.

**Author Contributions:** Conceptualization, L.L. and S.Z.; methodology, S.Z.; software, L.F.; validation, L.Z., J.Z. and F.L.; formal analysis, I.T.; investigation, S.Z.; resources, S.C.; data curation, L.L.; writing original draft preparation, S.Z.; writing—review and editing, L.Z.; visualization, S.C.; supervision, L.L.; project administration, L.L.; funding acquisition, L.L. All authors have read and agreed to the published version of the manuscript.

**Funding:** The work is supported by the Project funded by the National Natural Science Foundation of China (No. 11801029), the Fundamental Research Funds for the Central Universities (No. QNXM20220048) and the Open Fund of the State key laboratory of advanced metallurgy in the University of Science and Technology Beijing (No. K22-08) and the Australian Research Council (ARC) via the Discovery Project (DP180103858). Authors Liu and Zheng wish to acknowledge that this research is partially supported by the Natural Science Foundation of China (No. 11772046).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

**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.

### *Article* **Supervised Neural Network Procedures for the Novel Fractional Food Supply Model**

**Basma Souayeh 1,2,\*, Zulqurnain Sabir 3,4, Muhammad Umar 4 and Mir Waqas Alam 1,2**

	-<EMAIL_ADDRESS>
**Abstract:** This work presents the numerical performances of the fractional kind of food supply (FKFS) model. The fractional kinds of the derivatives have been used to acquire the accurate and realistic solutions of the FKFS model. The FKFSM system contains three types, special kind of the predator *L*(*x*), top-predator *M*(*x*) and prey populations *N*(*x*). The numerical solutions of three different cases of the FKFS model are provided through the stochastic procedures of the scaled conjugate gradient neural networks (SCGNNs). The data selection for the FKFS model is chosen as 82%, for training and 9% for both testing and authorization. The precision of the designed SCGNNs is provided through the achieved and Adam solutions. To rationality, competence, constancy, and correctness is approved by using the stochastic SCGNNs along with the simulations of the regression actions, mean square error, correlation performances, error histograms values and state transition measures.

**Keywords:** fractional order; food supply model; scaled conjugate gradient; artificial neural networks; numerical solutions; Adam method

### **1. Introduction**

There are various mathematical models that designate the natural phenomena based on the prey-predator investigations along with the collaborations of different species [1,2]. The functional response term in the prey-predator modelling has an important role to present that most of the prey affects the predators with the use of time. There are numerous functional responses species that have been reported in the literature, such as a ratio-dependent [3–5], Beddington–DeAngelis [6–8] and the Holling phase I to III [9,10]. One of the important models is food supply (FS), which is applied in the association of multiple prey or predators. The updated form of the FS system together with common qualitative investigations and numerous communications is presented in [11–13]. The mathematical modelling has an important role to present the dynamics of the nonlinear differential systems, e.g., SITR based coronavirus [14], dengue virus [15] and nervous stomach system [16].

In the FS chain, the role of the "Allee effects" is very important. The Allee effects defined in 1930 after the name by the famous scientist Allee. These effects allocate the progress to reduce the growing rate by using the small quantity of public. The Allee effects appear in the fishery, vertebrates, invertebrates, and plants. The Allee effects occasionally indicate the negative influences in the dispensation of population dynamics based on the fishery. The "Allee effects" have been divided into multiplication and addition [17–20]. Initially, Singh et al. described the double shape of "Allee effects" with the improved

**Citation:** Souayeh, B.; Sabir, Z.; Umar, M.; Alam, M.W. Supervised Neural Network Procedures for the Novel Fractional Food Supply Model. *Fractal Fract.* **2022**, *6*, 333. https:// doi.org/10.3390/fractalfract6060333

Academic Editors: Libo Feng, Lin Liu and Yang Liu

Received: 22 May 2022 Accepted: 14 June 2022 Published: 16 June 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Gower-Leslie system based on the prey predator, in which prey population shows the various junction associated with the suitable parameters. Vinoth et al. [21] formulated a mathematical model to investigate the dynamical FS system using the "Allee effect" based on the addition [22].

The aim of this work is to provide the numerical performances of the fractional kind of food supply (FKFS) model by using the stochastic procedures of the scaled conjugate gradient neural networks (SCGNNs). The stochastic solvers have been used to exploit the variety of applications in recent years. Few of them are the nonlinear dynamics of the coronavirus models [23], functional form of the singular models [24], the infectious-based HIV models [25], functional for of the delay differential system [26] and nonlinear model of the smoking [27]. The idea to implement the fraction kinds of the derivatives is to perform the accurate and realistic solutions. In fractional order models, the minute particulars based on the superfast transition and super slow evolution are examined that provides more detail of the dynamics of the system by using the fractional calculus, which is not easy to interpret by using the integer order counterparts. Additionally, the system dynamics for the index is performed by using the fractional calculus. The fractional order derivatives show much better performance as compared to the integer order with the availability of the situation. The fractional kind of the derivatives have been applied to authenticate the performance of the system using the applications of the real-world applications [28,29]. Moreover, the fractional derivatives have been extensively investigated to solve the number of applications based on the control networks, engineering, physical and mathematical systems. The implementation of the fractional calculus is performed broadly over the last 30 years by using the substantial operators, such as Weyl-Riesz [30], Caputo [31], Riemann-Liouville [32], Erdlyi-Kober [33], and Grnwald-Letnikov [34]. All these operators have their own worth and significance. However, the most widely definition of the Caputo derivative that works to solve homogeneous initial conditions as well as non-homogeneous initial conditions. The Caputo derivatives are considered easy to implement as compared to the other definitions. Bases on this fractional order applications, the authors are interested to develop the FKFS model and provide the numerical performances through the SCGNNs.

The remaining structure of the paper is given as: The FKFS system is constructed in the Section 2. The designed methodology based on the stochastic SCGNNs procedures is provided in the Section 3. The simulations of the results are provided in the Section 4. The concluding notes are given in the Section 5.

#### **2. Mathematical FKFS System with Insights**

In this section, the communication model is provided based on the two or more prey and predators. A differential FKFS system using the analysis of mutual qualitative along with the multiple relationships is given in [35,36]. Few researchers presented the multiple trophic-level of food supply systems through the structure of logistic prey *L*(*x*), Holling type or Lotka–Volterra predator *M*(*x*) and top-predator *N*(*x*) [37–43]. The mathematical form of the three species based on the FS system is presented as [44]:

$$\begin{cases} \frac{dL(\mathbf{x})}{dx} = a\_0 L(\mathbf{x}) - \frac{\rho\_0 L(\mathbf{x}) M(\mathbf{x})}{L(\mathbf{x}) + d\_0} - \frac{k\_1}{k\_2 + L(\mathbf{x})} - b\_0 L^2(\mathbf{x}), & L\_0 = i\_1 \\ \frac{dM(\mathbf{x})}{dx} = \frac{\rho\_1 L(\mathbf{x}) M(\mathbf{x})}{L(\mathbf{x}) + d\_1} - a\_1 M(\mathbf{x}) - \frac{\rho\_2 M(\mathbf{x}) N(\mathbf{x})}{M(\mathbf{x}) + d\_2}, & M\_0 = i\_2 \\ \frac{dN(\mathbf{x})}{dx} = c\_3 N^2(\mathbf{x}) - \frac{\rho\_3 N^2(\mathbf{x})}{M(\mathbf{x}) + d\_3} & M\_0 = i\_3 \end{cases} \tag{1}$$

where prey *L*(*x*) and species *M*(*x*) indicate the Volterra scheme that presents the population of the predator to decrease exponentially in the prey absence. The relationship of the species *N*(*x*) and the prey *M*(*x*) is provided by using the Leslie–Gower approach that represents the predator population reduces per capita accessibility [45,46]. *a*0 and *c*3 are the growth rates of *L*(*x*) and *N*(*x*), the environmental protection factors for *L*(*x*) are *d*0 and *d*1, while the reduction per capita of *M*(*x*) is *υ*2 2 described in *d*2, the term *a*1 shows the values of the *M*(*x*), which reduces the nonappearance of *L*(*x*), *b*0 provides the competition strength for *L*(*x*), the residual lessens for *N*(*x*) based on the food shortage *M*(*x*) is signified by *d*3, the maximum

presentations through the lessening of per capita of *L*(*x*) is represented by *ρ*0, *ρ*1, *ρ*2 and *ρ*3, the hyperbolic *k*1 *k*2+*L*(*x*) function shows the addictive form of the Allee effects, while *k*1 and *k*2 are the constant values of the Allee effects. If *k*1 < *k*2, then it means a weak Allee effect; otherwise *k*2 < *k*1, shows a strong Allee effect, the initial conditions are represented by *i*1, *i*2 and *i*3. The mathematical form of the FKFS system is given as:

$$\begin{cases} \frac{dL^v(\mathbf{x})}{dx^v} = a\_0 L(\mathbf{x}) - \frac{\rho\_0 L(\mathbf{x}) M(\mathbf{x})}{L(\mathbf{x}) + d\_0} - \frac{k\_1}{k\_2 + L(\mathbf{x})} - b\_0 L^2(\mathbf{x}), & L\_0 = i\_1, \\\frac{d^v M(\mathbf{x})}{dx^v} = \frac{\rho\_1 L(\mathbf{x}) M(\mathbf{x})}{L(\mathbf{x}) + d\_1} - a\_1 M(\mathbf{x}) - \frac{\rho\_2 M(\mathbf{x}) N(\mathbf{x})}{M(\mathbf{x}) + d\_2}, & M\_0 = i\_2. \\\frac{d^v N(\mathbf{x})}{dx^v} = c\_3 N^2(\mathbf{x}) - \frac{\rho\_3 N^2(\mathbf{x})}{M(\mathbf{x}) + d\_3}, & M\_0 = i\_3. \end{cases} \tag{2}$$

where *υ* shows the fractional order Caputo derivative to solve the fractional FS model given in Equation (2). The values of the fractional order derivative *υ* are taken between 0 and 1 to present the behavior of the fractional FS model. The fractional kinds of the derivative in the FS system (2) are incoroprated to observe the minute particulars, i.e., superslow evolution and superfast transients that is not easy to interpret by using the integer order counterparts as shown in the system (1). In recent few years, the fractional calculus have been implemented in various submission, such as anomalous heat transfer [47], pine wilt disease model with convex rate [48], patterns of the spatiotemporal using the systems based on the Belousov–Zhabotinskii reaction [49], quantitative approximation of soil animal substance content using the visible/near infrared spectrometry [50], predator-prey model with herd performance [51], Hepatitis B virus mathematical model [52] and biological based population growing model using the carrying volume [53].

The novel features of the proposed SCGNNs for solving the mathematical FKFS system are defined as:

#### **3. Designed SCGNNs Procedure**

This section of the study provides the procedure of the stochastic computing SCGNNs scheme for the mathematical form of the FKFS system as defined in the set of system (1). The workflow diagram is provided in Figure 1 for the mathematical FKFS model using the computing SCGNNs scheme based on the three blocks, the mathematical model, designed methodology and results performances. The design performances are given in two measures.

$$\begin{cases} \frac{d L^{\nu}(\infty)}{d \boldsymbol{\lambda}^{\nu}} = a\_{0} L(\infty) - \frac{\rho\_{0} L(\infty) M(\infty)}{L(\infty) + d\_{0}} - \frac{k\_{1}}{k\_{2} + L(\infty)} - b\_{0} L^{2}(\infty), & \boldsymbol{L}\_{0} = i\_{1}, \\\frac{d^{\nu} M(\infty)}{d \boldsymbol{\lambda}^{\nu}} = \frac{\rho\_{1} L(\infty) M(\infty)}{L(\infty) + d\_{1}} - a\_{1} M(\infty) - \frac{\rho\_{2} M(\infty) N(\infty)}{M(\infty) + d\_{2}}, & \boldsymbol{M}\_{0} = i\_{2}, \\\frac{d^{\nu} N(\infty)}{d \boldsymbol{\lambda}^{\nu}} = c\_{1} N^{2}(\infty) - \frac{\rho\_{3} N^{2}(\infty)}{M(\infty) + d\_{1}}, & \boldsymbol{M}\_{0} = i\_{1}. \end{cases}$$

**Figure 1.** Workflow of the SCGNNs construction to solve the FKFS model.

The significant procedures regarding to the generalization have been provided by using the Adam scheme, while the numerical procedures are implemented with the default parameter setting to generate the model dataset. The hidden neurons have been selected 15 in this study along with the data selection for the FKFS model as 82%, for training and 9% for both testing and authorization. The artificial intelligence abilities based supervised learning SCGNNs have been performed with best cooperation in the indices, including complexity, premature convergence, overfitting and underfitting cases. Additionally, these parameters of the networks are set after exhaustive simulation studies, experience, knowledge and care and small variations in these setting results in degraded performance of the networks.

The second phase of the stochastic SCGNNs is expressed by using the generic perception based on the solo neuron model as presented in Figure 2. The Figure 2a shows the single layered neural network structure, while the designed layer construction, a single input layer vector having 15 hidden numbers of neurons in the hidden layer along with the three outcomes in the outer layer as described in Figure 2b for solving the mathematical FKFS model. The stochastic based SCGNNs are applied by using the 'Matlab' software (nftool command) for the appropriate sections of hidden neurons, testing statistics, learning methods and verification statics. Whereas the implementation performances of the SCGNNs scheme to solve the mathematical FKFS model along with the parameter setting is provided in Table 1. The networks training is performed using the proposed stochastic SCGNNs scheme, where the backpropagation is oppressed to improve the Jacobian '*JB*' for the performance, i.e., MSE, to adjust the weight vectors along with the bias variables of *B*. The variation or modification of the decision variables with the use of scale conjugate gradient is given as:

$$\begin{array}{c} Jf = JB \times JB, \\ Je = JB \times e, \\ dB = \frac{-(JJ + I \times mn)}{fc} \end{array}$$

where *e* indicates the error, and *I* is the identity vector. The SCGNNs scheme's parameter setting is provided in Table 1 along with the slight disparity/change/modification may result in poor performance, i.e., premature convergence. Therefore, these settings will be unified with extensive attention, after directing thorough the numerical investigation and understanding.

**Table 1.** Parameter setting to execute the SCGNNs procedure.

#### **Table 1.** *Cont.*

$$\mathbf{(a)}$$

**Figure 2.** Generic and specific ANNs structure to solve the FKFS model. (**a**) Generic structure of single neuron; (**b**) Designed layer structure, a single input layer vector having 15 hidden numbers of neurons in the hidden layer along with the three outcomes in the outer layer.

#### **4. Results of the FKFS Model**

Three fractional order cases of the model have been presented by using the designed SCGNNs operator. The mathematical descriptions of these operators are given as:

**Case 1:** The updated form of Equation (2) based on the FKFS model by taking *υ* = 0.5, *a*0 = 1.5, *a*1 = 1, *b*0 = 0.06, *ρ*0 = 1, *ρ*1 = 2, *ρ*2 = 0.405, *ρ*3 = 1, *c*3 = 1.5, *k*1 = *k*2 = 0.1, *d*0 = 10, *d*1 = 10, *d*2 = 10, *d*3 = 20 and *i*1 = *i*2 = *i*3 = 1.2 is shown as:

$$\begin{cases} \frac{d^{15}L(\mathbf{x})}{dx^{0.5}} = 1.5L(\mathbf{x}) - \frac{L(\mathbf{x})M(\mathbf{x})}{L(\mathbf{x}) + 10} - \frac{0.1}{0.1 + L(\mathbf{x})} - 0.06L^2(\mathbf{x}), & L\_0 = 1.2, \\\frac{d^{05}M(\mathbf{x})}{dx^{0.5}} = \frac{2L(\mathbf{x})M(\mathbf{x})}{L(\mathbf{x}) + 10} - M(\mathbf{x}) - \frac{0.405M(\mathbf{x})N(\mathbf{x})}{M(\mathbf{x}) + 10}, & M\_0 = 1.2, \\\frac{d^{05}N(\mathbf{x})}{dx^{0.5}} = 1.5N^2(\mathbf{x}) - \frac{N^2(\mathbf{x})}{M(\mathbf{x}) + 20}, & N\_0 = 1.2. \end{cases} \tag{3}$$

**Case 2:** The updated form of Equation (2) based on the FKFS model by taking *υ* = 0.7, *a*0 = 1.5, *a*1 = 1, *b*0 = 0.06, *ρ*0 = 1, *ρ*1 = 2, *ρ*2 = 0.405, *ρ*3 = 1, *c*3 = 1.5, *k*1 = *k*2 = 0.1, *d*0 = 10, *d*1 = 10, *d*2 = 10, *d*3 = 20 and *i*1 = *i*2 = *i*3 = 1.2 is shown as:

$$\begin{cases} \frac{d^{0.7}L(\mathbf{x})}{dx^{0.7}} = 1.5L(\mathbf{x}) - \frac{L(\mathbf{x})M(\mathbf{x})}{L(\mathbf{x}) + 10} - \frac{0.1}{0.1 + L(\mathbf{x})} - 0.06L^2(\mathbf{x}), & L\_0 = 1.2, \\\frac{d^{0.7}M(\mathbf{x})}{dx^{0.7}} = \frac{2L(\mathbf{x})M(\mathbf{x})}{L(\mathbf{x}) + 10} - M(\mathbf{x}) - \frac{0.405M(\mathbf{x})N(\mathbf{x})}{M(\mathbf{x}) + 10}, & M\_0 = 1.2, \\\frac{d^{0.7}N(\mathbf{x})}{dx^{0.7}} = 1.5N^2(\mathbf{x}) - \frac{N^2(\mathbf{x})}{M(\mathbf{x}) + 20}, & N\_0 = 1.2. \end{cases} \tag{4}$$

**Case 3:** The updated form of Equation (2) based on the FKFS model by taking *υ* = 0.9, *a*0 = 1.5, *a*1 = 1, *b*0 = 0.06, *ρ*0 = 1, *ρ*1 = 2, *ρ*2 = 0.405, *ρ*3 = 1, *c*3 = 1.5, *k*1 = *k*2 = 0.1, *d*0 = 10, *d*1 = 10, *d*2 = 10, *d*3 = 20 and *i*1 = *i*2 = *i*3 = 1.2 is shown as:

$$\begin{cases} \frac{d^{0.5}L(\mathbf{x})}{dx^{0.9}} = 1.5L(\mathbf{x}) - \frac{L(\mathbf{x})M(\mathbf{x})}{L(\mathbf{x}) + 10} - \frac{0.1}{0.1 + L(\mathbf{x})} - 0.06L^2(\mathbf{x}), & L\_0 = 1.2, \\\frac{d^{0.9}M(\mathbf{x})}{dx^{0.9}} = \frac{2L(\mathbf{x})M(\mathbf{x})}{L(\mathbf{x}) + 10} - M(\mathbf{x}) - \frac{0.405M(\mathbf{x})N(\mathbf{x})}{M(\mathbf{x}) + 10}, & M\_0 = 1.2, \\\frac{d^{0.9}N(\mathbf{x})}{dx^{0.9}} = 1.5N^2(\mathbf{x}) - \frac{N^2(\mathbf{x})}{M(\mathbf{x}) + 20}, & N\_0 = 1.2. \end{cases} \tag{5}$$

Figures 3–7 illustrate the stochastic SCGNNs procedures for the FKFS mathematical system. Figure 3 shows the values of the STs along with the best performances of the FKFS mathematical system. The STs and MSE results based on the authentication, training and best curve measures have been demonstrated in Figure 3 using the stochastic SCGNNs procedures for the FKFS mathematical system. The obtained best measures of the FKFS model have been illustrated at iterations 81, 27 and 17 that have been performed as 7.58035 × 10−10, 1.72965 × 10 −9 and 4.49765 × 10−11. The second half of the Figure 3 shows the gradient values using the SCGNNs scheme for the FKFS mathematical system. The performances of the gradient are found as 9.35 × 10−8, 9.61 × 10−8 and 6.57 × 10−8. These depictions indicate the correctness and the convergence of the SCGNNs scheme for the FKFS mathematical system. The result assessments based on the training targets, training outputs, validations targets, validation outputs, test targets, test outputs, errors and fitness curves are illustrated in the 1st half of the Figure 4. While the EHs based on the training, validation, test and zero error have been drawn in the 2nd half of the Figure 4 for the FKFS mathematical system. The EHs performances are provided as 1.68 × 10−5, 5.79 × 10−6 and 1.05 × 10−7 for the FKFS mathematical system. Figure 5 represents the correlation performances based on the training, validation and testing in the mathematical form of the FKFS system. It is seen that the correlation measures are authenticated as 1 in the mathematical form of the FKFS system. These measures indicate the correctness of the stochastic SCGNNs procedure for the mathematical form of the FKFS model. The MSE convergence measures indicate the complexity values, training performances, validation measures, iterations, testing, and backpropagation are authenticated in Table 2 based on the mathematical form of the FKFS model.

**Table 2.** SCGNNs procedures for the mathematical form of the FKFS model.

**Figure 3.** MSE and STs for the mathematical form of the FKFS model. (**a**) MSE for C-1; (**b**) MSE for C-2; (**c**) MSE for C-3; (**d**) EHs for C-1; (**e**) EHs for C-2; (**f**) EHs for C-3.

**Figure 4.** Results valuations and EHs for the mathematical form of the FKFS model. (**a**) Result measures for C-1; (**b**) Result measures for C-2; (**c**) Result measures for C-3; (**d**) EHs for C-1; (**e**) EHs for C-2; (**f**) EHs for C-3.

**Figure 5.** Regression performances for the mathematical form of the FKFS model. (**a**) Regression for C-1; (**b**) Regression for C-2; (**c**) Regression for C-3.

**Figure 6.** Results overlapping for the mathematical form of the FKFS model. (**a**) Results of the logistic prey *L*(*x*); (**b**) Results of the Holling type *M*(*x*); (**c**) Results for the top-predator *N*(*x*).

Figures 6 and 7 indicate the comparative investigations based on the comparison of the solutions and AE performances to solve the FKFS system. Figure 6 shows the correctness of the SCGNNs scheme through the overlapping of the results for each class of the mathematical FKFS system. The AE values for each class of the mathematical FKFS system using the SCGNNs scheme are provided in Figure 7. The AE measures for the logistic prey *L*(*x*) are calculated as 10−5 to 10−6, 10−4 to 10−7 and 10−5 to 10−7 for 1st, 2nd and 3rd case of the mathematical FKFS system. The AE performances of the Holling type or Lotka–Volterra predator *M*(*x*) lie as 10−4 to 10−6, 10−4 to 10−5 and 10−5 to 10−8 for 1st, 2nd and 3rd case of the nonlinear FKFS system. The values of the AE for top-predator *N*(*x*) are calculated as 10−4 to 10−6, 10−5 to 10−6 and 10−5 to 10−8 for 1st, 2nd and 3rd case of the nonlinear FKFS system. These illustrations based on the AE authenticate the correctness of the stochastic SCGNNs LMB-NNs to solve the nonlinear FKFS system.

**Figure 7.** AE for the mathematical form of the FKFS model. (**a**) AE for the logistic prey *L*(*x*); (**b**) AE for the Holling type *M*(*x*); (**c**) AE for the top-predator *N*(*x*).

#### **5. Conclusions**

The motive of this work is to perform the solutions of the fractional food supply model. The fractional derivatives have been used to provide the realistic and accurate solutions of the food supply mathematical model. The fractional food supply mathematical system contains three categories, special kind of the predator *L*(*x*), top-predator *M*(*x*) and prey populations *N*(*x*). The efficient numerical performances of three different variations of the fractional food supply mathematical system have been provided by using the stochastic procedures based on the scaled conjugate gradient neural network scheme. The selection of the data for fractional food supply mathematical system is selected as 82%, for training and 9% for both testing and authorization along with the 15 numbers of neurons. The precision and accuracy of the designed SCGNNs have been provided through the achievements and reference solutions. The AE values have been calculated as 10−6 to 10−8, which shows the exactness of the scaled conjugate gradient neural network scheme for solving the fraction food supply system. The rationality, competence, constancy, and correctness has been approved by using the stochastic SCGNNs along with the simulations of the regression actions, mean square error, correlation performances, error histograms values and state transition measures. It is also observed that by taking the fractional order values close to 1, the solutions are performed better as compared to other values. These observations have been provided in the AE graphs to solve the model.

In upcoming studies, the proposed SCGNNs scheme have been implemented to present the solutions of the lonngren-wave systems, fluid dynamical models and fractional kinds of systems.

**Author Contributions:** Conceptualization, B.S. and Z.S.; methodology, B.S. and M.U.; software, B.S. and Z.S.; validation, B.S. and M.U.; formal analysis, Z.S. and M.W.A.; investigation, B.S. and Z.S.; resources, B.S., Z.S., M.U. and M.W.A.; data curation, B.S. and Z.S.; writing—original draft preparation, B.S. and M.W.A.; writing—review and editing, Z.S. and B.S.; visualization, M.U. and B.S.; supervision, B.S.; project administration, B.S.; funding acquisition, B.S.; All authors have read and agreed to the published version of the manuscript.

**Funding:** This work was supported by Al Bilad Bank Scholarly Chair for Food Security in Saudi Arabia, the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia [Project No. CHAIR37].

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data that support the findings of this study are available from the corresponding author upon reasonable request.

**Acknowledgments:** This work was supported by Al Bilad Bank Scholarly Chair for Food Security in Saudi Arabia, the Deanship of Scientific Research, Vice Presidency for Graduate Studies and Scientific Research, King Faisal University, Saudi Arabia [Project No. CHAIR37].

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

### *Article* **A Finite-State Stationary Process with Long-Range Dependence and Fractional Multinomial Distribution**

**Jeonghwa Lee**

Department of Statistics, Truman State University, Kirksville, MO 63501, USA<EMAIL_ADDRESS>
**Abstract:** We propose a discrete-time, finite-state stationary process that can possess long-range dependence. Among the interesting features of this process is that each state can have different long-term dependency, i.e., the indicator sequence can have a different Hurst index for different states. Furthermore, inter-arrival time for each state follows heavy tail distribution, with different states showing different tail behavior. A possible application of this process is to model overdispersed multinomial distribution. In particular, we define a fractional multinomial distribution from our model.

**Keywords:** long-range dependence; Hurst index; over-dispersed multinomial distribution

#### **1. Introduction**

Long-range dependence (LRD) refers to a phenomenon where correlation decays slowly with the time lag in a stationary process in a way that the correlation function is no longer summable. This phenomenon was first observed by Hurst [1,2] and since then it has been observed in many fields such as economics, hydrology, internet traffic, queueing networks, etc. [3–6]. In a second order stationary process, LRD can be measured by the Hurst index *H* [7,8],

$$H = \inf\{h : \limsup\_{n \to \infty} n^{-2h+1} \sum\_{k=1}^n cov(X\_1, X\_k) < \infty\}.$$

Note that *H* ∈ (0, 1), and if *H* ∈ (1/2, 1), the process possesses a long-memory property.

Among the well-known stochastic processes that are stationary and possess longrange dependence are fractional Gaussian noise (FGN) [9] and fractional autoregressive integrated moving average processes (FARIMA) [10,11].

Fractional Gaussian noise *Xj* is a mean-zero, stationary Gaussian process with covariance function:

$$\gamma(j) := \csc(X\_0, X\_j) = \frac{var(X\_0)}{2} (|j+1|^{2H} - 2|j|^{2H} + |j-1|^{2H})$$

where *H* ∈ (0, 1) is the Hurst parameter. The covariance function obeys the power law with exponent 2*H* − 2 for large lag,

$$
\gamma(j) \sim \text{var}(\mathbb{X}\_0) H(2H - 1) j^{2H - 2} \text{ as } j \to \infty.
$$

If *H* ∈ (1/2, 1), then the covariance function decreases slowly with the power law, and ∑*j γ*(*j*) = ∞, i.e., it has the long-memory property.

A FARIMA(*p*, *d*, *q*) process {*Xt*} is the solution of:

$$\phi(B)\nabla^d X\_t = \theta(B)\epsilon\_{t\star}$$

**Citation:** Lee, J. A Finite-State Stationary Process with Long-Range Dependence and Fractional Multinomial Distribution. *Fractal Fract.* **2022**, *6*, 596. https://doi.org/ 10.3390/fractalfract6100596

Academic Editors: Libo Feng, Yang Liu and Lin Liu

Received: 17 September 2022 Accepted: 11 October 2022 Published: 14 October 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2022 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where *p*, *q* are positive integers, *d* is real, *B* is the backward shift, *BXt* = *Xt*−1, and the fractional-differencing operator *d*, autoregressive operator *φ*, and moving average operator *θ* are, respectively,

$$\begin{aligned} \nabla^d &= (1 - B)^d = \sum\_{k=1}^{\infty} \frac{d(d-1)\cdots(d+1-k)}{k!}(-B)^k, \\ \phi(B) &= 1 - \phi\_1 B - \phi\_2 B^2 \cdots - \phi\_p B^p, \\ \theta(B) &= 1 - \theta\_1 B - \theta\_2 B^2 \cdots - \theta\_q B^q. \end{aligned}$$

where {*t*} is the white-noise process, which consists of iid random variables with the finite second moment. Here, the parameter *d* manages the long-term dependence structure, and by its relation to the Hurst index, *H* = *d* + 1/2, *d* ∈ (0, 1/2) corresponds to the long-range dependence in the FARIMA process.

Another class of stationary processes that can possess long-range dependence is from the countable-state Markov process [12]. In a stationary, positive recurrent, irreducible, aperiodic Markov chain, the indicator sequence of visits to a certain state is long-range dependent if and only if return time to the state has an infinite second moment, and this is possible only when the Markov chain has infinite state space. Moreover, if one state has the infinite second moment of return time, then all the other states also have the infinite second moment of return time, and all the states have the same rate of dependency; that is, the indicator sequence of each state is long-range dependence with the same Hurst index.

In this paper, we develop a discrete-time finite-state stationary process that can possess long-range dependence. We define a stationary process {*Xi*, *i* ∈ N} where the number of possible outcomes of *Xi* is finite, *S* = {0, 1, ··· , *m*} for any *m* ∈ N, and for *k* = 1, 2, ··· , *m*,

$$cov(I\_{\{X\_i = k\}}, I\_{\{X\_j = k\}}) = c'\_k |i - j|^{2H\_k - 2},\tag{1}$$

for any *i*, *j* ∈ N, *i* = *j*, and some constants *c k* ∈ R+, *Hk* ∈ (0, 1). This leads to:

$$
\cos(X\_{i\prime}X\_j) \sim c\_{k'}'|i-j|^{2H\_{k'}-2} \quad \text{as } |i-j| \to \infty,\tag{2}
$$

where *k* = *argmaxk*{*Hk*; *k* = 1, ··· , *m*}. If *Hk* = max{*Hk*; *k* = 1, ··· , *m*} ∈ (1/2, 1), (1.2) implies that as *n* → ∞, ∑*n i*=1 *cov*(*X*1, *Xi*) diverges with the rate of |*n*| 2*Hk* −1, and the process is said to have long-memory with Hurst parameter *Hk* . Furthermore, from (1.1), for *k* = {1, ··· , *m*}, the process {*I*{*Xi*=*k*}; *i* = 1, 2, ···} is long-range dependence if *Hk* ∈ (1/2, 1). In particular, if *Hi* = *Hj*, then the states "*i*" and "*j*" produce different levels of dependence. For example, if *Hi* < 1/2 < *Hj*, then the state "*j*" produces a long-memory counting process whereas state "*i*" produces a short-memory process.

A possible application of our stochastic process is to model the over-dispersed multinomial distribution. In the multinomial distribution, there are *n* trials, each trial results in one of the finite outcomes, and the outcomes of the trials are independent and identically distributed. When applying the multinomial model to real data, it is often observed that the variance is larger than what it is assumed to be, which is called over-dispersion, due to the violation of the assumption that trials are independent and have identical distribution [13,14], and there have been several ways to model an overdispersed multinomial distribution [15–18].

Our stochastic process provides a new method to model an over-dispersed multinomial distribution by introducing dependency among trials. In particular, the variance of the number of a certain outcomes among *n* trials is asymptotically proportional to the fractional exponent of *n*, from which we define:

$$\Upsilon\_k := \sum\_{i=1}^n I\_{\{X\_i = k\}} \text{ for } k = 1, 2, \dots, m\_n$$

and call the distribution of (*Y*1,*Y*2, ··· ,*Ym*) the fractional multinomial distribution.

The work in this paper is an extension of the earlier work of the generalized Bernoulli process [19], and the process in this paper is reduced to the generalized Bernoulli process if there are only two states in the possible outcomes of *Xi*, e.g., *S* = {0, 1}.