# Application of Intelligent Paradigm through Neural Networks for Numerical Solution of Multiorder Fractional Differential Equations

Authors: Naveed Ahmad Khan; Osamah Ibrahim Khalaf; Carlos Andrés Tavera Romero; Muhammad Sulaiman; Maharani A. Bakar
Journal: Computational Intelligence and Neuroscience
 (2022)

---

## Abstract

In this study, the intelligent computational strength of neural networks (NNs) based on the backpropagated Levenberg-Marquardt (BLM) algorithm is utilized to investigate the numerical solution of nonlinear multiorder fractional differential equations (FDEs). The reference data set for the design of the BLM-NN algorithm for different examples of FDEs are generated by using the exact solutions. To obtain the numerical solutions, multiple operations based on training, validation, and testing on the reference data set are carried out by the design scheme for various orders of FDEs. The approximate solutions by the BLM-NN algorithm are compared with analytical solutions and performance based on mean square error (MSE), error histogram (EH), regression, and curve fitting. This further validates the accuracy, robustness, and efficiency of the proposed algorithm.


---

## Body

## 1. Introduction

Mathematicians have regarded the theory of fractional calculus as a branch of pure mathematics for nearly three centuries. However, several researchers have recently discovered that noninteger derivatives and integrals are more useful for modelling phenomena with inherited and memory properties than integer orders [1–4]. Fractional differential equations (FDEs) are used to model various problems in science, engineering, economics, biological sciences, and applied mathematics [5–8]. FDEs are more complex than their integer order since the fractional operators are nonlocal and have weakly singular kernels [9–13]. The complications in integer order introduce significant computational difficulties for numerical methods to obtain solutions for such equations.Fractional differential equations have wide application in the fields of science and engineering. Some recent applications include fractional-order financial systems [14], electrical circuits [15], nuclear magnetic resonance [16], fractional-order Bloch system [17], fractional-order Lorenz system [18], hepatitis B disease in medicine [19], pollution levels in a lake [20], and fractional-order Chua’s system [21]. Due to the high usage of FDEs, several numerical and analytical methods have been proposed. Bhrawy [22–24] uses spectral methods based on Jacobi, Chebyshev, and Legendre polynomials over a bounded domain for an approximate solution of FDE’s. Atabakzadeh [25] and Tripathi [26] use the operational matrix of Caputo fractional-order derivatives for Chebyshev polynomials and fractional integration of the generalized hat basis functions to solve systems of FDEs. Baleanu [27] in 2013 uses modified generalized Laguerre collocation methods and the Tau method based on semi-infinite interval to calculate the approximate solution for linear and nonlinear FDEs. Ahmadian [28, 29] applied the Jacobi operational matrix to study a class of linear fuzzy FDEs. The spectral approximation method is used by Li [30] to compute the fractional derivative and integral and also presents the pseudo-spectral approximation technique for some classes of FDEs. Esmaeili [31] developed a numerical technique in which the properties of the Caputo derivative were used to reduce the fractional differential equation into a Volterra integral equation.Recently, the use of spectral methods to solve various types of differential and integral equations has gained interest due to their wide applicability in both finite and infinite domains [32–35]. These methods include Galerkin [33], collocation [36, 37], Tau [38], and Petrov Galerkin [39] classes. Researchers have widely used the Homotopy perturbation method (HPM) [40, 41], Legendre wavelets method (LWM) [42, 43], fractional-order Laguerre and Jacobi Tau methods [44, 45], Chebyshev Tau method (CTM) [45], variational iteration method (VIM) [46], differential transform method (DTM) [40], Bernoulli wavelets method (BWM) [47, 48], and Adomian decomposition method (ADM) [49] for the numerical solution of fractional differential equations.In recent times, stochastic numerical techniques based on artificial intelligence have been developed to solve stiff nonlinear problems arising in various fields. Such stochastic computing techniques use artificial neural networks to model approximate solutions. These numerical solvers have wide applications in various fields including petroleum engineering [50], heat transfer [51], civil engineering [52], wire coating dynamics [53], and diabetic retinopathy classification [54]. The abovementioned techniques inspire the authors to explore and incorporate soft computing architectures as an alternative, precise, and feasible way for solving nonlinear multiorder fractional differential equations. The main purpose of this article is to obtain approximate solutions for FDEs using artificial neural networks based on the Levenberg–Marquard algorithm. Some highlighted features of the given study are illustrated as follows:(i)
Novel applications of neuroheuristic techniques based on backpropagated Levenberg–Marquardt neural networks (BLM-NNs) are presented to obtain numerical solutions for different classes of nonlinear multi-order fractional differential equations.(ii)
The processes of training, validation, and testing are carried out by generating a reference solution or data set by using an analytical solution for different cases and examples of FDEs.(iii)
The performance of the proposed scheme is incorporated by fitting the approximate solutions with the reference solution. The absolute error between the targeted data and approximate solutions illustrates the worth and accuracy of the BLM-NN algorithm.(iv)
Convergence analysis based on mean square errors of the objective function, regression analysis, and histogram plots are employed to study the complexity, robustness, and correctness of the design scheme.(v)
The advantage of the proposed design is that it does not require any initial parameter settings. It has simple and smooth implementation with exhaustive applicability and stability.


