Precision algorithms in second-order fractional differential equations

For a long time, we have focused on the mathematical theory of fractional calculus (FOC). It was not until the mid and late 20th century that FOC theory was introduced and applied in some fields such as electrochemistry, rheology, signal processing, and power transmission theory. Recently, the research of fractional order systems has aroused the interest of many scholars [1]. In the field of automatic control, a new branch of fractional control has emerged.

The research of fractional order control should cover at least three aspects: (1) The establishment of a fractional-order system (FOS) model and its analysis based on describing the fractional-order object more accurately and concisely. (2) The fractional-order control strategy is selected to obtain better control performance. (3) Apply fractional arithmetic to process signals, data, etc. The above work is based on the calculation of FOC and FOS.

Based on the Tustin transformation theory and the analysis of the characteristics of the discrete generating function formula used for fractional operators, and using the Maclaurin series expansion of the binomial power function, this paper proposes an improved method based on the power series expansion (PSE) and the Tustin transformation. Fractional calculation method [2]. We apply it to the solution of the linear fractional system and derive the corresponding recurrence formula. The performance of the algorithm's operation speed and accuracy is analyzed through specific example simulations and comprehensive comparisons with several conventional methods. The algorithm proposed in this paper is suitable for the fractional operation of arbitrary signals, including discrete data sequences whose mathematical model is unknown and the solution of linear systems.

We assume that the FOC of function f(t) is D^a f (t). Where a is the order. Z transformation is denoted as F(z), DF(z) respectively. Considering Eqs. (1) and (3), we can get (5) DF(z)F(z)=(ω(z−1))a=(2T(1−z−11+z−1))a=(2T)a(1−z−1)a(1+z−1)a {{DF(z)} \over {F(z)}} = {\left( {\omega \left( {{z^{ - 1}}} \right)} \right)^a} = {\left( {{2 \over T}\left( {{{1 - {z^{ - 1}}} \over {1 + {z^{ - 1}}}}} \right)} \right)^a} = {\left( {{2 \over T}} \right)^a}{{{{(1 - {z^{ - 1}})}^a}} \over {{{\left( {1 + {z^{ - 1}}} \right)}^a}}} But (6) (1+z−1)aDF(z)=(2T)a(1−z−1)aF(z) {(1 + {z^{ - 1}})^a}DF(z) = {\left( {{2 \over T}} \right)^a}{(1 - {z^{ - 1}})^a}F(z) We carry out the Maclaurin series expansion of the binomial power function (1 − z⁻¹)^a (take n terms to approximate) and write it in matrix form (7) (1−z−1)a≈∑j=0n(−1)j(aj)z−j=[p0,p1,⋯,⋯][1,z−1⋯z−j⋯]T=PZ {\left( {1 - {z^{ - 1}}} \right)^a} \approx \sum\limits_{j = 0}^n {( - 1)^j}\left( {\matrix{ a \cr j \cr } } \right){z^{ - j}} = \left[ {{p_0},{p_1}, \cdots , \cdots } \right]{\left[ {1,{z^{ - 1}} \cdots {z^{ - j}} \cdots } \right]^T} = PZ (aj) \left( {\matrix{ a \cr j \cr } } \right) is the binomial coefficient. P is the weight coefficient matrix. It is not difficult to get the recurrence relationship of the element pj=(−1)j(aj) {p_j} = ( - {1)^j}\left( {\matrix{ a \cr j \cr } } \right) of P by mathematical induction (8) p0=1,pj=(1−a+1j)pj−1 {p_0} = 1,{p_j} = \left( {1 - {{a + 1} \over j}} \right){p_{j - 1}} The same can be obtained (9) (1−z−1)a≈[q0,q1,⋯,qj,⋯][1,z−1⋯z−j⋯]=QZ \matrix{ {{{\left( {1 - {z^{ - 1}}} \right)}^a} \approx \left[ {{q_0},{q_1}, \cdots ,{q_j}, \cdots } \right]} \hfill \cr {\left[ {1,{z^{ - 1}} \cdots {z^{ - j}} \cdots } \right] = QZ} \hfill \cr } (10) q0=1, qj=(−1+a+1j)qj−1 {q_0} = 1,\quad {q_j} = \left( { - 1 + {{a + 1} \over j}} \right){q_{j - 1}} We assume that the discrete sequences of f(t) and D^af (t) are f(kT) and Df (kT). Abbreviated as f_k and df_k. k = 0, 1, 2,⋯ ,n. We substitute formula (7) and formula (9) into formula (6) to obtain d f_k + q₁d f_k−1 + ⋯ + q_k−1d f₁ + q_kd f₀ = (2/T)^a (f_k + p₁ f_k−1 + ⋯ + p_k−1 f₁ + p_k f₀). Then (11) dfk=(2/T)a(fk+p1fk−1+⋯+pk−1f1+pkf0)−(q1dfk−1+⋯+qk−1df1+qkdf0) d{f_k} = (2/T{)^a}\left( {{f_k} + {p_1}{f_{k - 1}} + \cdots + {p_{k - 1}}{f_1} + {p_k}{f_0}} \right) - \left( {{q_1}d{f_{k - 1}} + \cdots + {q_{k - 1}}d{f_1} + {q_k}d{f_0}} \right) Formula (11) is the FOC solution formula of f(t). The specific recurrence steps are as follows.

