Next Article in Journal
Existence Results for Coupled Nonlinear Sequential Fractional Differential Equations with Coupled Riemann–Stieltjes Integro-Multipoint Boundary Conditions
Next Article in Special Issue
A Numerical Method for Simulating Viscoelastic Plates Based on Fractional Order Model
Previous Article in Journal
Right Fractional Sobolev Space via Riemann–Liouville Derivatives on Time Scales and an Application to Fractional Boundary Value Problem on Time Scales
Previous Article in Special Issue
Adaptive Fractional Image Enhancement Algorithm Based on Rough Set and Particle Swarm Optimization
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Numerical Algorithm for Calculating the Time Domain Response of Fractional Order Transfer Function

1
School of Information Engineering, Shenyang University, Shenyang 110044, China
2
School of Information Science and Engineering, Northeastern University, Shenyang 110819, China
*
Author to whom correspondence should be addressed.
Fractal Fract. 2022, 6(2), 122; https://doi.org/10.3390/fractalfract6020122
Submission received: 20 January 2022 / Revised: 13 February 2022 / Accepted: 17 February 2022 / Published: 19 February 2022

Abstract

:
This paper proposes new numerical algorithms for calculating the time domain responses of fractional order transfer functions (FOTFs). FOTFs are divided into two categories, explicit fractional order transfer functions (EFOTFs) and implicit fractional order transfer functions (IFOTFs). Transforming an EFOTF into an equivalent fractional order differential equation, its time domain response can be obtained by solving the equation by the difference method. IFOTF cannot be transformed into an equivalent equation, so its time domain response cannot be calculated by existing difference methods. A new numerical algorithm is designed for calculating a convolution and its inverse operation, the time domain response of IFOTF can be calculated based on the algorithm. Error analysis shows that the proposed numerical algorithms are of first-order accuracy. Four calculation examples are presented, and the results are consistent with the theoretical analysis.

1. Introduction

Fractional calculus is a theory about the integration and differentiation of non-integer order. It is a powerful tool to describe systems which have long-term memory and long-range spatial interaction. Fractional calculus has been applied in control theory in recent years, and there have been many theoretical achievements in controllability [1,2], observability [3,4], stability [5,6,7], robustness [8,9] and control methods [10,11,12]. At the same time, fractional order theory has also achieved great success in control engineering, such as the fractional order PID controller [13,14], air-based precision positioning system [15], fractional order positive position feedback compensator [16], and fractional order gray model [17,18]. More achievements can be found in [19].
As many novel fractional order transfer functions (FOTFs) have emerged in the above applications, an urgent problem is how to calculate their time domain responses. It is usually difficult to find the analytical expression of the time domain response; therefore, numerical algorithms are used to solve this problem. If a FOTF can be transformed into an equivalent fractional order differential equation, its time domain response can be calculated by solving the equation with a numerical algorithm; such a FOTF is called explicit fractional order transfer function (EFOTF) in this paper. There are many numerical algorithms for solving fractional order differential equations, such as the linear multi-step algorithm [20,21,22], separation of variables algorithm [23], predictor–corrector algorithm [24,25,26], matrix method [27] and some recently proposed algorithms [28,29,30]. In addition, the time domain response of EFOTF can also be calculated by MATLAB toolbox, FOTF [31].
If a FOTF cannot be transformed into an equivalent fractional order differential equation, its time domain response cannot be calculated by the difference algorithm; such a FOTF is called implicit fractional order transfer function (IFOTF) in this paper. Some frequency methods are effective for some simple IFOTFs, for example, the time domain response of the following IFOTF can be calculated by the method proposed in [32].
F ( s ) = 1 1 + s P T α ,
where P T is the transitional frequency, and α is a real number between 0 and 1. However, it is difficult to estimate the calculation error, and frequency methods are ineffective for more complex IFOTFs, such as
F ( s ) = 340 s 0.756 ( s 2 + 3.85 s + 5880 ) 1.15 .
The IFOTF is a model of ionic polymer metal composite (IPMC) proposed in [33]; more information can be found in [34,35]. Therefore, a time domain algorithm is needed, which can calculate the time response of the IFOTF directly.
The main contribution of this work is to design a numerical algorithm for the IFOTF. The time domain response can be calculated by the proposed algorithm without considering the specific form of the IFOTF. The error analysis is presented, and it is proved that the proposed algorithm is of first-order accuracy. For the IFOTFs which can only be analyzed by frequency methods before, the proposed algorithm can give their time domain characteristics directly. The IFOTF is analyzed in both the frequency domain and time domain, so the result is more reliable.
Section 1 is devoted to expound the background and the significance of this work. The definition and the notation are introduced in Section 2. The first algorithm is proposed for EFOTF in Section 3. The second algorithm is designed for IFOTF based on the first algorithm in Section 4. The error analysis is presented in Section 5. Section 6 consists of four calculation examples. Section 7 contains the conclusion and the future work.

2. Definition

