Int. J. Anal. Appl. (2022), 20:51 

 

 

Received: Jul. 14, 2022. 

2010 Mathematics Subject Classification. 90C30. 

Key words and phrases. nonlinear programming problem; fractional gradient-based system; residual power series method; 

optimal solution; Caputo fractional derivative. 

 

https://doi.org/10.28924/2291-8639-20-2022-51 © 2022 the author(s) 

ISSN: 2291-8639  

1 

 

Analytical Solution of Nonlinear Fractional Gradient-Based System Using Fractional Power 

Series Method 

 

Radwan Abu-Gdairi* 

Department of Mathematics, Faculty of Science, Zarqa University, Zarqa 13132, Jordan 

*Corresponding author: rgdairi@zu.edu.jo 

ABSTRACT. This paper adapted a reliable treatment technique, called the fractional residual power series, to 

the fractional gradient-based system in solving a class of nonlinear programming model in Caputo’s sense. 

The gradient-based system has been constructed to transform the nonlinear programming problem with 

equality constraints to unconstrained optimization problem, and then the fractional residual power series 

method is implemented to obtain the essential behavior of underlying problem. The proposed methods have 

been applied effectively to produce optimal solution in rapidly convergent fractional series representations 

without linearization, or any limitations. To confirm the performance of the proposed methods, some 

optimization problems are tested. Further, numerical comparisons with other existing methods are also given. 

The results exhibit that the FRPS method is easy, simple, effective, and fully compatible with the complexity 

of such models. 

 

1. INTRODUCTION 

Over the past years, the topic of optimization has had the interest of many scholars in various 

applications of science and technology and associated with several categories of optimization 

problems. Besides, many effective methods have been developed to find the optimal solution to these 

problems. For more details see [1-6]. The gradient-based approach is one such method that is used 

to solve nonlinear programming (NLP) problems. Transforming the optimization problem into a 

system of ordinary differential equations (ODEs) is the basic idea of this method, which is equipped 

https://doi.org/10.28924/2291-8639-20-2022-51


2 Int. J. Anal. Appl. (2022), 20:51 
 

with ideal conditions to reach the optimal solutions to this problem [7-13]. Fractional calculus is a 

generalization of the derivatives and integrals of an arbitrary system. Recently, the subject of 

fractional calculus has received the attention of scientists and engineers because of its important 

applications in various fields, whether science or engineering [14-19]. Many real-life problems in 

various fields of applied science have been modeled using fractional differential equations (FDEs), 

which are generalizations of ODEs. 

For describing the behavior of the unknowns of FDEs, many researchers usually implement 

some numerical or numerical analytical techniques instead. In this regard, some recent techniques 

are proposed for solving FDEs. Among them decomposition technique, symmetric perturbation 

technique, variable frequency technique, and partial differential transformation technique [20-26]. 

Further, more applications and promising approaches are utilized to treat the nonlinear fractional 

gradient-based systems of FDEs could be founded in the references [27-30]. Motivated by the existing 

techniques, the main contribution of this article is to transform equality constrained NLP problem to 

unconstrained optimization problem by identifying a penalty function, then construct a gradient-

based system of FDEs. Besides the, the fractional residual power series (FRPS) technique is applied 

to provide us the accordance between the optimal solution of the NLP problem and the power series 

solution of the obtained FDEs system. 

FRPS technique is one of the modern numeric-analytic techniques was initially proposed in 

[31] to investigate the sequential solution of fuzzy differential equations of both first and second 

degree. It has been used to generate accurate approximate solutions in terms of fractional series 

formulas for several kinds for linear and non-linear FDEs, Partial FDEs and Fuzzy FDEs (see [32-

36]). This scheme is used basically the residual-error function and employed the fractional 

differentiation in each stage in determining the coefficient of the suggested series expansion without 

linearity, division, or perturbation (see [37-42]). It does not require any converting while switching 

from the lower order to the higher order; as a result, the method can be applied directly to the given 

problem by choosing an initial guess approximation. FRPS is quick and easy calculation to find series 

solutions via utilizing software package. Also, different Taylor series method, FRPS needs easy 

computation state with high reliability and less time. 

The organization of this paper was to present a brief presentation of some basic and necessary 

definitions and properties in fractional calculus in section 2, in addition to the fact that the central 

problem in this paper was presented in section 3. Section 4 presents the details of the application of 



3 Int. J. Anal. Appl. (2022), 20:51 
 

the proposed technique and provides a clear and simple algorithm for the basic steps of this technique. 

Clarifying the applicability and efficiency of the proposed technique by comparing the results derived 

from it for some numerical examples with the results of the fourth degree Rung-Kutta method, done 

in section 5. Section 6 is designed to provide some concluding observations. 

 

2. PRELIMINARIES 

In this section, we present the definition of the Riemann-Liouville fractional integral operator, 

the Caputo partial derivative, and some of their properties [43-49]. Throughout this paper, the symbol 

𝑅 denotes the set of real numbers, 𝑁 the set of integers, and Γ(. ) is a gamma function. For more 

details, please see [45-49] and references therein. 

Definition 1. Let 𝛼 ≥ 0 and for each 𝛼, 𝑥 ∈ ℝ. The Riemann-Liouville fractional integral operator of 

order, for a function 𝑢(𝑥), is defined as follow  

𝐼𝛼 𝑢(𝑥) =
1

Γ(𝛼)
∫ (𝑥 − 𝑠)1−𝛼 𝑢(𝑠)𝑑𝑠

𝑥

0
, 𝑡 > 0,  

In particular, if  𝛼 = 0, then 𝐼𝛼 𝑢(𝑡) = 𝑢(𝑡). 

Definition 2. The Caputo fractional derivative of a function 𝑢(𝑥) with  𝑚 − 1 < 𝛼 ≤ 𝑚, 𝑚 ∈ ℕ, is 

defined as  

