# Discrete Waveform Relaxation Method for Linear Fractional Delay Differential-Algebraic Equations

Authors: Hongliang Liu; Yayun Fu; Bailing Li
Journal: Discrete Dynamics in Nature and Society
 (2017)

---

## Abstract

Fractional order delay differential-algebraic equations have the characteristics of time lag and memory and constraint limit. These yield some difficulties in the theoretical analysis and numerical computation. In this paper, we are devoted to solving them by the waveform relaxation method. The corresponding convergence results are obtained, and some numerical examples show the efficiency of the method.


---

## Body

## 1. Introduction

Fractional delay differential-algebraic equations are composed of fractional delay differential equations and algebraic equations, and they are more accurate in describing some scientific and engineering problems with memory function and algebraic restrictions. This kind of mathematical model is widely applied in many fields of electromagnetism, electrical systems, and biological materials. Nowadays, for the fractional differential equations, the main works focus on obtaining the analytic solution, approximately analytical solution, and numerical solution. For example, the Laplace transform method [1, 2], Fourier transform method [3, 4], and power method [5] are used to obtain the analytic solution. Furthermore, approximate analytical methods have been developed, such as variational iteration method [6, 7], Adomian decomposition method [8], homotopy analysis method [9], and multiple time scale method [10]. As for the studies of numerical methods, scholars discuss the linear multistep method [11], predictor corrector method [12, 13], Adomian method [14–16], spectrum method, and the finite element method [17, 18] for solving the fractional differential equations.The waveform relaxation (WR) method is proposed by Learsmee and so forth [19]. The main properties of the method are flexibility and good parallelism, which is used to solve differential equations [20, 21], fractional differential equations [22, 23], delay differential equations [24, 25], and differential-algebraic equations [26–28]. The waveform relaxation method based on multisplitting is proposed to improve parallelism [29, 30]. Furthermore, acceleration techniques are used to increase the speed of iteration, such as Krylov subspace acceleration, strange-type preconditioners, and windowing for the WR method [20, 21, 31, 32].The fractional differential-algebraic equations have received much attention; nevertheless, the numerical methods in this field are still young; a few studies have been considered on the convergence of the numerical methods, such as variational iteration method, Adomian decomposition method, and fractional differential transform method [33, 34]. Ding and Jiang applied the WR method to solve fractional differential-algebraic equations and obtain the convergence results [35]. But it is difficult to verify convergence conditions. The variational iteration method for the fractional delay integrodifferential-algebraic equations has been studied in [36], but it is not suitable for long-time numerical calculation. In order to overcome this limitation, according to the characteristics of the problems, the discrete WR method is employed to solve the linear fractional delay differential-algebraic equations by constructing the efficient iterative method and obtain the convergence results which is much easier to achieve.


## 2. Convergence