There are different fractional order definitions; Riemann–Liouville (RL) and Caputo definitions are used in this paper.
Definition 1
([36]). The RL fractional order derivative is defined as
t 0 RL D t α f ( t ) = 1 Γ ( n α ) d n d t n t 0 t f ( τ ) ( t τ ) α n + 1 d τ ,
where t 0 is an arbitrary real constant, α is a real constant bigger than zero, n = α , variable t [ t 0 , + ) , Γ is the Gamma function.
Definition 2
([37]). The Caputo fractional order derivative is defined as
t 0 C D t α f ( t ) = 1 Γ ( n α ) t 0 t f ( n ) ( τ ) ( t τ ) α n + 1 d τ .
The parameters in (3) and (4) have the same meaning.
In (3) and (4), t 0 denotes the initial time. If t 0 0 , it can always be converted to zero by variable substitution; therefore, t 0 is always equal to zero.
The transfer function refers to the ratio of the Laplace transform of the response (output) to the Laplace transform of the excitation (input) under the zero initial condition. Therefore, it is always assumed that the function has a zero initial condition when discussing the transfer function. Under the zero initial condition, Definitions 1 and 2 are equivalent [36]; they are written as 0 D t α f ( t ) , i.e.,
0 D t α f ( t ) = 0 RL D t α f ( t ) = 0 C D t α f ( t ) .
FOTFs are divided into two categories, EFOTF and IFOTF.
Definition 3.
The explicit fractional order transfer function (EFOTF) is defined as
F ( s ) = b m s β m + b m 1 s β m 1 + + b 1 s β 1 a n s α n + a n 1 s α n 1 + + a 1 s α 1 ,
where s is a complex variable, coefficients a 1 , , a n , b 1 , , b m are all real constants, orders α n > > α 1 0 , β m > > β 1 0 are all nonnegative real constants, and α n β m .
Definition 4.
The implicit fractional order transfer function (IFOTF) is defined as
F ( s ) = b m s β m + b m 1 s β m 1 + + b 1 s β 1 μ 2 / ν 2 a n s α n + a n 1 s α n 1 + + a 1 s α 1 μ 1 / ν 1 ,
where μ 1 , ν 1 , μ 2 , ν 2 are all positive integers, other parameters have the same meaning as the parameters in (6).
The topic of this paper is how to calculate the time domain response of (6) and (7). The details of the problem are described below. The time interval [ 0 , T ] is divided into N identical subintervals and the step size h = T / N . The value of the excitation u ( t ) at the grid point is written as u k ( 0 k N ) . The same notation is applied to the response y ( t ) , written as y k ( 0 k N ) . The problem is how to calculate the response y k of (6) and (7), under some given excitation u k .

3. The Numerical Algorithm for EFOTF

The numerator and denominator of (6) are denoted as V ( s ) and W ( s ) , and the Laplace transforms of excitation u ( t ) and response y ( t ) are written as U ( s ) and Y ( s ) . EFOTF (6) is written as
W ( s ) Y ( s ) = V ( s ) U ( s ) .
According to the convolution theorem, the inverse Laplace transform of (8) is
0 t w ( τ ) y ( t τ ) d τ = 0 t v ( τ ) u ( t τ ) d τ ,
where v ( t ) and w ( t ) are the inverse Laplace transforms of W ( s ) and V ( s ) , respectively. The trapezoidal rule is employed to calculate the two integrals in (9),
h j = 1 k 1 w j y k j + w 0 y k + w k y 0 2 = h j = 1 k 1 v j u k j + v 0 u k + v k u 0 2 + O ( h 2 ) .
The recurrence formula is obtained from (10),
y k = 2 w 0 j = 1 k 1 ( v j u k j w j y k j ) + v 0 u k + v k u 0 w k y 0 w 0 + O ( h ) .
Formula (11) can be used to calculate y k , as long as w j and v j are known. The following method is proposed for calculating w j and v j .
EFOTF (6) is equivalent to the following fractional order differential equation,
i = 1 n a i 0 D t α i y ( t ) = i = 1 m b i 0 D t β i u ( t ) ,
with zero initial condition, its solution is the response of (6).
The fractional linear multistep method [21,22] is employed to calculate the fractional order derivative in (12),
0 D t α i y ( t ) | t = k h = 1 h α i j = 0 k ω j ( α i , p ) y k j + O ( h p ) ,
where the coefficient ω j ( α i , p ) is calculated by the following recursive formula,
ω j ( α i , p ) = 0 , j < 0 , θ 0 α i , j = 0 , 1 θ 0 k = 1 p θ k 1 k 1 + α i j ω j k ( α i , p ) , j > 0 .
The calculation error of (13) is O ( h p ) , and the proof can be found in [21].
Plug (13) into (12). It yields
i = 0 n a i h α i j = 0 k ω j ( α i , p ) y k j = i = 0 m b i h β i j = 0 k ω j ( β i , p ) u k j + O ( h p ) .
Multiply both sides of (15) by step size h, transpose the order of summation,
h j = 0 k i = 0 n a i h α i ω j ( α i , p ) · y k j = h j = 0 k i = 0 m b i h β i ω j ( β i , p ) · u k j + O ( h p + 1 ) .
Set p = 1 in (16) and compare (10) and (16). It yields the recurrence formula for calculating w j and v j ,
w j = i = 0 n a i h α i ω j ( α i , 1 ) + O ( h ) , v j = i = 0 m b i h β i ω j ( β i , 1 ) + O ( h ) , when 0 < j < k , w j = 2 i = 0 n a i h α i ω j ( α i , 1 ) + O ( h ) , v j = 2 i = 0 m b i h β i ω j ( β i , 1 ) + O ( h ) , when j = 0 o r k .
The above calculation process is summarized as the numerical algorithm.
Algorithm 1.Calculating the time domain response of EFOTF (6).
1. 
Calculate w j , v j by (17).
2. 
Use (11) to calculate y k , the time domain response at the grid point.