First consider the general form of the differential equation of a special fractional-order linear system: (12) a1Da1y(t)+a2Da2y(t)+⋯+anDany(t)=u(t) {a_1}{D^{a1}}y(t) + {a_2}{D^{a2}}y(t) + \cdots + {a_n}{D^{an}}y(t) = u(t) u(t) can be formed by a certain function and its fractional derivative. Considering formula (6), formula (7), formula (9), and formula (11), the approximate value d^a_iy_k of D^a_iy(t) can be obtained. We rewrite it as (13) daiy(t)≈daiyk=(2T)ai[yk+∑j=1kpj(ai)yk−j]−∑j=1kqj(ai)daiyk−j {d^{{a_i}}}y(t) \approx {d^{{a_i}}}{y_k} = {\left( {{2 \over T}} \right)^{{a_i}}}\left[ {{y_k} + \sum\limits_{j = 1}^k p{j^{({a_i})}}{y_{k - j}}} \right] - \sum\limits_{j = 1}^k q{j^{({a_i})}}{d^{{a_i}}}{y_{k - j}} pj^(a_i) and qj^(a_i) can be directly recursively obtained by formula (8) and formula (10). We substitute it into Eq. (12) and arrange it as ∑i=1nai{(2T)ai[yk+∑j=1kpj(ai)yk−j]−∑j=1kqj(ai)daiyk−j}=uk \sum\limits_{i = 1}^n {a_i}\left\{ {{{\left( {{2 \over T}} \right)}^{{a_i}}}\left[ {{y_k} + \sum\limits_{j = 1}^k p{j^{({a_i})}}{y_{k - j}}} \right] - \sum\limits_{j = 1}^k q{j^{({a_i})}}{d^{{a_i}}}{y_{k - j}}} \right\} = {u_k} . Then (14) yk=1∑i=1nai(2T)ai[uk−∑i=1nai(2T)ai∑j=1kqj(ai)yk−j+∑i=1nai∑j=1kqj(ai)dj(ai)yk−j] {y_k} = {1 \over {\sum\limits_{i = 1}^n {a_i}{{\left( {{2 \over T}} \right)}^{{a_i}}}}}\left[ {{u_k} - \sum\limits_{i = 1}^n {a_i}{{\left( {{2 \over T}} \right)}^{{a_i}}}\sum\limits_{j = 1}^k q{j^{({a_i})}}{y_{k - j}} + \sum\limits_{i = 1}^n {a_i}\sum\limits_{j = 1}^k q{j^{({a_i})}}d{j^{({a_i})}}{y_{k - j}}} \right] Equations (13) and (14) are recursive algorithm formulas that constitute the output of the fractional linear system. Specific steps are as follows.

First, to verify the effectiveness of the FOC recursive solution algorithm proposed in this paper, here is a simple simulation curve of the 0.5-order derivative, 0.5-order integral, and the first-order integral of the sine function f(t) = sint (Figure 5).

Fig. 5

It can be seen that the algorithm curve in this paper is the same as the result of the conventional PSE expansion approximation of the Tustin transform operator as a whole. This is almost the same as the curve obtained by the most direct and simple G-L definition discrete method [11]. The G-L method is almost identical to the theory obtained by applying the Cauchy definition to remove the area around t = 0⁺. Therefore, the results in Figure 5 show that the algorithm in this paper is feasible.