𝐷𝛼 𝑢 =  {

1

𝛤(𝑚 − 𝛼)
∫ (𝑥 − 𝑠)𝑚−𝛼−1𝑢(𝑚)(𝑠)𝑑𝑠 , 𝑚 − 1 < 𝛼 < 𝑚,    𝑥 > 0

𝑥

0

.

𝑢(𝑚)(𝑥)                           ,           𝛼 = 𝑚   .

 

 

On the other hand, the operator 𝐷𝛼 has some basic properties such as, for any real number 𝐴, 

then we have 𝐷𝛼 𝐴 = 0, and for 𝑢(𝑥) = 𝑥𝑘 ,   𝑘 ≥ −1, we have 𝐷𝛼 𝑥𝑘 =
𝛤(𝑘+1)

𝛤(𝑘+1−𝛼)
𝑥 𝑘−𝛼. Moreover, the 

Caputo fractional derivative has the linearity property, this means, for each constant 𝛾 and 𝜇 we 

have 𝐷𝛼 [𝛾𝑢(𝑥) + 𝜇𝑣(𝑥)] = 𝛾𝐷𝛼 𝑢(𝑥) + 𝜇𝐷𝛼 𝑣(𝑥). 

Lemma 1. For 𝑚 − 1 < 𝛼 ≤ 𝑚 and 𝑢(𝑥) ∈ 𝐶𝑚 [0, ∞), 𝑚 ∈ ℕ, then we have  

1. 𝐷𝛼 𝐼𝛼 𝑢(𝑥) = 𝑢(𝑥), 

2. 𝐼𝛼 𝐷𝛼 𝑢(𝑥) = 𝑢(𝑥) − ∑ 𝑢(𝑘)𝑚−1𝑖=0 (0
+)

𝑥𝑘

𝑘!
, 𝑥 > 0. 

 

3. OPTIMIZATION PROBLEM 

The second part presents the details of the optimization problem to be studied in this 

paper. We consider the following NLP problem with equality constrains 



4 Int. J. Anal. Appl. (2022), 20:51 
 

min 𝑢(𝑥)  𝑠. 𝑡.  𝑤(𝑥) = 0, 𝑥 ∈ ℝ𝑛 , (2) 

where 𝑥 ∈ ℝ𝑛 is decision variable, 𝑢(𝑥) is a vector-valued function of a real variable, and 𝑤 =

(𝑤1, 𝑤2, … , 𝑤𝑝)
𝑇

: ℝ𝑛 → ℝ𝑝 (𝑝 ≤ 𝑛) are twice continuously differentiable function such that whose 

gradient ∇𝑤(𝑥) has full rank. We assume that a feasible region of (2) is nonempty bounded set. 

The penalty function can be defined as  

𝑃(𝑥) =  𝑢(𝑥) +  𝜙(𝑥), (3) 

where 𝜙(𝑥) is the penalty term defined on ℝ𝑛 and has the following property 

𝜙(𝑥) = {
0,                                𝑤(𝑥) = 0
𝑛𝑜𝑛𝑛𝑒𝑔𝑎𝑡𝑖𝑣𝑒, 𝑤(𝑥) ≠ 0.

 
(4) 

One can be defined the penalty term 𝜙(𝑥) as 

𝜙(𝑥) =  
1

𝜇
‖𝑤(𝑥)‖𝜇

𝜇
=  

1

𝜇
∑ |𝑤𝑖 (𝑥)|

𝜇
𝑝

𝑖=1
, 

(5) 

where 𝜇 > 0 is a constant.  

A well-known penalty function for the problem (2) has been defined as 

𝑃(𝑥, 𝜂𝑚) = 𝑢(𝑥) +
𝜂𝑚
𝜇

 ∑ |𝑤𝑖 (𝑥)|
𝜇

𝑝

𝑖=1
, 

(6) 

where the penalty parameter 𝜂𝑚 satisfying the inequality 0 < 𝜂𝑚 < 𝜂𝑚+1 for all 𝑚,        𝜂𝑚 →

∞. It is worth mentioning that the penalty parameter 𝜂𝑚 can be chosen randomly. Thus, 𝜂𝑚can be 

chosen in positive constant depends on the difficulty of minimizing the penalty function at every 

iteration. It can be shown that under some suitable conditions the solution of the equality constrains 

NLP problem (2) are solution of the following unconstrained optimization problem  

min 𝑃(𝑥, 𝜂𝑚) = 𝑢(𝑥) +
𝜂𝑚
𝜇

 ∑ |𝑤𝑖 (𝑥)|
𝜇

𝑝

𝑖=1
     𝑠. 𝑡   𝑥 ∈ ℝ𝑛 . 

(7) 

We assume that the unconstrained optimization problem (7) for each 𝑚, has a solution and 

we denote it by 𝑥𝑚. The main results that connect the solutions of the equality constrained NLP 

problem (2) and unconstrained problem (7) present in the following theorem.  

Theorem 1. Let {xm} be a sequence generated by the penalty method of the unconstrained problem 

(7). Then any limit point of the sequence is a solution to the equality constrained NLP problem (2).  

The author in [12] showed that the unconstrained optimization problem (7) can be described by the 

following gradient based dynamic system of ODEs 

𝑑𝑥(𝑡)

𝑑𝑡
= −∇𝑥 (𝑢(𝑥) +

𝜂𝑚
𝜇

 ∑ |𝑤𝑖 (𝑥)|
𝜇

𝑝

𝑖=1
) , 𝑥(𝑡0) = 𝑏 ∈ ℝ

𝑛 
(8) 



5 Int. J. Anal. Appl. (2022), 20:51 
 

where  ∇𝑥 is the gradient vector with respect to the 𝑥 ∈ ℝ
𝑛. We can describe the system (8) by an 

approach based on fractional gradient based dynamic system by the following system of FDEs 

𝐷𝑡
𝛼 𝑥(𝑡) = −∇𝑥 (𝑢(𝑥) +

𝜂𝑚

𝜇
 ∑ |𝑤𝑖 (𝑥)|

𝜇𝑝
𝑖=1 ) , 𝑥(𝑡0) = 𝑏 ∈ ℝ

𝑛, 0 < 𝛼 ≤ 1.    (9) 

The system (9) has an equilibrium point 𝑥𝑒, if 𝑥𝑒 ∈ ℝ
𝑛 is satisfies the right-hand side of the system. 

A more convenient form of the system (9) can be expressed as follows 

𝐷𝑡
𝛼 𝑥𝑖 (𝑡) = 𝑓𝑖 (𝑡, 𝜂𝑚, 𝑥1, 𝑥2, … , 𝑥𝑛), 𝑖 = 1,2, … , 𝑛, 𝑥(𝑡0) = 𝑏 ∈ ℝ

𝑛 , 0 < 𝛼 ≤ 1.   (10) 

The stable equilibrium point of the fractional system (10) is acquired with the RPS technique. 

 

4. FRPS TECHNIQUE  

This subsection devoted to applying the RPS method to derive analytic solution of the system 

of FDEs (10). We begin by propose the definition of fractional power series. 

Definition 3. A power series (PS) expansion at 𝑡 = 𝑡0 of the following form 

∑ 𝑐𝑖 (𝑡 − 𝑡0)
𝑗𝛼

∞

𝑖=0

= 𝑐0 + 𝑐1(𝑡 − 𝑡0)
𝛼 + 𝑐2(𝑡 − 𝑡0)

2𝛼 + ⋯, 
 

(11) 

for 0 ≤ 𝑛 − 1 < 𝛼 ≤ 𝑛, 𝑛 ∈ ℕ, 𝑡 ≥ 𝑡0 where 𝑐𝑗’s are constants called the fractional power series 

(FPS). 

Theorem 2. Suppose that 𝑥(𝑡) has a FPS representation at 𝑡 = 𝑡0 of the form 

𝑥(𝑡) = ∑ 𝑐𝑗 (𝑡 − 𝑡0)
𝑗𝛼

∞

𝑖=0

, 𝑡0 ≤ 𝑡 ≤ 𝑡0 + 𝑅. 
 

(12) 

If 𝑥(𝑡) ∈ 𝐶[𝑡0, 𝑡0 + 𝑅), and 𝐷𝑡
𝑖𝛼 𝑥(𝑡) ∈ 𝐶(𝑡0, 𝑡0 + 𝑅), for 𝑖 = 0,1,2, …, then  the coefficients 𝑐𝑖 in Eq. 

(6) will take the form 𝑐𝑖 =
𝐷𝑖𝛼𝑥(𝑡0)

𝛤(𝑖𝛼+1)
,  where 𝐷𝑖𝛼 = 𝐷𝛼 ∙ 𝐷𝛼 ∙∙∙ 𝐷𝛼 (𝑖-times) and 𝑅 is the convergent 

radius. 

Theorem 3. If 𝐾 ∈ (0,1), such that ‖𝑥𝑖+1(𝑡)‖ ≤ 𝐾‖𝑥𝑖 (𝑡)‖ ∀𝑖 ∈ ℕ and 0 < 𝑡 < Τ < 1, then the series 

of numerical solutions converges to an exact solution  

Proof. We notice that ∀ 0 < 𝑡 < Τ < 1, 

‖𝑥(𝑡) − 𝑥𝑖 (𝑡)‖ = ‖ ∑ 𝑥𝑖 (𝑡)

∞

𝑚=𝑖+1

‖ ≤ ∑ ‖𝑥𝑖 (𝑡)‖

∞

𝑚=𝑖+1

≤ ‖𝑏‖ ‖ ∑ 𝐾𝑚
∞

𝑚=𝑖+1

‖  

                                   =
𝐾𝑘+1

1−𝐾
‖𝑏‖ → 0 𝑎𝑠 𝑘 → ∞.  

 



6 Int. J. Anal. Appl. (2022), 20:51 
 

Now, to utilize the FRPS technique to solve the system (10), we substitute the FPS 

expansion (12) among the recurrent fractional differentiation of truncation residual functions. 

Suppose that the solution of the system (10) about the initial point 𝑡 = 𝑡0 takes the form 

𝑥𝑖 (𝑡) = ∑ 𝑐𝑖𝑗

∞

𝑗=0

(𝑡 − 𝑡0)
𝑗𝛼 , 𝑖 = 1,2, … , 𝑛. 

 

(13) 

Differentiate the expansion form (13) for  𝑡 > 0 within the interval of convergence, we get 

𝐷𝑡
𝛼 𝑥𝑖 (𝑡) = ∑

𝑐𝑖𝑗 Γ(𝑗𝛼 + 1)

Γ((𝑗 − 1)𝛼 + 1)

∞

𝑗=1

(𝑡 − 𝑡0)
(𝑗−1)𝛼 , 𝑖 = 1,2, … , 𝑛. 

 

(14) 

Hence,  

𝐷𝑡
2𝛼 𝑥𝑖 (𝑡) = ∑

𝑐𝑖𝑗 Γ(𝑗𝛼 + 1)

Γ((𝑗 − 2)𝛼 + 1)

∞

𝑗=2

(𝑡 − 𝑡0)
(𝑗−2)𝛼 , 𝑖 = 1,2, … , 𝑛, 

 

(15) 

where 𝐷𝑡
2𝛼 = 𝐷𝑡

𝛼 𝐷𝑡
𝛼. Therefore, we suppose that the approximate solution of 𝑥𝑖 (𝑡) can be 

constructed by the following series  

𝑥𝑖𝑘 (𝑡) = ∑ 𝑐𝑖𝑗

𝑘

𝑗=0

(𝑡 − 𝑡0)
𝑗𝛼 , 𝑖 = 1,2, … , 𝑛. 

 

(16) 

Since 𝑥𝑖 (𝑡) satisfy the initial condition 𝑥𝑖 (𝑡0) = 𝑏𝑖 , 𝑖 = 1,2, … , 𝑛, then 𝑐𝑖0 = 𝑏𝑖 and series (16) can be 

rewritten by 

𝑥𝑖𝑘 (𝑡) = 𝑏𝑖 + ∑ 𝑐𝑖𝑗

𝑘

𝑗=1

(𝑡 − 𝑡0)
𝑗𝛼 , 𝑘 = 1,2,3, . . . , 𝑖 = 1,2, … , 𝑛. 

 

(17) 

According the RPS technique, we define the residual function, 𝑅𝑒𝑠𝑥𝑖 (𝑡), about 𝑡 = 𝑡0 for the system 

(10) as follows 

𝑅𝑒𝑠𝑥𝑖 (𝑡) = 𝐷𝑡
𝛼 𝑥𝑖 (𝑡) + ∇𝑥𝑖 𝑃(𝑥𝑖 , 𝜂), 𝑖 = 1,2, … , 𝑛. (18) 

by assuming the penalty parameter 𝜂𝑚 = 𝜂 (constant). The 𝑘-th residual function, 𝑅𝑒𝑠𝑥𝑖,𝑘 (𝑡), can 

be defined by   

𝑅𝑒𝑠𝑥𝑖,𝑘(𝑡) = 𝐷𝑡
𝛼 𝑥𝑖𝑘 (𝑡) + ∇𝑥𝑖𝑘 𝑃(𝑥𝑖𝑘 , 𝜂), 𝑘 = 1,2,3, … , 𝑖 = 1,2, … , 𝑛. (19) 

Note that, 𝑅𝑒𝑠𝑥𝑖 (𝑡) = 0 and 𝐿𝑖𝑚𝑘→∞𝑅𝑒𝑠𝑥𝑖,𝑘 (𝑡) = 𝑅𝑒𝑠𝑥𝑖 (𝑡), 𝑖 = 1,2, … , 𝑛, for all 𝑡 ≥ 0. As a matter 

of fact, it yields the following algebraic system 

𝐷𝑡
(𝑘−1)𝛼

𝑅𝑒𝑠𝑥𝑖,𝑘 (𝑡0) = 0, 𝑘 = 1,2,3, … , 𝑖 = 1,2, … , 𝑛. 
(20) 



7 Int. J. Anal. Appl. (2022), 20:51 
 

Through the algebraic system (20) we can obtained the coefficient 𝑐𝑗 . 𝑗 = 1,2,3, … , 𝑘 by applying the 

following steps. 

Algorithm 1. Algorithm of finding the coefficients 𝒄𝒊𝒋 of the k-th truncated series (16). 

Step 1: Substitute the 𝑘-th residual expansion (19) into (18), where 

 

(21) 

 

𝑥𝑖𝑘 (𝑡) = 𝑐𝑖0 + ∑ 𝑐𝑖𝑗

𝑘

𝑗=1

(𝑡 − 𝑡0)
𝑗𝛼 , 𝑐𝑖𝑗 =

𝐷𝑗𝛼 𝑏𝑖
𝛤(𝑗𝛼 + 1)

. 

Step 2: Find the formula  𝐷𝑡
(𝑘−1)𝛼

(𝑅𝑒𝑠𝑥𝑖,𝑘 (𝑡0)) , 𝑖 = 1,2, … , 𝑛. 

Step 3: The coefficient 𝑐𝑖𝑗 can be obtained by solve the algebraic system 

(22) 𝐷𝑡
(𝑘−1)𝛼

𝑅𝑒𝑠𝑥𝑖,𝑘 (𝑡0) = 0, 𝑘 = 1,2,3, … , 𝑖 = 1,2, … , 𝑛. 

Step 4: Substitute the obtained values of 𝑐𝑖𝑗 back into Eq. (17). 
 

 

5. NUMERICAL IMPLEMENTATION AND RESULTS 

This section is designed to apply the proposed technique, FRPS, and evaluate its performance 

and accuracy in deriving some numerical solutions for a number of examples and comparing the 

results obtained with the analytical solutions of the examples. We used the Mathematica software 

package to perform numerical and symbolic calculations. 

Example 1. Consider the following NLP problem   

min 𝑢(𝑥) = 100(𝑥1
2 − 𝑥2)

2 + (𝑥1 − 1),   
 

subject to 
 

𝑤(𝑥) = 𝑥1(𝑥1 − 4) − 2𝑥2 + 12 = 0. 

 

 

 

(23) 

According to (6), the correspond penalty function for the problem (23) at 𝜇 = 2 is given by 

𝑃(𝑥1, 𝑥2, 𝜂) = 100(𝑥1
2 − 𝑥2)

2 + (𝑥1 − 1) +
1

2
𝜂(𝑥1

2 − 4𝑥1 − 2𝑥2 + 12)
2, 

(24) 

where the penalty variable 𝜂 ∈ ℝ+, 𝜂 → ∞. Hence, we get the following correspond system of FDEs 

𝐷𝑡
𝛼 𝑥1(𝑡) = −400(𝑥1

2 − 𝑥2)𝑥1 − 2(𝑥1 − 1) 

                                                      −𝜂(2𝑥1 − 4)(𝑥1
2 − 4𝑥1 − 2𝑥2 + 12),  

 

𝐷𝑡
𝛼 𝑥2(𝑡) = 200(𝑥1

2 − 𝑥2) + 2𝜂(𝑥1
2 − 4𝑥1 − 2𝑥2 + 12),  

𝑥1(0) = 𝑥2(0) = 0,  (25) 



8 Int. J. Anal. Appl. (2022), 20:51 
 

 where 0 < 𝛼 ≤ 1.  

  To apply the FRPS technique to solve the system (25), we suppose that  

𝑥1(𝑡) = ∑ 𝑐1𝑗

∞

𝑗=0

(𝑡 − 𝑡0)
𝑗𝛼 , 𝑥2(𝑡) = ∑ 𝑐2𝑗

∞

𝑗=0

(𝑡 − 𝑡0)
𝑗𝛼 . 

 

(26) 

   Thus, the k-th truncated series is 

𝑥1𝑘 (𝑡) = ∑ 𝑐1𝑗

𝑘

𝑗=1

(𝑡 − 𝑡0)
𝑗𝛼 , 𝑥2𝑘 (𝑡) = ∑ 𝑐2𝑗

𝑘

𝑗=1

(𝑡 − 𝑡0)
𝑗𝛼 . 

 

(27) 

   Consequently, the 𝑘-th residual function, 𝑅𝑒𝑠𝑥𝑖,𝑘 (𝑡), can be defined by  

𝑅𝑒𝑠𝑥1,𝑘 (𝑡) = 𝐷𝑡
𝛼 𝑥1𝑘 (𝑡) + 400(𝑥1

2 − 𝑥2)𝑥1 + 2(𝑥1 − 1) 

                                         +𝜂(2𝑥1 − 4) (𝑥1
2 − 4𝑥1 − 2𝑥2 + 12),   

 

𝑅𝑒𝑠𝑥2,𝑘 (𝑡) = 𝐷𝑡
𝛼 𝑥2𝑘 (𝑡) − 200(𝑥1

2 − 𝑥2) − 2𝜂(𝑥1
2 − 4𝑥1 − 2𝑥2 + 12), (28) 

    Finally, Algorithm 1 applied to find the coefficients 𝑐𝑖𝑗, then getting the approximate solution of 

the system (25).  

The efficiency of FRPS approach is introduced via establishing some numerical comparisons 

for the obtained results and the results obtained by Runge-Kutta approach and listed in Table 1, for 

problem (23), from this table we noted that the FRPS solutions are compatible with the exact 

solution more that the RK4 solutions. The geometric behavior of the FRPS approximate solutions 

are shown against the exact and RK4 solutions as in Figure 1. Clearly, the figure indicates that the 

exact and FRPS approximate solutions are in good agreement for different values of 𝛼, comparing 

with the exact and RK4 approximate solutions over the domain of 𝛼.  

Table 1. Comparison of 𝑥1(𝑡) and 𝑥2(𝑡) for problem (23) between FRPS with RK4 solutions at 

𝜂 = 200 and 𝛼 = 1. 

 

𝑡 𝐹𝑅𝑃𝑆 𝑥1(𝑡) 𝐹𝑅𝑃𝑆 𝑥2(𝑡) 𝑅𝐾4 𝑥1(𝑡) 𝑅𝐾4 𝑥2(𝑡) Absolute error  

𝑥1(𝑡) 

Absolute 

error  𝑥2(𝑡) 

0.0000  0.000000 0.000000 0.000000  0.000000  0.000000 0.000000 

0.0005  1.971356 3.872061 1.970899  3.871887  0.000457 0.000174 

0.0010  1.978467 3.908321 1.978274  3.907993  0.000193 0.000328 

0.0013  1.98151 3.923481 1.981384  3.922554  0.000126 0.000927 

0.0015  1.984216 3.932718 1.983132  3.930578  0.001084 0.00214 

0.0020  1.9864 3.937884 1.984252  3.935654  0.002148 0.00223 



9 Int. J. Anal. Appl. (2022), 20:51 
 

 

  

(a)                                                                                   (b) 

Figure 1. Comparison of 𝑥1(𝑡) and 𝑥2(𝑡) for problem (23) between FRPS with RK4 solutions at 

𝜂 = 200 and 𝛼 = 1: (a) 𝑥1(𝑡), (b) 𝑥2(𝑡), solidline: FRPS, dashline: RK4. 

Example 2. Consider the following NLP problem: 

min 𝑢(𝑥) = (𝑥1 − 1)
2 + (𝑥1 − 𝑥2)

2 + (𝑥2 − 𝑥3)
2 + (𝑥3 − 𝑥4)

4 + (𝑥4 − 𝑥5)
4,  

subject to:  

𝑤1(𝑥) = 𝑥1 + 𝑥2
2 + 𝑥3

3 − 2 − 3√2 = 0  

𝑤2(𝑥) = 𝑥2 − 𝑥3
2 + 𝑥4 + 2 − 2√2 = 0,  

𝑤3(𝑥) = 𝑥1𝑥5 − 2 = 0   (29) 

   We define the penalty function for the problem (29) as  

𝑃(𝑥1, 𝑥2, 𝑥3, 𝑥4, 𝑥5, 𝜂) = (𝑥1 − 1)
2 + (𝑥1 − 𝑥2)

2 + (𝑥2 − 𝑥3)
2 + (𝑥3 − 𝑥4)

4 +

(𝑥4 − 𝑥5)
4 +

𝜂

2
[(𝑥1 + 𝑥2

2 + 𝑥3
3 − 2 − 3√2)

2
+  (𝑥2 − 𝑥3

2 + 𝑥4 + 2 − 2√2)
2

+

 (𝑥1𝑥5 − 2)
2],  

 

 

 

(30) 

where the penalty variable 𝜂 ∈ ℝ+, 𝜂 → ∞. Depending on this penalty function, the correspond 

FDEs system can be written as  

𝐷𝑡
𝛼 𝑥𝑖 (𝑡) = ∇𝑥 𝑖 𝑢(𝑥) + 𝜂 ∑ ∇𝑥 𝑖 𝑤𝑖 (𝑥)

5

𝑖=1
𝑤𝑖 (𝑥), 𝑥𝑖 (0) = 2, 𝑖 = 1,2,3,4,5,  

 

(31) 

where 0 < 𝛼 ≤ 1. We adapt the FRPS technique to the FDEs system (31) with penalty variable 𝜂 =

300 at fractional order derivative 𝛼 = 0.9, we acquired solutions as shown in Table 2. Figure 2 

present the obtained FRPS solutions 𝑥1(𝑡) and 𝑥4(𝑡) at various fractional derivative order 𝛼. 

Obviously, from Figure 2, the curves-FRPS approximate solutions consistent with each other and 

approach the exact curve with increasing fractional values to the integer-order value 𝛼 =  1.  



10 Int. J. Anal. Appl. (2022), 20:51 
 

Table 2. The obtained solution of NLP problem (29) by FRPS technique at 𝜂 = 300 and 

𝛼 = 0.9. 

𝑡 𝑥1(𝑡) 𝑥2(𝑡) 𝑥3(𝑡) 𝑥4(𝑡) 𝑥5(𝑡) 

0.000 2.00000 2.00000 2.00000 2.00000 2.00000 

1.000 1.19855 1.36387 1.47429 1.64587 1.67875 

2.000 1.19951 1.36981 1.47237 1.63658 1.67801 

3.000 1.20181 1.36473 1.48531 1.63782 1.67947 

4.000 1.21414 1.36588 1.48431 1.63842 1.68371 

5.000 1.26941 1.36947 1.48997 1.64889 1.68753 

6.000 1.27542 1.37842 1.49568 1.64998 1.68947 

 

  

(a)                                                                                     (b) 

 

Figure 2. The FRPS solutions: (a) 𝑥1(𝑡), (b) 𝑥4(𝑡); 𝛼 = 0.9 solid, 𝛼 = 0.8 dashed, and 𝛼 = 0.7 

dot-dashed, at 𝜂 = 300. 

Example 3. Consider the following NLP problem: 

𝑚𝑖𝑛. 𝑢(𝑡) = (𝑥1 − 20)
2 + (𝑥2 + 20)

2,  

subject to:  

 𝑤(𝑡) =
𝑥1

2

100
+

𝑥2
2

4
− 1 = 0, 

 

(32) 

   The penalty function of this problem can be written as  

𝑃(𝑥1, 𝑥2, 𝜂) = (𝑥1 − 20)
2 + (𝑥2 + 20)

2 +
𝜂

2
(

𝑥1
2

100
+

𝑥2
2

4
− 1)

2

, 
 

(33) 

where the penalty variable 𝜂 ∈ ℝ+, 𝜂 → ∞. Therefore, we can write the correspond FDEs 

system as  



11 Int. J. Anal. Appl. (2022), 20:51 
 

𝐷𝑡
𝛼 𝑥1(𝑡) = 2𝑥1 − 40 + 𝜂 (

𝑥1
3

5000
+

𝑥1𝑥2
2

200
−

𝑥1
50

), 
 

𝐷𝑡
𝛼 𝑥2(𝑡) = 2𝑥2 + 40 + 𝜂 (

𝑥1
2𝑥2

200
+

𝑥2
3

8
−

𝑥2
2

), 
 

𝑥1(0) = 𝑥2(0) = 0, (34) 

where 0 < 𝛼 ≤ 1.  

The FRPS utilize to construct the solution of this system and we get the following numerical 

solutions of the NLP problem (32) as shown in Table 3. It is very clear that the results obtained for 

example 3 indicated that the proposed method is simple, and its performance is very effective 

comparing with Runge-Kutta method.  

Table 3. Comparison of 𝑥1(𝑡) and 𝑥2(𝑡)  for problem (32) between FRPS with RK4 solutions 

𝜂 = 106 and 𝛼 = 1. 

 

𝑡 𝐹𝑅𝑃𝑆 𝑥1(𝑡) 𝐹𝑅𝑃𝑆 𝑥2(𝑡) 𝑅𝐾 𝑥1(𝑡) 𝑅𝐾 𝑥2(𝑡) Error 

𝑥1(𝑡) 

Error 𝑥2(𝑡) 

0.0  0.000000 0.000000 0.000000   0.000000  0.000000 0.000000 

1.0  5.850119 −0.47496  5.847465   −0.477547  0.002654 0.002584 

2.0  6.744776 −0.53968  6.741904   −0.542553  0.002872 0.002873 

4.0  8.460532 −0.58119  8.457551   −0.584133  0.002981 0.002948 

6.0  8.735192 −0.61163  8.732010   −0.615305  0.003182 0.003675 

7.0  8.980451 −0.63668  8.977187   −0.640466  0.003264 0.003782 

8.0 9.092537 −0.663578 9.090838 −0.661673 0.001699 0.001905 

10.0 9.304238 −0.691245 9.303232  −0.680066 0.001006 0.018393 

    

 

5.  CONCLUSIONS 

In this a novel paper, the corresponding system of FDEs for NLP was designed and analyzed 

under the meaning of Caputo derivative and then solved the target system via a modern efficient 

technique, named FRPS technique.  The benefit of utilizing the present technique is that it gives 

accurate convergence approximate solution to an exact solution with needs a small size of 

computation to the optimal solutions to the NLP problems. Three attractive NLP are considered to 

validate the applicability and superiority of the FRPS scheme. Simulation data and graphical 



12 Int. J. Anal. Appl. (2022), 20:51 
 

representations are discussed and indicated that the obtained results by FRPS approach are in good 

agreement with each other and with exact solution at various values of fractional derivative order 𝛼, 

as well as from the figures the results emphasized that obtained FRPS solutions more accurate from 

the obtained Runge-Kutta solutions versus the exact solution. Consequently, the presented technique 

has the ability to handle both linear and nonlinear fractional biological phenomena. In future works, 

the FRPS can be extended to generate accurate approximate solutions for systems of fractional 

partial differential equations.  

 

Conflicts of Interest: The author declares that there are no conflicts of interest regarding the 

publication of this paper. 

 

References 

[1] W. Sun, Y.X. Yuan, Optimization Theory and Methods: Nonlinear Programming, Springer-Verlag, New 

York (2006). 

[2] N. Özdemir, F. Evirgen, Solving NLP Problems With Dynamic System Approach Based on Smoothed 

Penalty Function, Selcuk J. Appl. Math. 10 (2009), 63–73. https://hdl.handle.net/20.500.12462/4852. 

[3] H. Yamashita, A Differential Equation Approach to Nonlinear Programming, Math. Program. 18 (1980), 

155–168. https://doi.org/10.1007/bf01588311. 

[4] M. Al-Smadi, O.A. Arqub, D. Zeidan, Fuzzy Fractional Differential Equations Under the Mittag-Leffler 

Kernel Differential Operator of the ABC Approach: Theorems and Applications, Chaos Solitons Fractals. 

146 (2021), 110891. https://doi.org/10.1016/j.chaos.2021.110891. 

[5] A.V. Fiacco, G.P. Mccormick. Nonlinear Programming: Sequential Unconstrained Minimization 

Techniques, John Wiley, New York (1968). 

[6] C.A. Botsaris, Differential Gradient Methods, J. Math. Anal. Appl. 63 (1978), 177–198. 

https://doi.org/10.1016/0022-247x(78)90114-2. 

[7] M. Al‐Smadi, Fractional Residual Series for Conformable Time‐Fractional Sawada–Kotera–Ito, Lax, and 

Kaup–Kupershmidt Equations of Seventh Order, Math Methods Appl. Sci. (2021). 

https://doi.org/10.1002/mma.7507. 

[8] A.A. Brown, M.C. Bartholomew-Biggs, ODE Versus SQP Methods for Constrained Optimization, J. 

Optim. Theory Appl. 62 (1989), 371–386. https://doi.org/10.1007/bf00939812. 

[9] J. Schropp, I. Singer, A Dynamical Systems Approach to Constrained Minimization, Numer. Funct. Anal. 

Optim. 21 (2000), 537–551. https://doi.org/10.1080/01630560008816971. 

https://hdl.handle.net/20.500.12462/4852
https://doi.org/10.1007/bf01588311
https://doi.org/10.1016/j.chaos.2021.110891
https://doi.org/10.1016/0022-247x(78)90114-2
https://doi.org/10.1002/mma.7507
https://doi.org/10.1007/bf00939812
https://doi.org/10.1080/01630560008816971


13 Int. J. Anal. Appl. (2022), 20:51 
 

[10] S. Wang, X.Q. Yang, K.L. Teo, A Unified Gradient Flow Approach to Constrained Nonlinear Optimization 

Problems, Comput. Optim. Appl. 25 (2003), 251–268. https://doi.org/10.1023/a:1022973608903. 

[11] M. Al-Smadi, O.A. Arqub, Computational Algorithm for Solving Fredholm Time-Fractional Partial 

Integrodifferential Equations of Dirichlet Functions Type With Error Estimates, Appl. Math. Comput. 342 

(2019), 280–294. https://doi.org/10.1016/j.amc.2018.09.020. 

[12] L. Jin, L.W. Zhang, X.T. Xiao, Two Differential Equation Systems for Equality-Constrained Optimization, 

Appl. Math. Comput. 190 (2007), 1030–1039. https://doi.org/10.1016/j.amc.2006.11.041. 

[13] M. Al-Smadi, A. Freihat, H. Khalil, S. Momani, R. Ali Khan, Numerical Multistep Approach for Solving 

Fractional Partial Differential Equations, Int. J. Comput. Methods. 14 (2017), 1750029. 

https://doi.org/10.1142/s0219876217500293. 

[14] N. Andrei, Gradient Flow Algorithm for Unconstrained Optimization, ICI Technical Report, 2004. 

[15] L. Jin, A Stable Differential Equation Approach for Inequality Constrained Optimization Problems, Appl. 

Math. Comput. 206 (2008), 186–192. https://doi.org/10.1016/j.amc.2008.09.007. 

[16] N. Özdemir, F. Evirgen, A Dynamic System Approach to Quadratic Programming Problems With Penalty 

Method, Bull. Malays. Math. Sci. Soc. (2) 33 (2010), 79–91. 

[17] S. Hasan, M. Al-Smadi, A. El-Ajou, et al. Numerical Approach in the Hilbert Space to Solve a Fuzzy 

Atangana-Baleanu Fractional Hybrid System, Chaos Solitons Fractals. 143 (2021) 110506. 

https://doi.org/10.1016/j.chaos.2020.110506. 

[18] A. Qazza, R. Hatamleh, M. Al-hawari, Dirichlet Problem in the Simply Connected Domain, Bounded by 

Unicursal Curve, Int. J. Appl. Math. 22 (2009), 599-614. 

[19] A. Qazza, R. Hatamleh, Dirichlet Problem in the Simply Connected Domain, Bounded by the Nontrivial 

Kind, Adv. Differ. Equ. Control Processes, 17 (2016), 177-188. 

[20] R. Saadeh, M. Alaroud, M. Al-Smadi, et al. Application of Fractional Residual Power Series Algorithm to 

Solve Newell–Whitehead–Segel Equation of Fractional Order, Symmetry. 11 (2019), 1431. 

 https://doi.org/10.3390/sym11121431. 

[21] R. Edwan, R. Saadeh, S. Hadid, M. Al-Smadi, S. Momani, Solving Time-Space-Fractional Cauchy Problem 

with Constant Coefficients by Finite-Difference Method, in: D. Zeidan, S. Padhi, A. Burqan, P. Ueberholz 

(Eds.), Computational Mathematics and Applications, Springer Singapore, Singapore, 2020: pp. 25–46. 

https://doi.org/10.1007/978-981-15-8498-5_2. 

[22] A. Burqan, A. El-Ajou, R. Saadeh, M. Al-Smadi, A New Efficient Technique Using Laplace Transforms 

and Smooth Expansions to Construct a Series Solution to the Time-Fractional Navier-Stokes Equations, 

Alexandria Eng. J. 61 (2022), 1069–1077. https://doi.org/10.1016/j.aej.2021.07.020. 

[23] R. Saadeh, A. Qazza, A. Burqan, A New Integral Transform: ARA Transform and Its Properties and 

Applications, Symmetry. 12 (2020), 925. https://doi.org/10.3390/sym12060925. 

https://doi.org/10.1023/a:1022973608903
https://doi.org/10.1016/j.amc.2018.09.020
https://doi.org/10.1016/j.amc.2006.11.041
https://doi.org/10.1142/s0219876217500293
https://doi.org/10.1016/j.amc.2008.09.007
https://doi.org/10.1016/j.chaos.2020.110506
https://doi.org/10.3390/sym11121431
https://doi.org/10.1007/978-981-15-8498-5_2
https://doi.org/10.1016/j.aej.2021.07.020
https://doi.org/10.3390/sym12060925


14 Int. J. Anal. Appl. (2022), 20:51 
 

[24] M. Al-Smadi, N. Djeddi, S. Momani, et al. An attractive numerical algorithm for solving nonlinear Caputo-

Fabrizio fractional Abel differential equation in a Hilbert space, Adv. Differ. Equ. 2021 (2021), 271. 

https://doi.org/10.1186/s13662-021-03428-3. 

[25] S.A. El-Wakil, A. Elhanbaly, M.A. Abdou, Adomian Decomposition Method for Solving Fractional 

Nonlinear Differential Equations, Appl. Math. Comput. 182 (2006), 313-324. 

https://doi.org/10.1016/j.amc.2006.02.055. 

[26] M. Al-Smadi, H. Dutta, S. Hasan, S. Momani, On Numerical Approximation of Atangana-Baleanu-Caputo 

Fractional Integro-Differential Equations Under Uncertainty in Hilbert Space, Math. Model. Nat. Phenom. 

16 (2021), 41. https://doi.org/10.1051/mmnp/2021030. 

[27] A. Khan, A. Khan, M. Sinan, Ion Temperature Gradient Modes Driven Soliton and Shock by Reduction 

Perturbation Method for Electron-Ion Magneto-Plasma, Math. Model. Numer. Simul. Appl. 2 (2022), 1-

12. https://doi.org/10.53391/mmnsa.2022.01.001. 

[28] F. Evirgen, Conformable Fractional Gradient Based Dynamic System for Constrained Optimization 

Problem, Acta Phys. Pol. A. 132 (2017), 1066–1069. https://doi.org/10.12693/aphyspola.132.1066. 

[29] F. Evirgen, M. Yavuz, An Alternative Approach for Nonlinear Optimization Problem With Caputo - Fabrizio 

Derivative, ITM Web Conf. 22 (2018), 01009. https://doi.org/10.1051/itmconf/20182201009. 

[30] A. Alvarado-Sánchez, E. Rangel-Cortes, A. Hernández-Hernández, Bi-Dimensional Crime Model Based on 

Anomalous Diffusion With Law Enforcement Effect, Math. Model. Numer. Simul. Appl. 2 (2022), 26–40. 

https://doi.org/10.53391/mmnsa.2022.01.003. 

[31] O. Abu Arqub. Series Solution of Fuzzy Differential Equations Under Strongly Generalized Differentiability, 

J. Adv. Res. Appl. Math. 5 (2013), 31-52. https://doi.org/10.5373/jaram.1447.051912. 

[32] M. Alaroud, M. Al-Smadi, R.R. Ahmad, et al. Computational Optimization of Residual Power Series 

Algorithm for Certain Classes of Fuzzy Fractional Differential Equations, Int. J. Differ. Equ. 2018 (2018), 

8686502. https://doi.org/10.1155/2018/8686502. 

[33] S. Hasan, M. Al-Smadi, H. Dutta, et al. Multi-Step Reproducing Kernel Algorithm for Solving Caputo–

Fabrizio Fractional Stiff Models Arising in Electric Circuits, Soft Comput. 26 (2022), 3713–3727. 

https://doi.org/10.1007/s00500-022-06885-4. 

[34] S. Momani, N. Djeddi, M. Al-Smadi, S. Al-Omari, Numerical Investigation for Caputo-Fabrizio Fractional 

Riccati and Bernoulli Equations Using Iterative Reproducing Kernel Method, Appl. Numer. Math. 170 

(2021), 418-434. https://doi.org/10.1016/j.apnum.2021.08.005. 

[35] M. Al-Smadi, O. Abu Arqub, S. Momani, Numerical Computations of Coupled Fractional Resonant 

Schrödinger Equations Arising in Quantum Mechanics Under Conformable Fractional Derivative Sense, 

Physica Scripta, 95 (2020), 075218. https://doi.org/10.1088/1402-4896/ab96e0. 

[36] A. Elsaid. Fractional Differential Transform Method Combined With the Adomian Polynomials, Appl. 

Math. Comput. 218 (2012), 6899-6911. https://doi.org/10.1016/j.amc.2011.12.066. 

https://doi.org/10.1186/s13662-021-03428-3
https://doi.org/10.1016/j.amc.2006.02.055
https://doi.org/10.1051/mmnp/2021030
https://doi.org/10.53391/mmnsa.2022.01.001
https://doi.org/10.12693/aphyspola.132.1066
https://doi.org/10.1051/itmconf/20182201009
https://doi.org/10.53391/mmnsa.2022.01.003
https://doi.org/10.5373/jaram.1447.051912
https://doi.org/10.1155/2018/8686502
https://doi.org/10.1007/s00500-022-06885-4
https://doi.org/10.1016/j.apnum.2021.08.005
https://doi.org/10.1088/1402-4896/ab96e0
https://doi.org/10.1016/j.amc.2011.12.066


15 Int. J. Anal. Appl. (2022), 20:51 
 

[37] A. Freihet, S. Hasan, M. Al-Smadi, et al. Construction of Fractional Power Series Solutions to 

Fractional Stiff System Using Residual Functions Algorithm, Adv. Differ. Equ. 2019 (2019), 95. 

https://doi.org/10.1186/s13662-019-2042-3. 

[38] M. Al-Smadi, O. Abu Arqub, S. Hadid. An Attractive Analytical Technique for Coupled System of 

Fractional Partial Differential Equations in Shallow Water Waves With Conformable Derivative, Commun. 

Theor. Phys. 72 (2020), 085001. https://doi.org/10.1088/1572-9494/ab8a29. 

[39] M. Al-Smadi, O. Abu Arqub, S. Hadid, Approximate solutions of nonlinear fractional Kundu-Eckhaus and 

coupled fractional massive Thirring equations emerging in quantum field theory using conformable residual 

power series method, Physica Scripta, 95 (2020), 105205. https://doi.org/10.1088/1402-4896/abb420. 

[40] M. Al-Smadi, A. Freihat, O. Abu Arqub, et al. A Novel Multistep Generalized Differential Transform 

Method for Solving Fractional-Order Lü Chaotic and Hyperchaotic Systems, J. Comput. Anal. Appl. 19 

(2015), 713–724. 

[41] M. Alabedalhadi, M. Al-Smadi, S. Al-Omari, et al. New Optical Soliton Solutions for Coupled Resonant 

Davey-Stewartson System With Conformable Operator, Opt. Quant. Electron. 54 (2022), 392. 

https://doi.org/10.1007/s11082-022-03722-8. 

[42] I. Podlubny, Fractional Differential Equation, Academic Press, San Diego, 1999. 

[43] M. Alabedalhadi, M. Al-Smadi, S. Al-Omari, et al. Structure of Optical Soliton Solution For Nonliear 

Resonant Space-Time Schrödinger Equation In Conformable Sense With Full Nonlinearity Term, Physica 

Scripta, 95 (2020), 105215. https://doi.org/10.1088/1402-4896/abb739. 

[44] Z. Altawallbeh, M. Al-Smadi, I. Komashynska, A. Ateiwi, Numerical Solutions of Fractional Systems of 

Two-Point BVPs by Using the Iterative Reproducing Kernel Algorithm, Ukr. Math J. 70 (2018), 687–701. 

https://doi.org/10.1007/s11253-018-1526-8. 

[45] A.A. Kilbas, H.M. Srivastava and J.J. Trujillo, Theory and Applications of Fractional Differential Equations, 

Elsevier Science B.V., Amsterdam, 2006. 

[46] S. Hasan, M. Al-Smadi, A. Freihet, et al.Two Computational Approaches for Solving a Fractional Obstacle 

System in Hilbert Space, Adv. Differ. Equ. 2019 (2019), 55. https://doi.org/10.1186/s13662-019-1996-

5. 

[47] S. Momani, O. Abu Arqub, A. Freihat, et al. Analytical Approximations for Fokker-Planck Equations of 

Fractional Order in Multistep Schemes, Appl. Comput. Math. 15 (2016), 319-330. 

[48] M. Al‐Smadi, O. Abu Arqub, M. Gaith, Numerical Simulation of Telegraph and Cattaneo Fractional‐Type 

Models Using Adaptive Reproducing Kernel Framework, Math. Meth. Appl. Sci. 44 (2020), 8472-8489. 

https://doi.org/10.1002/mma.6998. 

[49] R. Abu-Gdairi, M. Al-Smadi, G. Gumah, An Expansion Iterative Technique for Handling Fractional 

Differential Equations Using Fractional Power Series Scheme, J. Math. Stat. 11 (2015), 29-38. 

https://doi.org/10.3844/jmssp.2015.29.38. 

https://doi.org/10.1186/s13662-019-2042-3
https://doi.org/10.1088/1572-9494/ab8a29
https://doi.org/10.1088/1402-4896/abb420
https://doi.org/10.1007/s11082-022-03722-8
https://doi.org/10.1088/1402-4896/abb739
https://doi.org/10.1007/s11253-018-1526-8
https://doi.org/10.1186/s13662-019-1996-5
https://doi.org/10.1186/s13662-019-1996-5
https://doi.org/10.1002/mma.6998
https://doi.org/10.3844/jmssp.2015.29.38