4. The Numerical Algorithm for IFOTF

There is no equation equivalent to IFOTF (7), so its time domain response cannot be calculated by Algorithm 1. A new numerical algorithm is proposed for IFOTF in this section.
The numerator and denominator of (7) are denoted as V ^ ( s ) and W ^ ( s ) , and the Laplace transforms of excitation u ( t ) and response y ( t ) are written as U ( s ) and Y ( s ) . IFOTF (7) is written as
W ^ ( s ) Y ( s ) = V ^ ( s ) U ( s ) .
Take the ν 1 ν 2 th power of both sides of (18),
W ^ ν 1 ν 2 ( s ) Y ν 1 ν 2 ( s ) = V ^ ν 1 ν 2 ( s ) U ν 1 ν 2 ( s ) ,
where W ^ ν 1 ν 2 ( s ) and V ^ ν 1 ν 2 ( s ) are polynomials. Consider U ν 1 ν 2 ( s ) as the excitation and Y ν 1 ν 2 ( s ) as the response. An equivalent EFOTF is obtained,
Y ν 1 ν 2 ( s ) U ν 1 ν 2 ( s ) = i = 1 n a i s α i μ 1 i = 1 m b i s β i μ 2 .

4.1. Numerical Convolution

According to the convolution theorem, the inverse Laplace transform of U ν 1 ν 2 ( s ) is
L 1 { U ν 1 ν 2 ( s ) } = u ( t ) u ( t ) u ( t ) ν 1 ν 2 .
The notation ∗ denotes a convolution,
u ( t ) u ( t ) = 0 t u ( t τ ) u ( τ ) d τ .
The value of u ( j h ) is written as u j , and h is the step size. Assume that the first-order error is contained in u j , i.e.,
u j = u ( j h ) + O ( h ) , ( 0 j N ) .
Apply the trapezoidal rule to calculate the value of (22),
0 k h u ( k h τ ) u ( τ ) d τ = h j = 0 k u k j u j h 2 ( u 0 u k + u k u 0 ) + O ( h 2 ) j = 0 k u j + O ( h 2 ) = h j = 0 k u k j u j + O ( h ) , ( 0 k N ) .
Although u j contains O ( h ) , (24) is still of first-order accuracy. Therefore, the right side of (21) can be calculated by applying (24) repeatedly. The result remains of first-order accuracy.
The generating polynomial is introduced to simplify the expression. The generating polynomial of u j is defined as
u ¯ ( x ) = j = 0 N u j x j ,
where u j is defined as (23). The square of u ¯ ( x ) is
u ¯ 2 ( x ) = k = 0 2 N j = 0 k u k j u j x k .
Comparing (24) and (26), it yields that the result of (24) is the coefficient of (26) multiplied by step size h. The right side of (21) can be calculated by this method. It is expressed as,
u ( t ) u ( t ) u ( t ) ν 1 ν 2 | t = k h = h ν 1 ν 2 1 S k u ¯ ν 1 ν 2 ( x ) + O ( h ) , ( 0 k N ) ,
where S k represents selecting the coefficient of x k .
The result of (27) is the value of L 1 { U ν 1 ν 2 ( s ) } at grid point. It is the excitation of (20). Because (20) is an EFOTF, its response can be calculated by Algorithm 1. The result is written as R k ( 0 k N ) . The result is also the value of L 1 { Y ν 1 ν 2 ( s ) } at grid point, so it yields the equation,
y ( t ) y ( t ) y ( t ) ν 1 ν 2 | t = k h = R k + O ( h ) , ( 0 k N ) .

4.2. Numerical Inversion of Convolution