Consider theCaputo fractional delay differential-algebraic system(1)MD0αaCxt+Axt+Byt+Cxt-τ+Dyt-τ=Ut,t∈0,T,Ext+Nyt+Fxt-τ+Gyt-τ=Vt,t∈0,T,x0=y0=0,t∈0,T,where D∗αx(t) denotes the Caputo derivative of order α, n1-1<α≤n1,n1∈N, φ(t), ψ(t) are smooth vector functions, and φ(0)=0, ψ(0)=0, M∈Rn×n, N∈Rd×d are nonsingular matrix, A,C∈Rn×n, B,D∈Rn×d, E,F∈Rd×n, G∈Rd×d, and x(t)∈Cn[0,T], y(t)∈Rd, n+d=s, U(t)∈Rn, V(t)∈Rd are input function.Applying the WR method for the coefficient matrix of (1), we get the iteration form(2)M1D0αaCxk+1t+A1xk+1t+B1yk+1t+C1xk+1t-τ+D1yk+1t-τ=M2D0αaCxkt+A2xkt+B2ykt+C2xkt-τ+D2ykt-τ+Ut,E1xk+1t+N1yk+1t+F1xk+1t-τ+G1yk+1t-τ=E2xkt+N2ykt+F2xkt-τ+G2ykt-τ+Vt,xk+1t=xkt=⋯=x0t=φt,t∈-τ,0,yk+1t=ykt=⋯=y0t=ψt,t∈-τ,0,k=0,1,…,where A=A1-A2, B=B1-B2, C=C1-C2, D=D1-D2, E=E1-E2, F=F1-F2, G=G1-G2, M=M1-M2, N=N1-N2, and M1, N1 are nonsingular matrix. We assume that the step is h, τ=mh,m≥1, tp-j=tp-jh, p≥1, xp-j(k) is approximating value of x(k)(tp-j), and yp-j(k) is approximating value of y(k)(tp-j), j=0,1,…,p.Based on the fractionalCaputo derivative of Grünwald-Letnikov discrete format (3)D0αaCxt≐1hα∑j=0pωjxp-j,where ωj=(-1)jαj=(-1)jΓ(j-α)/Γ(-α)Γ(j+1) and α0=Γ(-α)/Γ(-α)Γ(1)=1, we obtain the iterative format (4)M1+hαA1xpk+1+∑j=1pωjM1xp-jk+1+hαv¯n=M2+hαA2xpk+∑j=1pωjM2xp-jk+hαu¯n,E1xpk+1+N1ypk+1+F1xp-mk+1+G1yp-mk+1=E2xpk+N2ypk+F2xp-mk+G2yp-mk+Vtp,xk+1t=xkt=⋯=x0t=φt,t∈-τ,0,yk+1t=ykt=⋯=y0t=ψt,t∈-τ,0,p=1,2,…,N,k=0,1,…,where(5)v¯n=B1ypk+1+C1xp-mk+1+D1yp-mk+1,u¯n=B2ypk+C2xp-mk+D2yp-mk+Utp,where N is the number of the time steps and Nh=T, the initial value x0(0)=0. Introduce (6)Xk=X1kT,X2kT,…,XpkTT,Yk=Y1kT,Y2kT,…,YpkTT.Next, we will discuss two cases.First Case. When p<m, (4) can be written as(7)H1hαI⊗B1I⊗E1I⊗N1Xk+1Yk+1=H2hαI⊗B2I⊗E2I⊗N2XkYk+Z^Q^,where Hi=Mα⊗Mi+hα(I⊗Ai),i=1,2, and(8)Mα=w0w1w0⋮⋱⋱wp-1⋱w1w0;introduce Z^=[z^1T,z^2T,…,z^iT,…,z^pT]T,Q^=[q^1T,q^2T,…,q^iT,…,q^pT]T, where (9)z^i=hαUti-Cφti-m-Dψti-m,i=1,2,…,p,q^i=Vti-Fφti-m-Gψti-m,i=1,2,…,p.Letδ1(k+1)=[(X(k+1)-X(k))T,(Y(k+1)-Y(k))T]T,k=0,1,….From (7), we have(10)δ1k+1=H1hαI⊗B1I⊗E1I⊗N1-1·H2hαI⊗B2I⊗E2I⊗N2δ1k.Let the iterative matrix beR1 in (10); the convergence condition of the iterative process (4) is ρ(R1)<1. Moreover, we have the following theorem.Theorem 1.
For the small steph and p<m, if (11)ρM1-1M2<1,ρN1-1N2<1,then the discrete waveform relaxation iteration method (4) is convergent.Proof.
The iteration process (4) convergence is decided by R¯1. Firstly, the left coefficient matrix of (7) can be written as (12)H1hαI⊗B1I⊗E1I⊗N1=Mα⊗M10I⊗E1I⊗N1+hαI⊗A1I⊗B100,because of the following equation: (13)Mα⊗M10I⊗E1I⊗N1-1=Mα⊗M1-10W1I⊗N1-1,where W1=-(I⊗N1)-1(I⊗E1)(Mα⊗M1)-1; then, (14)H1hαI⊗B1I⊗E1I⊗N1-1=Mα⊗M1-10W1I⊗N1-1-1+Ohα,and we can get(15)R1=Mα⊗M1-10W1I⊗N1-1+Ohα·Mα⊗M20I⊗E2I⊗N2+hαI⊗A2I⊗B200=Mα⊗M1-1Mα⊗M20W~1I⊗N1-1I⊗N2+Ohα,where W~1=W1(Mα⊗M2)+(I⊗N1)-1(I⊗N2); we have(16)ρMα⊗M1-1Mα⊗M2=ρM1-1M2.Based on the matrix spectral radius definition of R1 expression in (15), if (17)ρM1-1M2<1,ρN1-1N2<1,and h is small enough, we have ρ(R1)<1.
Therefore, whenp<m, the convergence is proved.Second Case. When p>m, we assume p=m+c; then, system (4) can be written as(18)H-1N¯1K¯1L¯1Xk+1Yk+1=H-2N¯2K¯2L¯2XkYk+S^T^,where (19)H-i=Mα⊗Mi+hαI⊗Ai+P¯ξ⊗Ci,i=1,2,N¯i=hαI⊗Bi+P¯ξ⊗Di,i=1,2,K¯i=I⊗Ei+P¯ξ⊗Fi,i=1,2,L¯i=I⊗Ni+P¯ξ⊗Gi,i=1,2,P¯ξ=000⋮⋱⋱100⋮⋱⋱⋱0⋯1⋯00,where (20)P¯ξ=ξij=1,i=m,m+1,…,N,j=i-m+1,0,others;introduce S^=[s^1T,s^2T,…,s^iT,…,s^pT]T,T^=[t^1T,t^2T,…,t^iT,…,t^pT]T, where s^i,t^i are (21)s^i=hαUti-Cφti-m-Dψti-m,i=1,2,…,m-1,hαUti,i=m,m+1,…,p,t^i=Vti-Fφti-m-Gψti-m,i=0,…,m-1,Vti,i=m,m+1,…,p.Letδ~1(k+1)=[(X(k+1)-X(k))T,(Y(k+1)-Y(k))T]T,k=0,1,…, and from (18) we have(22)δ~1k+1=H-1N¯1K¯1L¯1-1·H-2N¯2K¯2L¯2δ~1k.Let the iterative matrix beR2 in (22); then, the convergence condition of the iterative process (22) is ρ(R2)<1; we have the following theorem.Theorem 2.
For the small steph, if it satisfies (23)ρM1-1M2<1,ρN1-1N2<1,then the discrete waveform relaxation iteration method (18) is convergent.Proof.
The iteration process (18) convergence is decided by R2. Firstly, the left coefficient matrix of (18) can be written as (24)H-1N¯1K¯1L¯1=Mα⊗M10K¯1L¯1+hαI⊗A1+P¯ξ⊗C1I⊗B1+P¯ξ⊗D100,and because (25)L¯1-1=N1-10N1-1⋮⋱⋱-N1-1G1N1-10N1-1⋱⋱⋱0-N1-1G1N1-1⋯0N1-1=I⊗N1-1+00L~10,where L~1=-P¯ξ⊗(-N1-1G1N1-1), hence (26)Mα⊗M10K¯1L¯1-1=Mα⊗M1-10W2L¯1-1,where W2=-L1-1(I⊗E1+P¯ξ⊗F1)(Mα⊗M1)-1; then,(27)H-1N¯1K¯1L¯1-1=Mα⊗M1-10W^2I⊗N1-1+Ohα, where W-2=W2+L~1; we have(28)R2=Mα⊗M1-10W2L¯1-1+Ohα·Mα⊗M20K¯2L¯2+hαI⊗A2+P¯ξ⊗C2I⊗B2+P¯ξ⊗D200=Mα⊗M1-1Mα⊗M20W~2L¯1-1L¯2+Ohα,where W~2=W-2(Mα⊗M2)+L¯1-1L¯2, and because (29)L¯1-1L¯2=N1-1N20N1-1N2⋮⋱⋱-N1-1G1N1-1G20N1-1N2⋱⋱⋱0-N1-1G1N1-1G2⋯0N1-1N2=I⊗N1-1I⊗N2+00L~20,where L~2=-P¯ξ⊗(N1-1G1N2-1G2), therefore, we have(30)ρL¯1-1L¯2=ρN1-1N2;from (8), we have (31)ρMα⊗M1-1Mα⊗M2=ρM1-1M2.
Based on the matrix spectral radius definition and combiningR2 expression of (28) and if (32)ρM1-1M2<1,ρN1-1N2<1,and h is enough small, we have ρ(R2)<1.
Therefore, the convergence is proved.