The unit step response is obtained through a specific FOS model. Based on the accuracy and speed, a comprehensive analysis and comparison between the method in this paper and the commonly used algorithms are made. The calculation and simulation tool is Matlab7. Let the FDEs of the simulation example be (15) D2y(t)+0.5D0.5y(t)+y(t)=D0.5u(t)+u(t) {D^2}y(t) + 0.5{D^{0.5}}y(t) + y(t) = {D^{0.5}}u(t) + u(t) For a class of FOS transfer function models that partial fractions can expand, the problem of obtaining theoretical numerical solutions has been well solved. We apply it to this example to easily find the analytical solution of its unit step output as (16) ystep(t)=L−1[G(s)U(s)]=L−1[∑k=14cks−1s0.5−λk]=∑k=14ckt0.5E0.5,0.5+1(λkt0.5) {y_{step}}(t) = {L^{ - 1}}\left[ {G(s)U(s)} \right] = {L^{ - 1}}\left[ {\sum\limits_{k = 1}^4 {{{c_k}{s^{ - 1}}} \over {{s^{0.5}} - {\lambda _k}}}} \right] = \sum\limits_{k = 1}^4 {c_k}{t^{0.5}}{E_{0.5,0.5 + 1}}\left( {{\lambda _k}{t^{0.5}}} \right) c_k, λ_k is the pole and residue of G (s^0.5) respectively. Take the error function (17) E=‖ystep*−ystep‖1/‖ystep*‖1 E = \left\| {y\mathop {step}\limits^* - ystep} \right\|1/\left\| {y\mathop {step}\limits^* } \right\|1 ystep* y\mathop {step}\limits^* is the comparison reference value, which is calculated by Eq. (16). ystep is the compared value. Considering that ystep* y\mathop {step}\limits^* reaches the steady-state value in about 20s, the length of the comparison data is 200 (step T = 0.1s). The precision index function is J = E^*/E. E^* is the reference value. Here we take the E value of the full memory G-L definition method (T = 0.1s). The index value has standard unity meaning. The simulation curves of each algorithm are shown in Figures 6–8. Table 1 shows the algorithm error and accuracy index values.

Fig. 6

Fig. 7

Fig. 8

The error is larger for step T = 0.1s (E = 0.862), and the smaller T is, the closer to the theoretical curve. Although high enough accuracy can be achieved when T = 0.001s is used, it is very difficult and too slow in terms of engineering realization. Considering the actual application of the project, the sampling period in the following simulations is selected to be easy to implement in the project.

Figure 6 shows that using the recursive discrete approximation algorithm given in this article to obtain the system's step response can obtain a relatively satisfactory result. The response curve of the memory length N = 100 almost coincides with the theoretical solution curve. Research shows (see Table 1) that the algorithm in this paper with an overall error of 10-2 orders of magnitude and full memory can achieve fairly high accuracy. The accuracy when N = 20 is equivalent to the result of the smaller step size T = 0.05s of the full memory G-L method. This verifies the conclusion that the response distortion of the Tustin discrete operator mentioned in the previous analysis is smaller than the Euler discrete operator. But when N > 100, the effect of accuracy improvement is not obvious [12]. This shows and verifies that the coefficients of the binomial power function can converge quickly. The above-mentioned excellent characteristics can be more clearly reflected in the comparison curve between the full memory algorithm in this paper and the G-L definition discrete method shown in Figure 7. The algorithm in this paper quickly converges to the theoretical value at the initial moment, and the accuracy advantage is also very significant.

It can be seen from Figure 8 that the initial response of the comparison curve between the algorithm in this paper and the conventional PSE-Tustin method is not much different. After 5 s, the difference is obvious, but this phenomenon is not reflected in Figure 5. The reason may be that different FOC discrete approximation methods exhibit different accuracy due to the difference in signal behavior. For this reason, we apply the fractional-order Al-Alaoui discrete operator to this example using the conventional PSE method and compare and analyze the results to show (Figure 9 and Table 1). The response curve shape of PSE-Al-Alaoui is similar to that of PSE-Tustin. In addition, the difference between the number of expansion terms of PSE-Tustin and the obvious influence on the accuracy is that the number of expansion terms has a less obvious influence on the accuracy. This shows that the convergence of the algorithm is better than that of PSE-Tustin. Overall, the results of the algorithm in this paper are the best. The PSE-Tustin method is better than the PSE-Al-Alaoui method.

Fig. 9

This paper proposes an improved PSE-Tustin recursive algorithm. This algorithm is suitable for the fractional operation of arbitrary signals, including discrete data sequences whose mathematical models are unknown and the solution of linear systems. A simulation example verifies the effectiveness of the algorithm. We show the characteristics and advantages of the algorithm in this paper by comparing the accuracy with several commonly used fractional-order algorithms.