The solution of (28) is written as y j ( 0 j N ) , and its generating polynomial is
y ¯ ( x ) = j = 0 N y j x j .
Calculate the left side of (28) by the same method as (27) and ignore the small quantity term O ( h ) . The following equation is obtained:
h ν 1 ν 2 1 S k y ¯ ν 1 ν 2 ( x ) = R k , ( 0 k N ) .
Let x = 0 and k = 0 in (30). It yields
y 0 = R 0 h 1 ν 1 ν 2 1 ν 1 ν 2 .
It is supposed that y 0 , , y k 1 are calculated, and thus, y k can be calculated by the following method. The polynomial y ¯ ( x ) is written as the sum of three parts,
y ¯ ( x ) = y ¯ l ( x ) + y ¯ e ( x ) + y ¯ h ( x ) = j = 0 k 1 y j x j + y k x k + j = k + 1 N y j x j .
According to the binomial theorem,
y ¯ ν 1 ν 2 ( x ) = y ¯ l ( x ) + y ¯ e ( x ) ν 1 ν 2 + j = 1 ν 1 ν 2 ν 1 ν 2 j y ¯ h j ( x ) y ¯ l ( x ) + y ¯ e ( x ) ν 1 ν 2 j = j = 0 ν 1 ν 2 ν 1 ν 2 j y ¯ e j ( x ) y ¯ l ν 1 ν 2 j ( x ) + j = 1 ν 1 ν 2 ν 1 ν 2 j y ¯ h j ( x ) y ¯ l ( x ) + y ¯ e ( x ) ν 1 ν 2 j = y ¯ l ν 1 ν 2 ( x ) + ν 1 ν 2 y ¯ l ν 1 ν 2 1 ( x ) y ¯ e ( x ) + j = 2 ν 1 ν 2 ν 1 ν 2 j y ¯ e j ( x ) y ¯ l ν 1 ν 2 j ( x ) + j = 1 ν 1 ν 2 ν 1 ν 2 j y ¯ h j ( x ) y ¯ l ( x ) + y ¯ e ( x ) ν 1 ν 2 j .
In (33), only the first and the second terms contain x k , B k 1 denotes the coefficient of x k contained in the first term, and B k 2 denotes the coefficient of x k contained in the second term. Plug B k 1 and B k 2 into (30). It yields,
B k 1 + B k 2 = R k h 1 ν 1 ν 2 .
Expand the second term of (33),
ν 1 ν 2 y ¯ l ν 1 ν 2 1 ( x ) y ¯ e ( x ) = ν 1 ν 2 y k x k j = 0 k 1 y j x j ν 1 ν 2 1 = ν 1 ν 2 y k y 0 ν 1 ν 2 1 x k + ν 1 ν 2 y k x k j = 1 k 1 y j x j ν 1 ν 2 1 .
The second term of (35) only contains the term whose order is higher than k, so B k 2 is equal to the coefficient of the fist term in (35),
B k 2 = ν 1 ν 2 y 0 ν 1 ν 2 1 y k .
Lemma 1 is cited from [38] for expanding the first term of (33).
Lemma 1
([38]). Let c be a positive integer. For all x 1 , x 2 , , x d ,
( x 1 + x 2 + + x d ) c = c c 1 c 2 c d x 1 c 1 x 2 c 2 x d c d ,
where the summation extends over all nonnegative integer solutions c 1 , c 2 , , c d of c 1 + c 2 + + c d = c , and the multinomial coefficients are defined by
c c 1 c 2 c d = n ! c 1 ! c 2 ! c d ! .
Expand the first term of (33) by Lemma 1. The coefficient of x k is
B k 1 = ν 1 ν 2 c 0 c 1 c k 1 y 0 c 0 y 1 c 1 y k 1 c k 1 .
The summation extends over all nonnegative integer solutions c 1 , c 2 , , c k 1 , which satisfy the following equations,
j = 0 k 1 c j = ν 1 ν 2 , j = 0 k 1 j c j = k .
According to (39), B k 1 can be calculated directly with y 0 , y 1 , , y k 1 , which are calculated. Plug (36) and (39) into (34). It yields the formula for calculating y k ,
y k = R k h 1 ν 1 ν 2 B k 1 ν 1 ν 2 y 0 ν 1 ν 2 1 .
The above calculation process is summarized as the numerical algorithm.
Algorithm 2.Calculating the time domain response of IFOTF (7).
1. 
Transform IFOTF (7) into EFOTF (20).
2. 
Calculate the time domain excitation of (20) by (27).
3. 
Solve the time domain response of (20) by Algorithm 1; obtain equation (30).
4. 
Calculate y 0 by (31), and calculate y k ( 0 k N ) by (41).

5. Error Analysis

Two formulas, (11) and (17), are employed in Algorithm 1; the parameters in (11) are calculated by (17). It yields that Algorithm 1 is of first-order accuracy by plugging (17) into (11).
Algorithm 2 consists of three parts. In the first part, IFOTF (7) is transformed into EFOTF (20), the excitation of (20) is calculated by (27), which is of first-order accuracy. The generating polynomial is introduced solely for the sake of convenience in expression; it does not produce calculation error. Therefore, the result of this part is of first-order accuracy.
In the second part, the response of EFOTF (20) is calculated by Algorithm 1; it is of first-order accuracy. It is worth noting that the excitation of EFOTF (20) is used in this part; it contains the first-order error. This error does not affect the calculation accuracy of (11); therefore, the result of the second part is of first-order accuracy.
In the last part, (28) is solved by (41), and solution y k is the time domain response of IFOTF (7). Use p to denote the order of error, i.e.,
y ( k h ) = y k + O ( h p ) .
Plug y k into (30). The calculation is of pth-order accuracy if p < 1 ; otherwise, the calculation is of first-order accuracy.
Assume p < 1 and plug (42) into the left side of (30). The result is of pth-order accuracy, i.e.,
h ν 1 ν 2 1 S k y ¯ ν 1 ν 2 ( x ) = R k + O ( h p ) .
It is contradictory to (28) because (28) and (43) have the same solution. It yields p 1 .
The solution y k is calculated by (41). The numerator of (41) contains R k , which is calculated by Algorithm 1. Algorithm 1 is of first-order accuracy, so the calculation accuracy of y k cannot be higher than the first order, i.e., p 1 . Combined with the preceding analysis, it yields p = 1 ; Algorithm 2 is of first-order accuracy.

6. Calculation Examples

Four examples are presented in this section, the calculation results are used to validate the proposed algorithms.