## 2. Solution Methodology

In the field of artificial intelligence (AI), supervised machine learning refers to a collection of algorithms that describe a predictive model based on data set with known outcomes. The method is learned through the uses of an efficient teaching algorithm, such as artificial neural networks, which use optimization procedures to minimise the error function. The infrastructure of the proposed BLM-NN algorithm is based on two fundamental steps. In the first step, a data set of 1201 points is generated by using an analytical solution from 0 to 6 with a 0.005 step size. In the next step, the Levenberg–Marquardt framework of fitting tool “nftool” from the neural network toolbox of MATLAB R2018a is used to approximate the solutions with75% training, 15% validation, and 15% testing. The suggested structure of the BLM-NN algorithm with 60 neurons is shown in Figure 1. A summary of the working procedure of the design scheme is presented through the flow chart in Figure 2.Figure 1
Structure of a supervised neural network.Figure 2
A complete overview of the working procedure of the BLM-NN algorithm.The performance of a design scheme is measured through the performance indicators in terms of mean square error (MSE) of fitness function, regressionR2, error histograms, and absolute errors (AE). The mathematical formulation of the MSE, R2, and AE is given as follows:(1)MSE=1k∑j=1kujt−u^jt2,R2=1−∑j=1ku^jt−u¯jt2∑j=1kujt−u¯jt2,and(2)AE=ujt−u^jt,j=1,2,…,k.Here, uj, u¯j, and u^j denote the reference, approximate, and mean of the solution at jth input and k is the number of mesh points. The desired value of MSE and AE for perfect fitting is equal to zero, while the value of R2 is one.


## 3. Numerical Experimentation