## 3. Illustrative Examples

In this section, some illustrative examples are given to show the efficiency of the WR method for solving the linear fractional delay differential-algebraic equations.Example 1.
Consider the initial value problem of the linear fractional delay differential-algebraic equations(33)-510150000D0αaCxt+-0.10000.30.50.0010-0.5xt+0.120.2000.1-0.5000.8xt-0.5=sin⁡tsin⁡t0,t∈0,T,xt=0,0,0T,t∈-0.5,0,where α=0.5, x(t)=(x1(t),x2(t),x3(t))T; introduce x~(t)=(x1(t),x2(t))T,y(t)=x3(t); system (33) can be written as(34)MD0αaCx~t+Ax~t+Byt+Cx~t-0.5+Dyt-0.5=Ut,t∈0,T,Ex~t+Nyt+Fx~t-0.5+Gyt-0.5=Vt,t∈0,T,x~t=0,0T,yt=0,t∈-0.5,0,and we have (35)M=-5115,A=-0.1000.3,B=00.5,C=0.120.200.1,D=0-0.5,E=0.0010,N=-0.5,F=00,G=0.8,Ut=sin⁡tsin⁡tT,Vt=0.
Splitting the coefficient matrix, we can construct theJacobi waveform relaxation method and deduce (36)M1=-5005,M2=01-10,A1=-0.100-0.3,A2=0000,B1=00.5,B2=00,C1=0.12000.1,C2=00.200,D1=0-0.5,D2=00,E1=0.0010,E2=00,F1=F2=00,N1=-0.51.1,N2=-0.051.1,G1=0.8,G2=0.
According to the iterative format (4), we take time step h=0.1 and use the Jacobi waveform relaxation method to solve the initial value problem (33). Then, we take time step h=0.005 and take the numerical solution of (33) obtained by the implicit Gru¨nwald-Letnikov method as the approximate true solution. Then, the error of the numerical solution obtained by WR method is shown in Table 1.Table 1
The error of numerical solutions.
T
T
=
1
T=2
T
=
3
T
=
4
x~1(t)
3.0981
e
-
003
1.1308
e
-
003
2.3020e-003
3.8279
e
-
003
x~2(t)
4.1365
e
-
003
1.1548e-003
3.5419e-003
5.4337
e
-
003
y1(t)
1.9550e-005
3.4345
e
-
004
1.6464
e
-
003
5.455
e
-
003The numerical solutions are plotted in Figure1, where ∗, ⋄, and ○, respectively, are numerical solutions of x~1(t), x~2(t), and y1(t) by using the waveform relaxation.Figure 1
Numerical solutions by using the two methods for solving (33).Table1 and Figure 1 show that the waveform relaxation method for obtaining the solution is accurate.Table2 shows that the iterative allowable error is smaller, and the number of corresponding iterations is higher. If the iteration tolerance is the same, then the condition of iterations number is relatively stable.Table 2
Between the iteration allowable error and the number of iterations.
T
1.0
×
10
-
4
1.0
×
10
-
6
1.0
×
10
-
8
T
=
1
5
7
10
T
=
2
4
7
10
T
=
3
4
7
10
T
=
4
5
8
10Example 2.
Consider the initial value problem of the linear fractional delay differential-algebraic equations(37)5100260000000000D0αaCxt+4210.5150.81.20.10.3620.20.227xt+31.2411.452610320114xt-0.5=tsin⁡t1-cos⁡t,t∈0,T,xt=0,0,0,0T,t∈-0.5,0,where α=0.5, x(t)=(x1(t),x2(t),x3(t),x4(t))T; introduce x~(t)=(x1(t),x2(t))T, y(t)=(x3(t),x4(t))T; system (37) can be written as(38)MD0αaCx~t+Ax~t+Byt+Cx~t-0.5+Dyt-0.5=Ut,t∈0,T,Ex~t+Nyt+Fx~t-0.5+Gyt-0.5=Vt,t∈0,T,x~t=0,0T,yt=0,0T,t∈-0.5,0,and we have (39)M=5126,A=4215,B=10.50.82,C=31.21.45,D=4126,E=0.10.30.20.2,N=6227,F=1001,G=3214,Ut=tsin⁡tT,Vt=1-cos⁡t.Splitting the coefficient matrix, we can construct the Gauss-Seidel waveform relaxation method and deduce (40)M1=5026,M2=0100,A1=4015,A2=0200,B1=100.81.2,B2=00.500,C1=301.45,C2=01.200,D1=4026,D2=0100,E1=0.10.30.20.2,E2=0000,N1=6027,N2=0200,F1=1001,F2=0000,G1=3014,G2=0200.
According to the iterative format (4), we take time step h=0.1 and use theGauss-Seidel waveform relaxation method to solve the initial value problem (37). Then, we take time step h=0.005 and take the numerical solution of (37) obtained by the implicit Gru¨nwald-Letnikov method as the approximate true solution. Then, the error of the numerical solution obtained by WR method is shown in Table 3.
The numerical solutions are plotted in Figure2, where ∗, ⋄, and ○, respectively, are numerical solutions of x~1(t), x~2(t), and y1(t) by using the waveform relaxation.
Table3 and Figure 2 show that the waveform relaxation method for obtaining the solution is accurate.
Table4 shows that the iterative allowable error is smaller, and the number of corresponding iterations is higher. If the iteration tolerance is the same, then the condition of iterations number is relatively stable.Table 3
The error of numerical solutions.
T
T
=
1
T=2
T
=
3
T
=
4
x~1(t)
4.1571
e
-
003
4.9982
e
-
003
4.7782
e
-
003
3.4877
e
-
003
x~2(t)
1.2799
e
-
003
1.2010
e
-
003
3.2493
e
-
003
3.0521
e
-
003
y1(t)
4.1742
e
-
003
4.1408
e
-
003
1.8243
e
-
003
1.2053
e
-
003
y2(t)
1.0734
e
-
003
1.2718
e
-
003
7.4072
e
-
004
1.1129
e
-
004Table 4
Between the iteration allowable error and the number of iterations.
T
1.0
×
10
-
4
1.0
×
10
-
6
1.0
×
10
-
8
T
=
1
4
5
7
T
=
2
3
5
7
T
=
3
4
6
7
T
=
4
4
6
7Figure 2
Numerical solutions by using the two methods for solving (34).


## 4. Conclusion

In this paper, the discrete waveform relaxation method has been extended to solve linear fractional delay differential equations. Some illustrative examples show that this method gives very good approximation to the exact solution. The discrete waveform relaxation method is a promising method to handle linear fractional delay differential-algebraic equations.



---

*Source: 6306570-2017-06-15.xml*