6.1. Example 1

An EFOTF is expressed as
Y ( s ) = 1 s 0.7 + s 0.5 U ( s ) ,
where Y ( s ) and U ( s ) are the Laplace transforms of response and excitation, respectively. The inverse Laplace transform of U ( s ) is
u ( t ) = Γ ( 1.8 ) Γ ( 1.1 ) t 0.1 + Γ ( 1.8 ) Γ ( 1.3 ) t 0.3 .
It can be verified that the time domain response is
y ( t ) = t 0.8 .
Use Algorithm 1 to calculate the time domain response in [ 0 , 10 ] . Select step size h = 0.01 ; the analytical and numerical solutions are shown in Figure 1. There are two curves in Figure 1, and because the calculation error is very small, the two curves coincide.
Use (46) to test the accuracy of the numerical solution. The computation error is defined as | y ( j h ) y j | , and the computation errors under different step sizes are reported in Table 1. The error decreases along with decreasing the step size h, and they are proportional to each other. This result is consistent with the theoretical analysis, Algorithm 1 is of first-order accuracy.

6.2. Example 2

Algorithm 1 can also be used to solve fractional order differential equations. This example was previously completed by the author of this paper; only the problem and result are cited, and the calculation process can be found in [39].
The equation is proposed in [26],
0 C D t 1.455 y ( t ) = t 0.1 E 1 , 1.545 ( t ) E 1 , 1.445 ( t ) e t y ( t ) 0 C D t 0.555 y ( t ) + e 2 t 0 C D t 1 y ( t ) 2 ,
where E 1 , 1.545 ( t ) and E 1 , 1.445 ( t ) are two-parameter Mittag–Leffler functions. The analytical solution is y ( t ) = e t under the initial condition y ( 0 ) = 1 , y ( 0 ) = 1 .
Solve (47) by Algorithm 1 and the PECE-type Adams algorithm proposed in [26]. The maximal error and the execution time are reported in Table 2; the result of the PECE-type Adams algorithm is cited from [26], and the result of Algorithm 1 is cited from [39]. Considering the difference of computer hardware, it is of little significance to compare the execution time. The maximal error is an important parameter for testing an algorithm. It is obvious that Algorithm 1 is more accurate; its maximal error is proportional to the step size.

6.3. Example 3

An IFOTF is expressed as
Y ( s ) = 1 ( 4 s + 1 ) 0.5 U ( s ) ,
where Y ( s ) and U ( s ) are the Laplace transforms of response and excitation, respectively. The inverse Laplace transform of U ( s ) is
u ( t ) = t 2 .
The time domain response is
y ( t ) = 0.5 e t 4 0 D t 0.5 t 2 e t 4 .
The analytical expression of (50) can be verified by the formula introduced in [40]. An IFOTF is expressed as
Y ( s ) = 1 ( τ s + 1 ) ν U ( s ) ,
and the analytical expression of the time domain response is
y ( t ) = τ ν e t τ 0 D t ν u ( t ) e t τ .
Use Algorithm 2 to calculate the time domain response in [ 0 , 10 ] . Select step size h = 0.01 ; the analytical and numerical solutions are shown in Figure 2. There are two curves in Figure 2; because the computation error is very small, the two curves coincide.
Use (50) to test the accuracy of the numerical solution. The computation error is defined as | y ( j h ) y j | , computation errors under different step sizes are reported in Table 3. The error is proportional to step size, the result is consistent with the theoretical analysis, Algorithm 2 is of first-order accuracy.

6.4. Example 4

IFOTF (2) is a model for the ionic polymer metal composite (IPMC) actuator [33], which is mentioned at the beginning of the paper. Its time domain response is calculated in this example. Rewrite its expression
Y ( s ) = 340 s 0.756 ( s 2 + 3.85 s + 5880 ) 1.15 U ( s ) ,
where Y ( s ) and U ( s ) are the Laplace transforms of response and excitation, respectively. The expression of U ( s ) is
U ( s ) = Γ ( 8 ) ( s + 1 ) 8 ,
the inverse Laplace transform of U ( s ) is
u ( t ) = t 7 e t .
The expression of Y ( s ) is
Y ( s ) = 340 Γ ( 8 ) ( s + 1 ) 8 s 0.756 ( s 2 + 3.85 s + 5880 ) 1.15 .
Use Algorithm 2 to calculate the time domain response in [ 0 , 20 ] . The time domain response is calculated by another method to verify the result of Algorithm 2. The inverse Laplace transform of (56) is calculated by the numerical algorithm proposed in [41]. Select step size h = 0.01 ; the results of the two algorithms are plotted in Figure 3. The two curves almost coincide, and it verifies that the two results are credible. It is worthy noting that (56) is used in the numerical inversion of the Laplace transform, but not used in Algorithm 2. That means that the expression of the excitation is not necessary for Algorithm 2; it is a more effective method in practical application.
The difference is defined as the result of Algorithm 2 minus the result of the numerical inversion of the Laplace transform. The difference is reported in Table 4; the difference decreases along with decreasing the step size. It shows that the two results are credible, although neither of them is an exact solution.

7. Discussion