In order to illustrate the performance of the BLM-NN algorithm, we have considered various examples of nonlinear multiorder fractional differential equations. All calculations and evaluations for this research are performed on HP laptop EliteBook 840 G2 with an intel(R) Core (TM) i5-5300 CPU @ 2.30 GHz, 8.00 GB RAM, 64 bit operating in Microsoft Windows 10 Education edition, running the R2018a version of MATLAB.Example 1.
Consider the following nonlinear fractional differential equation [55]:(3)Dυux+u2x=Γυ+2x+6x3−vΓ4−υ+xυ+1+x32,0<υ≤2,with the following equation:(4)u0=u′0=0.
The exact solution of (3) is ux=xυ+1+x3. Four fractional orders are considered i.e., Case I υ=1.2, Case II υ=1.4, Case III υ=1.6, and Case IV υ=1.8.
In order to find approximate solutions for various orders of (3), the BLM-NN algorithm is executed using “nftool” in the MATLAB package. The performance and convergence of the mean square error (MSE) of the objective function are shown in Figure 3. It can be seen that the best validated performance for υ=1.2,1.4,1.6, and 1.8 are 1.0846e−10, 8.5718e−11, 9.7898e−11, and 1.7456e−10 which are attained at 1000 epoch. Table 1 demonstrates the approximate solution for each case of Example 1. Further, the fitting of approximate solutions with analytical solutions is plotted in Figure 4. The absolute errors between targeted data and obtained solutions for multiple orders of (3) are illustrated in Figures 4 and 5, respectively. The values of AE for each case lie around 10−5 to 10−6, 10−5 to 10−7, 10−5 to 10−7, and 10−5 to 10−8, respectively. Table 2 represents the measure of convergence for each testing, validation, training, gradient, mu, and complexity analysis in terms of time taken by the system to achieve the desired results. It can be seen that the values for the gradient for each case lie around 10−4 to 10−7, while the maximum time taken by the system is 5 seconds. The training state of operators during the process of optimization for Example 1 is shown in Figure 6. The accuracy and efficiency of the proposed algorithm is shown by the results of regression as dictated in Figure 7.Figure 3
Performance analysis of design schemes based on MSE for different cases of Example1. (a) Case I, (b) case II, (c) case III, and (d) case IV.
(a)(b)(c)(d)Table 1
Approximate solutions obtained by the proposed algorithm for different cases of multiorder fractional differential equations.
xExample 1Example 2Example 3Case ICase IICase IIICase IV00000020.50.342 6380.314 4650.289 9380.268 5870.041 666 671.875122220.333 333 331.51.55.815 0616.021 1786.244 7056.487 1141.1250.875212.594 7913.278 0314.062 8714.964 42.666 666 6702.523.132 0324.641 8726.455 3928.633 645.208 333 33−1.125338.211 5840.966 6144.398 6448.674 029−2.53.558.613 0163.094 1768.851 2876.247 6414.291 666 7−4.125485.112 1391.857 62100.758 3112.502 921.333 333 3−64.5118.481 9128.083 1141.054158.57730.375−8.1255159.493 2172.591 3190.6632215.597541.666 666 7−10.55.5208.915226.198 3250.5035284.683455.458 333 3−13.1256267.514 9289.716 2321.4856366.946772−16Figure 4
Comparison of the analytical solution with the approximate solution obtained by the BLM-NN algorithm for different orders of FDE as illustrated in Example1. (a) Case I, (b) case II, (c) case III, and (d) case IV.
(a)(b)(c)(d)Figure 5
Error histogram between target values and approximated values for multiple orders of equation (4). (a) Case I, (b) case II, (c) case III, and (d) case IV.
(a)(b)(c)(d)Table 2
Statistical analysis of performance measures including MSE, gradient, mu, number of iterations, and time taken by the system for calculating the results of Examples1, 2, and 3.
ExampleCaseHidden neuronsMean square errorGradientMuEpochsRegressionTime (s)TrainingValidationTesting1I607.50E-111.08E-101.06E-101.00E-071.00E-10100015II605.84E-118.57E-111.13E-101.00E-071.00E-10100015III607.64E-119.79E-111.15E-105.82E-041.00E-10100015IV601.26E-101.75E-101.66E-101.40E-041.00E-091000152605.18E-141.12E-131.61E-131.24E-081.00E-13890123605.52E-121.58E-132.53E-121.51E-061.00E-10100015Figure 6
Training state of design scheme for all cases of Example1. (a) Case I, (b) case II, (c) case III, and (d) case IV.
(a)(b)(c)(d)Figure 7
Regression analysis of cases I, II, III, and IV of Example1.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)(m)(n)(o)(p)Example 2.
Consider the following nonlinear multiterm nonhomogenous fractional differential equation as [56](5)D2.2ux+Dα2ux+Dα1ux+u3x=2Γ1.8x0.8+2Γ4−α2x3−α2+2Γ4−α1x3−α1+x333,subjected to the following equation:(6)y0=y′0=y′′0=0,where α1 and α2 are 0.75 and 1.25, respectively. The exact solution of (5) is ux=x3/3. The approximate solution obtained by the proposed algorithm for (5) are shown in Figure 8(a). In addition, Figure 8(b) shows the accuracy of the solutions in terms of residual errors. It shows the accuracy of the solutions as the errors are approaching zero. Further, to validate the efficiency, absolute errors in solutions of BLM-NN are dictated through Table 3. It can be observed that the results of the design scheme overlap the exact solutions with minimum absolute errors as compared to the Haar wavelet collocation method [56] and the Bernoulli wavelet operational matrix method [57].Figure 8
(a) Approximate solutions and (b) absolute errors in our solutions for Example2.
(a)(b)Table 3
Comparison of absolute errors in solutions obtained by the BLM-NN algorithm with the Bernoulli wavelet operational matrix method and the Haar wavelet collocation method.
tN = 08N = 16N = 32N = 64BLM-NNHWCMBWOMHWCMBWOMHWCMBWOMHWCMBWOM1.00E-012.27E-041.10E-036.55E-052.00E-041.83E-056.96E-055.05E-061.53E-056.40E-082.00E-014.76E-041.70E-031.33E-045.00E-043.69E-051.17E-041.02E-053.69E-054.82E-083.00E-016.92E-042.50E-031.93E-048.00E-045.34E-051.75E-041.48E-055.42E-052.64E-084.00E-018.72E-044.00E-032.43E-049.00E-046.73E-052.74E-041.87E-055.91E-053.02E-105.00E-011.02E-035.30E-032.83E-041.40E-037.88E-053.52E-042.20E-059.10E-052.98E-086.00E-011.13E-035.90E-033.15E-041.20E-038.79E-053.87E-042.46E-058.28E-055.82E-087.00E-011.21E-035.30E-033.38E-041.70E-039.48E-053.58E-042.66E-051.14E-042.01E-088.00E-011.26E-035.80E-033.54E-041.90E-039.96E-053.96E-042.81E-051.26E-045.93E-089.00E-011.28E-038.00E-033.63E-041.60E-031.03E-045.36E-042.91E-051.12E-042.10E-06Example 3.
Consider the following nonlinear multiterm fractional differential equation as follows:(7)aD2ux+bxDα2ux+cxDux+exDα1ux+kxux=fx,x∈0,T,where 0<α1≤1,1<α2≤2 and fx are defined as follows:(8)fx=a−bxΓ3−α2x2−α2−cxx−exΓ3−α1x2−α1+kx2−12x2,with initial conditions as follows:(9)y0=2,y′0=0.for a=1,bx=x1/2,cx=x1/3,ex=x1/4,kx=x1/5,α1=0.5,andα2=1.5; the exact solution of (7) is yt=2−1/2x2. Figure 9 depicts the comparison of exact and numerical solutions obtained by the design algorithm for (7). The results calculated by the BLM-NN algorithm are compared with those obtained by the generalized block pulse operational matrix method [58] as shown in Table 4. The absolute errors lie around 10−7 to 10−8. The values of the performance function in terms of mean square error are shown in Table 2. The results in terms of computational complexity and absolute errors show the accuracy of the proposed algorithm in calculating solutions to fractional differential equations.Figure 9
(a) Approximate solutions and (b) absolute errors in our solutions for Example3.
(a)(b)Table 4
Comparison of absolute errors in the solutions of the BLM-NN algorithm with the generalized block pulse operational matrix method for different step sizes.
xh = 0.1h = 0.05h = 0.025h = 0.0125h = 0.00625BLM-NNGBPOM0.59.81E-082.40E-036.08E-041.52E-043.82E-059.55E-061.56.02E-072.00E-035.06E-041.27E-043.17E-057.95E-062.54.38E-071.60E-034.41E-041.10E-042.76E-056.90E-063.52.69E-071.60E-034.00E-041.03E-052.58E-056.46E-064.59.03E-071.60E-034.00E-059.99E-042.50E-056.25E-06Example 4.
Consider the following system of fractional differential equation:(10)DCυu1x=u2x,(11)DCυu2x=−u1x−u2x+xυ+1+πcscπυx1−υΓ−υ−1Γ2−υ+πxcscπυΓ−υ−1,with the following equation:(12)u10=0,u20=0.
The exact solutions of (10) and (11) are given as follows:(13)u1x=x1+υ,u2x=πυυ+1cscπυΓ1−υx.
We have solved this problem by considering different cases based on orders of derivative, i.e., Case Iυ=1/4, Case II υ=1/2, Case III υ=2/3, and Case IV υ=9/10.
Approximate solutions obtained by the BLM-NN algorithm foru1x and u2x are dictated in Table 5. The comparison or fitting of analytical solutions with approximate solutions is plotted in Figure 10. It can be seen that high-overlapping solutions with a minimum absolute error are obtained. Figure 11 represents the error histograms for different cases. The values of absolute errors for each case of (10) and (11) lie around 10−5 to 10−6, 10−5 to 10−7, 10−4 to 10−6, 10−5 to 10−6, 10−5 to 10−7, 10−6 to 10−8, 10−6 to 10−8, and 10−6 to 10−10, respectively. The smoothness of the algorithm has been detected from the convergence of the mean square error of the objective function. Figure 12 dictates that validated performance for each case of u1x and u2x are 1.3684e−12, 1.4825e−12, 3.123e−12, 4.9879e−11, 3.1169e−12, 1.46e−12, 1.1018e−12, and 2.661e−12, respectively. Further, details of performance indices are provided in Table 5. The values of the gradient for each case are 1.84e−07, 3.49e−06, 1.28e−05, 1.43e−07, 2.02e−06, 4.23e−07, 8.03e−07, and 3.33e−07. From Table 6, it can be seen that the values of mu for each case at 1000 epochs lie around 10−11 to 10−13. Regression analysis shown in Figure 13 further validates the efficiency and correctness of the technique.Table 5
Approximate solutions obtained by the proposed algorithm for different cases of the system of FDEs given in Example4.
xCase ICase IICase IIICase IVu1xu2xu1xu2xu1xu2xu1xu2x0.0000000000.50.420 4480.566 5020.353 5530.664 670.314 980.752 2880.267 9430.913 6781.011.133 00311.329 3411.504 57511.827 3551.51.660 0231.699 5051.837 1171.994 0111.965 5562.256 8632.160 5952.741 0332.02.378 4142.266 0062.828 4272.658 6813.174 8023.009 1513.732 1323.654 712.53.143 5842.832 5083.952 8473.323 3514.605 0393.761 4395.702 7724.568 3883.03.948 2223.399 0095.196 1523.988 0216.240 2514.513 7268.063 6265.482 0653.54.787 2383.965 5116.547 94.652 6918.068 2645.266 01410.807 66.395 7434.05.656 8544.532 01285.317 36210.079 376.018 30213.928 817.309 424.56.554 1395.098 5149.545 9425.982 03212.265 566.770 5917.422 238.223 0985.07.476 7445.665 01511.180 346.646 70214.620 097.522 87721.283 59.136 7755.58.422 7396.231 51712.898 647.311 37217.137 128.275 16525.508 7510.050 456.09.390 5076.798 01914.696 947.976 04219.811 569.027 45330.094 5210.964 13Figure 10
Comparison of analytical solution with approximate solution foru1a−d and u2e−h for multiple orders of Example 4. (a) Case I, (b) case II, (c) case III, (d) case IV, (e) case I, (f) case II, (g) case III, and (h) case IV.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 11
Error histogram between target values and approximated values for multiple orders of equations (10) and (11). (a) Case I, (b) case II, (c) case III, (d) case IV, (e) case I, (f) case II, (g) case III, and (h) case IV.
(a)(b)(c)(d)(e)(f)(g)(h)Figure 12
Convergence of performance in terms of mean square error foru1x (a–d) and u2x (e–h) for multiple orders of Example 4. (a) Case I, (b) case II, (c) case III, (d) case IV, (e) case I, (f) case II, (g) case III, and (h) case IV.
(a)(b)(c)(d)(e)(f)(g)(h)Table 6
Statistical analysis of performance measures including MSE, gradient, mu, number of iterations, and time taken by the system for obtaining the results of Example4.
Performance measuresCase ICase IICase IIICase IVu1xu2xu1xu2xu1xu2xu1xu2xHidden neurons6060606060606060Training8.67E-121.04E-111.64E-119.39E-133.88E-118.32E-133.83E-111.74E-12Validation1.37E-123.12E-121.48E-121.46E-123.12E-121.10E-124.99E-112.66E-12Testing8.52E-134.94E-131.06E-111.37E-128.76E-121.32E-124.29E-112.66E-12Gradient1.84E-073.49E-061.28E-051.43E-072.02E-064.23E-078.03E-073.33E-07Mu1.00E-131.00E-131.00E-131.00E-111.00E-121.00E-121.00E-111.00E-11Epochs100010005021000635100010001000Regression11111111Time (s)6 s5 s2 s5 s3 s5 s5 s5 sFigure 13
Regression analysis of cases I, II, III, and IV of Example4. (a) Case I, (b) case II, (c) case III, (d) case IV, (e) case I, (f) case II, (g) case III, and (h) case IV.
(a)(b)(c)(d)(e)(f)(g)(h)


## 4. Conclusion

In this paper, we have designed an integrated soft computing technique based on supervised learning. The computational strength of neural networks is utilized by the backpropagated Levenberg–Marquardt (BLM) algorithm to find approximate solutions for nonlinear multi-order fractional differential equations. The working procedure of BLM-NN algorithms is categorized into two steps in which the reference solution is generated by using analytical solutions. Furthermore, the data set is used by the BLM algorithm for validation, testing, and training of approximate solutions. Multiple figures, in terms of approximate solutions, curve fitting of analytical solutions and output data, error histograms, and regression and convergence of performance, are plotted to validate the efficiency of the design scheme. The tabulated data and figures dictate the accuracy, efficiency, and robustness of the design paradigm.In the future, the authors would like to extend the concept of soft computing based on neural networks to solve the mathematical models represented by partial differential equations and partial fractional differential equations.



---

*Source: 2710576-2022-01-19.xml*