Two numerical algorithms are proposed in this paper. The novelty is that the time domain response can be calculated without considering the form of FOTF. Although frequency methods are also effective for some simple IFOTFs, they are ineffective for more complex IFOTFs. The proposed algorithms are useful tools for analyzing FOTFs, and also for verifying the results of the frequency methods. The FOTF is analyzed in both the frequency domain and time domain; thus, the result is more reliable. Four examples are presented in the paper; the calculation results are consistent with the theoretical analysis. Examples 1 and 3 are simple problems for which the analytical expressions of the time domain responses are provided. The computation error can be obtained with the analytical expression, so the two examples can be used as benchmark problems for testing numerical algorithms; they are helpful for researchers to test their own numerical algorithms. In Example 3, Algorithm 1 is compared with the PECE-type Adams algorithm, and the result shows that Algorithm 1 is more accurate. In Example 4, the time domain response of (2) is calculated by Algorithm 2. The IFOTF (2) is the model of IPMC; the authors of this paper have not found other numerical algorithms to calculate its time domain response. If the parameters of Algorithm 1 are changed, a more accurate solution can be obtained; however, Algorithm 2 is of first-order accuracy. In future work, further studies will be carried out to improve the calculation accuracy of Algorithm 2.

Author Contributions

Conceptualization, L.B. and D.X.; methodology, L.B.; writing—original draft preparation, L.B.; writing—review and editing, L.B. and D.X. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Data Availability Statement

Not applicable.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Liu, B.; Su, H.S.; Wu, L.C.; Li, X.L.; Lu, X. Fractional-order controllability of multi-agent systems with time-delay. Neurocomputing 2021, 424, 268–277. [Google Scholar] [CrossRef]
  2. Ahmad, I.; Rahman, G.U.; Ahmad, S.; Alshehri, N.A.; Elagan, S.K. Controllability of a damped nonlinear fractional order integrodifferential system with input delay. Alex. Eng. J. 2022, 61, 1956–1966. [Google Scholar] [CrossRef]
  3. Shamardan, A.B.; Moubarak, M.R.A. Controllability and observability for fractional control systems. J. Fract. Calc. 1999, 15, 25–34. [Google Scholar]
  4. Hassanzadeh, I.; Tabatabaei, M. Calculation of controllability and observability matrices for special case of continuous-time multi-order fractional systems. ISA Trans. 2018, 82, 62–72. [Google Scholar] [CrossRef]
  5. Zhang, X.F.; Zhao, Z.L. Normalization and stabilization for rectangular singular fractional order T-S fuzzy systems. Fuzzy Sets Syst. 2020, 381, 140–153. [Google Scholar] [CrossRef]
  6. Zhang, X.F. Relationship between integer order systems and fractional order systems and its two applications. IEEE-CAA J. Autom. Sin. 2018, 5, 639–643. [Google Scholar] [CrossRef]
  7. Zhang, X.F.; Chen, Y.Q. Admissibility and robust stabilization of continuous linear singular fractional order systems with the fractional order α: The 0 < α < 1 case. ISA Trans. 2018, 82, 42–50. [Google Scholar]
  8. Chen, Y.Q.; Ahn, H.S.; Xue, D.Y. Robust controllability of interval fractional order linear time invariant systems. Signal Process. 2006, 86, 2794–2802. [Google Scholar] [CrossRef]
  9. Zhang, J.X.; Yang, G.H. Robust Adaptive Fault-Tolerant Control for a Class of Unknown Nonlinear Systems. IEEE Trans. Ind. Electron. 2017, 64, 585–594. [Google Scholar] [CrossRef]
  10. Zhang, X.F.; Huang, W.K. Adaptive neural network sliding mode control for nonlinear singular fractional order systems with mismatched uncertainties. Fractal Fract. 2020, 4, 50. [Google Scholar] [CrossRef]
  11. Zhang, J.X.; Yang, G.H. Low-complexity tracking control of strict-feedback systems with unknown control directions. IEEE Trans. Autom. Control 2019, 64, 5175–5182. [Google Scholar] [CrossRef]
  12. Zhang, J.X.; Yang, G.H. Fuzzy adaptive output feedback control of uncertain nonlinear systems with prescribed performance. IEEE Trans. Cybern. 2018, 48, 1342–1354. [Google Scholar] [CrossRef]
  13. Mandić, P.D.; Šekara, T.B.; Lazarević, M.P.; Bošković, M. Dominant pole placement with fractional order PID controllers: D-decomposition approach. ISA Trans. 2017, 67, 76–86. [Google Scholar] [CrossRef]
  14. Mandić, P.D.; Lazarević, M.P.; Šekara, T.B. D-Decomposition technique for stabilization of furuta pendulum: Fractional approach. Bull. Pol. Acad. Sci. Tech. Sci. 2016, 64, 189–196. [Google Scholar] [CrossRef] [Green Version]
  15. Krijnen, M.E.; van Ostayen, R.; HosseinNia, S.H. The application of fractional order control for an air-based contactless actuation system. ISA Trans. 2018, 82, 172–183. [Google Scholar] [CrossRef] [Green Version]
  16. Marinangeli, L.; Alijani, F.; HosseinNia, S.H. Fractional-order positive position feedback compensator for active vibration control of a smart composite plate. J. Sound Vib. 2018, 412, 1–16. [Google Scholar] [CrossRef]
  17. Yang, Y.; Xue, D.Y. A novel actual load forecasting by interval grey modeling methodology based on the fractional calculus. ISA Trans. 2018, 82, 200–209. [Google Scholar] [CrossRef]
  18. Yang, Y.; Xue, D.Y. Continuous fractional-order grey model and electricity prediction research based on the observation error feedback. Energy 2016, 115, 722–733. [Google Scholar] [CrossRef]
  19. Sun, H.G.; Yong, Z.; Baleanu, D.; Chen, W.; Chen, Y.Q. A new collection of real world applications of fractional calculus in science and engineering. Commun. Nonlinear Sci. Numer. Simul. 2018, 64, 213–231. [Google Scholar] [CrossRef]
  20. Lubich, C. On the stability of linear multistep methods for Volterra convolution equations. IMA J. Numer. Anal. 1983, 3, 439–465. [Google Scholar] [CrossRef]
  21. Lubich, C. Discretized fractional calculus. SIAM J. Math. Anal. 1986, 17, 704–719. [Google Scholar] [CrossRef]
  22. Lubich, C. Fractional linear multistep methods for Abel-Volterra integral equations of the first kind. IMA J. Numer. Anal. 1987, 7, 97–106. [Google Scholar] [CrossRef]
  23. Daftardar-Gejji, V.; Bhalekara, S. Boundary value problems for multi-term fractional differential equations. J. Math. Anal. Appl. 2008, 345, 754–765. [Google Scholar] [CrossRef] [Green Version]
  24. Diethelm, K.; Ford, N.J.; Freed, A.D. A predictor-corrector approach for the numerical solution of fractional differential equations. Nonlinear Dyn. 2002, 29, 3–22. [Google Scholar] [CrossRef]
  25. Diethelm, K. An investigation of some nonclassical methods for the numerical approximation of Caputo-type fractional derivatives. Numer. Algorithms 2008, 47, 361–390. [Google Scholar] [CrossRef]
  26. Diethelm, K. The Analysis of Fractional Differential Equations; Springer: Berlin, Germany, 2010; pp. 229–248. [Google Scholar]
  27. Podlubny, I. Matrix approach to discrete fractional calculus. Fract. Calc. Appl. Anal. 2000, 3, 359–386. [Google Scholar]
  28. Wei, Y.Q.; Liu, D.Y.; Driss, B.; Chen, Y.M. An improved pseudo-state estimator for a class of commensurate fractional order linear systems based on fractional order modulating functions. Syst. Control. Lett. 2018, 118, 29–34. [Google Scholar] [CrossRef]
  29. Chen, Y.M.; Liu, L.Q.; Liu, D.Y.; Boutat, D. Numerical study of a class of variable order nonlinear fractional differential equation in terms of Bernstein polynomials. Ain Shams Eng. J. 2018, 9, 1235–1241. [Google Scholar] [CrossRef] [Green Version]
  30. Liu, D.Y.; Tian, Y.; Boutat, D.; Laleg-Kirati, T. An algebraic fractional order differentiator for a class of signals satisfying a linear differential equation. Signal Process. 2015, 116, 78–90. [Google Scholar] [CrossRef] [Green Version]
  31. Xue, D.Y. Fractional-Order Control Systems-Fundamentals and Numerical Implementations; De Gruyter: Berlin, Germany, 2017; pp. 1–30. [Google Scholar]
  32. Mansouri, R.; Bettayeb, M.; Djennoune, S. Approximating of high order integer systems by fractional order reduced parameter models. Math. Comput. Model. 2010, 51, 53–62. [Google Scholar] [CrossRef]
  33. Caponetto, R.; Dongola, G.; Fortuna, L.; Graziani, S.; Strazzeri, S. A fractional model for IPMC actuators. In Proceedings of the 2008 IEEE Instrumentation and Measurement Technology Conference, Victoria, BC, Canada, 20 June 2008. [Google Scholar]
  34. Caponetto, R.; Graziani, S.; Pappalardo, F.; Umana, E.; Xibilia, M.G.; Giamberardino, P. A scalable fractional order model for IPMC actuators. IFAC Proc. Vol. 2012, 45, 593–596. [Google Scholar] [CrossRef]
  35. Caponetto, R.; Graziani, S.; Pappalardo, F.; Xibilia, M.G. Parametric Control of IPMC Actuator Modeled as Fractional Order System. Adv. Sci. Technol. 2013, 79, 63–68. [Google Scholar]
  36. Podlubny, I. Fractional Differential Equations: An Introduction to Fractional Derivatives, Fractional Differential Equations, to Methods of Their Solution and Some of Their Applications; Academic Press: Cambridge, MA, USA, 1998; pp. 41–119. [Google Scholar]
  37. Caputo, M. Linear Models of dissipation whose Q is almost frequency independent. Geophys. J. R. Astr. Soc. 1967, 13, 529–539. [Google Scholar] [CrossRef]
  38. Brualdi, R.A. Introductory Combinatorics, 5th ed.; Pearson Education: New York, NY, USA, 2009; pp. 127–160. [Google Scholar]
  39. Xue, D.Y.; Bai, L. Numerical algorithms for Caputo fractional-order differential equations. Int. J. Control 2017, 90, 1201–1211. [Google Scholar] [CrossRef]
  40. Sabatier, J.; Farges, C. Analysis of fractional models physical consistency. J. Vib. Control 2017, 23, 895–908. [Google Scholar] [CrossRef]
  41. Valsa, J.; Brancik, L. Approximate formulae for numerical inversion of Laplace transforms. Int. J. Numer. Model. 1998, 11, 153–166. [Google Scholar] [CrossRef]
Figure 1. The time domain response of EFOTF (44).
Figure 1. The time domain response of EFOTF (44).
Fractalfract 06 00122 g001
Figure 2. The time domain response of IFOTF (48).
Figure 2. The time domain response of IFOTF (48).
Fractalfract 06 00122 g002
Figure 3. The time domain response of IFOTF (2).
Figure 3. The time domain response of IFOTF (2).
Fractalfract 06 00122 g003
Table 1. Computation error of the time domain response of EFOTF (44).
Table 1. Computation error of the time domain response of EFOTF (44).
Step Size t = 2 t = 4 t = 6 t = 8 t = 10
h = 0.1 8.2728 × 10 3 8.4603 × 10 3 8.4762 × 10 3 8.3949 × 10 3 8.3039 × 10 3
h = 0.05 4.7671 × 10 3 4.7630 × 10 3 4.6350 × 10 3 4.5479 × 10 3 4.4700 × 10 3
h = 0.01 1.1865 × 10 3 1.1491 × 10 3 1.0730 × 10 3 1.0384 × 10 3 1.0109 × 10 3
h = 0.005 6.3320 × 10 4 6.0817 × 10 4 5.6145 × 10 4 5.4124 × 10 4 5.2541 × 10 4
h = 0.001 1.4169 × 10 4 1.3431 × 10 4 1.2169 × 10 4 1.1655 × 10 4 1.1261 × 10 4
Table 2. Comparison of the errors and the execution times.
Table 2. Comparison of the errors and the execution times.
Standard PECE-Type Adams AlgorithmAlgorithm 1
Step SizeMaximal ErrorExecution TimeMaximal ErrorExecution Time
h = 1 / 200 0.3904101.2 s0.012500.35667 s
h = 1 / 400 0.2193368.4 s0.006810.43396 s
h = 1 / 800 0.11641358.0 s0.003580.99232 s
h = 1 / 1600 0.06005017.4 s0.001852.54765 s
Table 3. Computation error of the time domain response of IFOTF (48).
Table 3. Computation error of the time domain response of IFOTF (48).
Step Size t = 2 t = 4 t = 6 t = 8 t = 10
h = 0.1 3.9420 × 10 2 6.3127 × 10 2 1.2113 × 10 1 1.4719 × 10 1 1.6519 × 10 1
h = 0.05 1.9809 × 10 2 3.1679 × 10 2 6.0705 × 10 2 7.3727 × 10 2 8.2712 × 10 2
h = 0.01 3.9728 × 10 3 6.3512 × 10 3 1.2162 × 10 2 1.4766 × 10 2 1.6560 × 10 2
h = 0.005 1.9870 × 10 3 3.1764 × 10 3 6.0825 × 10 3 7.3841 × 10 3 8.2810 × 10 3
h = 0.001 3.9748 × 10 4 6.3543 × 10 4 1.2167 × 10 3 1.4770 × 10 3 1.6564 × 10 3
Table 4. Difference between the results of the two algorithms.
Table 4. Difference between the results of the two algorithms.
Step Size t = 4 t = 8 t = 12 t = 16 t = 20
h = 0.2 3.55 × 10 1 4.70 × 10 1 2.26 × 10 2 8.80 × 10 2 5.97 × 10 2
h = 0.1 1.76 × 10 1 2.36 × 10 1 1.11 × 10 2 4.41 × 10 2 3.00 × 10 2
h = 0.05 8.76 × 10 2 1.18 × 10 1 5.55 × 10 3 2.21 × 10 2 1.51 × 10 2
h = 0.02 3.48 × 10 2 4.71 × 10 2 2.31 × 10 3 8.96 × 10 3 6.12 × 10 3
h = 0.01 1.73 × 10 2 2.35 × 10 2 1.24 × 10 3 4.56 × 10 3 3.14 × 10 3
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Bai, L.; Xue, D. Numerical Algorithm for Calculating the Time Domain Response of Fractional Order Transfer Function. Fractal Fract. 2022, 6, 122. https://doi.org/10.3390/fractalfract6020122

AMA Style

Bai L, Xue D. Numerical Algorithm for Calculating the Time Domain Response of Fractional Order Transfer Function. Fractal and Fractional. 2022; 6(2):122. https://doi.org/10.3390/fractalfract6020122

Chicago/Turabian Style

Bai, Lu, and Dingyü Xue. 2022. "Numerical Algorithm for Calculating the Time Domain Response of Fractional Order Transfer Function" Fractal and Fractional 6, no. 2: 122. https://doi.org/10.3390/fractalfract6020122

APA Style

Bai, L., & Xue, D. (2022). Numerical Algorithm for Calculating the Time Domain Response of Fractional Order Transfer Function. Fractal and Fractional, 6(2), 122. https://doi.org/10.3390/fractalfract6020122

Article Metrics

Back to TopTop