In order to see the difference between DAEs and different classes of DDAEs more clearly, we summarize the results about the strangeness-free formulation and the consistency initial condition for DAEs as well as for DDAEs in Tables5.1,5.2.
System type Regular, strangeness-free formulation restricted to [(i−1)τ,iτ], 16i 6` DAEs
•Eˆ1(t) 0
‚
˙ x(t)−
•Aˆ1(t) Aˆ2(t)
‚ x(t)=
•fˆ1(t) fˆ2(t)
‚ , where£Eˆ1T AˆT2⁄
is pointwise nonsingular.
Causal DDAEs
•Eˆ1(t) 0
‚
˙ x(t)=
•Aˆ1(t) Aˆ2(t)
‚ x(t)+
•Bˆ0,1(t) Bˆ0,2(t)
‚
x(t−τ) +
µ
P
i=1
• 0 Bˆi,2(t)
‚
x(i)(t−τ)+
•fˆ1(t) fˆ2(t)
‚ , where£Eˆ1T AˆT2⁄
is pointwise nonsingular.
General retarded DDAEs
•Eˆi,1(t) 0
‚ x(t)˙ =
•Aˆi,1(t) Aˆi,2(t)
‚ x(t)+
•Bˆi,1(t) 0
‚
x(t−τ)+
•gˆi,1(t) ˆ gi,2(t)
‚ , whereh
EˆTi,1 AˆTi,2i
is pointwise nonsingular.
General neutral DDAEs
•Eˆi,1(t) 0
‚
˙ x(t)=
•Aˆi,1(t) Aˆi,2(t)
‚ x(t)+
•Bˆi,1(t) Bˆi,2(t)
‚
x(t−τ)+
•gˆi,1(t) ˆ gi,2(t)
‚ , where
hEˆTi,1 AˆTi,2i
is pointwise nonsingular.
General advanced DDAEs The strangeness-free formulation takes the form (5.14).
The underlying DDE takes the form Pµ
α=0Aˆαi (t)x(α)(t)=Pη+µ−1
α=0 Bˆiα(t)x(α)(t−τ)+gˆi(t),
whereAˆµi is pointwise nonsingular and there exists at least one functionBˆiα,i≥µ+1, that is not identically zero.
Table 5.1:Strangeness-free formulation for different classes of DDAEs.
5.3. Regularization of DDAEs by operational arrays 74 System type Consistency of initial conditions
DAEs −Aˆ2(0)x(0)=fˆ2(0).
Causal DDAEs
8
>>
<
>>
:
0 =Aˆ2(0)φ(0)+
µ
P
i=1
Bˆi,2(0)φ(i)(−τ)+fˆ2(0),
0 =
µ
P
i=0
Bˆi,3(t)φ(i)(t−τ)+fˆ3(t), for allt∈(0,τ).
General retarded DDAEs
(Aˆ1,2(0)φ(0)+gˆ1,2(0) =0,
Bˆ1,3(t)φ(t−τ)+gˆ1,3(t) =0, for allt∈(0,τ).
General neutral DDAEs
(Aˆ1,2(0)φ(0)+Bˆ1,2(0)φ(−τ)+gˆ1,2(0) =0,
Bˆ1,3(t)φ(t−τ)+Bˆ1,4(t) ˙φ(t−τ)+gˆ1,3(t) =0, for allt∈(0,τ).
General advanced DDAEs 0=
η
P
i=0
Bˆi1,µ+2(t)φ(i)(t−τ)+fˆµ+21 (t), for allt∈(0,τ).
Table 5.2:Consistency of initial conditions for different classes of DDAEs.
ii) The differentiation operation that maps the DDAE (1.2a) into the equation d
dt(E(t) ˙x(t)−A(t)x(t))= d dt
¡B(t)x(t−τ)+f(t)¢ .
Making use of the differentiation operation and combining equations with their deriva-tives to derive a strangeness-free formulation or an underlying ODE is a well-known technique in the theory of DAEs, [27, 28]. In this section, under Assumptions 5.1, 5.13, we study the numerical reformulation for the DDAE (1.2a) in order to achieve the DDAE (5.20) pointwise. Consider one arbitrary pointt ∈(t0,tf). For k ∈N, dif-ferentiating the DDAE (1.2a)k times we obtain thederivative arraysordifferentiation inflated system
M(t)z(t)=P(t)z(t−τ)+g(t), (5.22) where
M := 2 6 6 6 6 6 6 4
−A E
−A˙ E˙−A E
−A¨ E¨−2 ˙A 2 ˙E−A E
... . .. . ..
−A(k) E(k)−k A(k−1) . . . kE˙−A E 3 7 7 7 7 7 7 5 ,
P := 2 6 6 6 6 6 6 4
B 0
B˙ B 0
B¨ 2 ˙B B 0
... . .. ... ... B(k) kB(k−1) . . . kB˙ B 0
3 7 7 7 7 7 7 5
, z:= 2 6 6 6 4
x
˙ x ... x(k+1)
3 7 7 7 5
, g:= 2 6 6 6 4
f f˙ ... f(k+1)
3 7 7 7 5 .
In the following by the subscriptjτ, j∈Z, we denote the evaluation of a function at the pointt+jτ. In particular, we will writeτinstead of 1τ. Assuming that the shift index
κis well-defined and that it is a constant for allt ∈I, the shift-inflated system (5.4) is then given by
E0τx˙0τ = A0τx0τ+B0τx−τ+f0τ, Eτx˙τ = Aτxτ+Bτx0τ+fτ, E2τx˙2τ = A2τx2τ+B2τxτ+f2τ,
...
Eκτx˙κτ = Aκτxκτ+Bκτx(κ−1)τ+fκτ.
(5.23)
In order to build the derivative array of system (5.23), one needs to rewrite (5.23) as 2
6 6 6 6 6 4
E0τ Eτ
E2τ . ..
Eκτ 3 7 7 7 7 7 5 2 6 6 6 6 6 4
˙ x0τ
˙ xτ
˙ x2τ
...
˙ xκτ
3 7 7 7 7 7 5
= 2 6 6 6 6 6 4
A0τ Bτ Aτ
B2τ A2τ . .. ...
Bκτ Aκτ 3 7 7 7 7 7 5 2 6 6 6 6 6 4
x0τ xτ x2τ
... xκτ
3 7 7 7 7 7 5 +
2 6 6 6 6 6 4
B0τx−τ+f0τ fτ f2τ
... fκτ
3 7 7 7 7 7 5 .
(5.24) Let us suppose that for the DAE (5.24), the strangeness index ˆµ(t) is well-defined and becomes a constant in a sufficiently small neighborhoodB(t) oft. Now we can build the derivative arrays withk=µˆ(t) for the DAE (5.24). However, to reduce the cost in the determination of the shift indexκ, it would be better to build the derivative arrays (with k =µ(t)) for the equations of system (5.23) instead of building the derivative arraysˆ for the entire system (5.24). Following this way, we have the so-calleddouble-inflated system
2 6 6 6 6 6 4
M0τ
−Pτ Mτ
−P2τ M2τ
. .. . ..
−Pκτ Mκτ 3 7 7 7 7 7 5 2 6 6 6 6 6 4
z0τ zτ z2τ
... zκτ
3 7 7 7 7 7 5
= 2 6 6 6 6 6 4
P0τ 0 0 ... 0
3 7 7 7 7 7 5
z−τ+ 2 6 6 6 6 6 4
g0τ gτ g2τ
... gκτ
3 7 7 7 7 7 5
. (5.25)
We denote the matrix coefficients of (5.25) as
M := 2 6 6 6 6 6 4
M0τ
−Pτ Mτ
−P2τ M2τ . .. . ..
−Pκτ Mκτ 3 7 7 7 7 7 5
, P := 2 6 6 6 6 6 4
P0τ
0 0 ... 0
3 7 7 7 7 7 5
, G:= 2 6 6 6 6 6 4
g0τ
gτ g2τ
... gκτ
3 7 7 7 7 7 5
. (5.26)
Here one should not confuse P with the matrix polynomial P in Subsection 4.3.2.
Within this section, we do not consider any matrix polynomial. We discuss now how to derive the regular, strangeness-free DDAE (5.20) from the double-inflated system (5.25).
Remark 5.18. It should be noted that due to the non-causality of the DDAE (1.2a), the differential equations of the regular, strangeness-free DDAE (5.20) must be selected from the double-inflated system (5.25), instead of from the original DDAE (1.2a). This is in contrast to both cases of non-delayed DAEs and of causal DDAEs. We illustrate
5.3. Regularization of DDAEs by operational arrays 76 this fact in the next example.
Example 5.19. Consider the IVP consisting of the DDAE
•1 0 0 0
‚ •x(t)˙
˙ y(t)
‚
=
•0 0 1 0
‚ •x(t) y(t)
‚ +
•0 0 0 −1
‚ •x(t−τ) y(t−τ)
‚ +
• et
−et+t−τ
‚
, (5.27)
for t∈[0,∞),τ=1, with an initial functionφ(t) :=
•et t
‚
, for t∈[−τ, 0].
By directly checking, we obtainκ=1and the regular, strangeness-free DDAE(5.35)is
•0 −1
0 0
‚ •x(t)˙
˙ y(t)
‚
=
•0 0 1 0
‚ •x(t) y(t)
‚ +
•0 0 0 −1
‚ •x(t−τ) y(t−τ)
‚ +
• −1
−et+t−τ
‚ .
Clearly, the differential equation 0=£
0 −1⁄
•x(t˙ )
˙ y(t)
‚
−1 cannot be selected from the original DDAE(5.27).
From Remark5.18, we see that it is necessary to select all the constraints forx(t) and ˙x(t) contained in (5.25). It is worth to note that x(t) and ˙x(t) are only present inz0τ but not inzτ, . . . ,zκτ. In the following, for notational convenience, we will use Matlab notation, [91].
Let the matrixU be such that its columns span the space corangeM(:, (2n+1) :end), i. e.,
UTM(:, (2n+1) :end)=0. (5.28) By scaling (5.25) withUT, we obtain the system
UTM(:, 1 : 2n)
•x(t)
˙ x(t)
‚
=UTPz−τ+UTG, (5.29) that contains all the constraints forx(t) and ˙x(t) in (5.25).
We further notice that in the DDAE (5.20) onlyx(t−τ) occurs even thoughz−τcontains not onlyx(t−τ) but also its derivatives ˙x(t−τ), ¨x(t−τ), . . . . Thus, the matrixU should be chosen to satisfy the additional condition
UTP(:, (n+1) :end)=0. (5.30) In the case thatUTP(:, (n+1) :end)6=0, it implies that the DDAE (1.2a) is of advanced type, and hence Assumption5.13is violated.
Denote by
M˜:=UTM(:, (n+1) : 2n), N˜:=UTM(:, 1 :n), ˜P:=UTP(:, 1 :n), ˜G:=UTG, (5.31a) and let ˆmbe the number of rows of the matrix functionM˜. ThusM˜,N˜, ˜P are matri-ces of size ˆmbyn. We consider the following spaces and matrices
Z2 basis of corange( ˜M)=ker( ˜MT), T2 basis of ker(Z2TN˜), Z1 basis of range( ˜MT2).
(5.31b)
On the one hand, the algebraic constraints of the DDAE (5.20) are contained in the resulting system obtained by scaling system (5.29) withZ2T as
Z2TN˜x(t)=Z2TP˜x(t−τ)+Z2TG˜. (5.32) On the other hand, the set of differential equations in the DDAE (5.20) are contained in the resulting system obtained by scaling system (5.29) withZ1T as
Z1TM˜x(t˙ )+Z1TN˜x(t)=Z1TP˜x(t−τ)+Z1TG˜. (5.33) Thus, the fact that the regular, strangeness-free DDAE (5.20) is contained in the cou-pled system (5.32), (5.33) implies that
rank
•Z1TM˜ Z2TN˜
‚
=n.
Furthermore, we still need to remove some equations in (5.32) and (5.33), due to the fact thatZ2TN˜ andZ1TM˜may not have full row rank. Therefore, we consider the fol-lowing spaces and matrices
Y2 basis of range(Z2TN˜),
Y1 basis of range(Z1TM˜). (5.34) Scaling (5.32) (resp., (5.33)) withY2T (resp.Y1T), we see that the DDAE (5.20) becomes
•Y1TZ1TM˜ 0
‚
˙ x(t)+
•Y1TZ1TN˜ Y2TZ2TN˜
‚ x(t)=
•Y1TZ1TP˜ Y2TZ2TP˜
‚
x(t−τ)+
•Y1TZ1TG˜ Y2TZ2TG˜
‚
, (5.35)
where
•Y1TZ1TM˜ Y2TZ2TN˜
‚
is nonsingular.
In summary, the regularization of the IVP (1.2) by using operational arrays is pro-posed in Algorithm5.2below.
Algorithm 5.2Regularization of the IVP (1.2)
1: Input:The IVP (1.2).
2: Return:The regular, strangeness-free DDAE (5.20) pointwise.
3: Letκ=0.
4: Construct the double-inflated system (5.25) with the coefficientsM, P, Gas in (5.26).
5: Let the matrixUbe such that the conditions (5.28), (5.30) hold.
6: Compute the matricesM˜,N˜,Z2,T2,Z1,Y1,Y2as in (5.31) and (5.34).
7: if
rank
•Z1TM˜ Z2TN˜
‚¶
=n,
then(5.35) is exactly the regular, strangeness-free DDAE (5.20).
8: else κ:=κ+1, go back to 4.
9: end if
5.3. Regularization of DDAEs by operational arrays 78 Let us demonstrate Algorithm5.2by studying the following example.
Example 5.20. Consider the IVP consisting of the DDAE
•1 0 0 0
‚ •x(t)˙
˙ y(t)
‚
=
•0 0 1 0
‚ •x(t) y(t)
‚ +
•0 1 0 0
‚ •x(t−τ) y(t−τ)
‚ +
•1−et10−τ
−t
‚
, (5.36)
for t∈[0,∞),τ=1, with an initial functionφ(t) :=
• t e10t
‚
, for t∈[−τ, 0].
Applying Algorithm5.2to(5.36), we proceed as follows.
Withκ=0we obtain
U =
2 6 6 6 4
0 −1 0
1 0 0
0 0 0
0 0 1
3 7 7 7 5
,M˜= 2 4
0 0
−1 0
−1 0 3 5, N˜=
2 4
−1 0
0 0
0 0
3 5,
Z2 = 2 4
0.7071 0.7071 0.5000 −0.5000
−0.5000 0.5000 3 5,T2=
• 0 1
‚
, Z1=[ ]3,0,
Z2TN˜ =
•0.7071 0 0.7071 0
‚
, Z1TM˜=[ ]0,2. Here by[ ]i,j we denote the empty matrix of the size i by j . Thus,rank
•Z1TM˜ Z2TN˜
‚¶
=1<2, which implies that the shift index is bigger than0 (equiv-alently, the DDAE(5.36)is noncausal).
Withκ=1we obtain
M˜ = 2 6 6 6 4
0 0
−1 0
−1 0
0 0
3 7 7 7 5
,N˜= 2 6 6 6 4
0 −0.7071
0 0
0 0
−1 0 3 7 7 7 5
, ˜P = 2 6 6 6 4
0 0
0 0
0 −1
0 0
3 7 7 7 5 ,
Z2 = 2 6 6 6 4
0.7071 0.7071 0 0.5000 −0.5000 0
−0.5000 0.5000 0
0 0 1.0000
3 7 7 7 5
, T2=[ ]2,0, Z1=[ ]4,0,
Z2TN˜ = 2 4
0 0.5000 0 0.5000 1.0000 0
3
5, Z1TM˜=[ ]0,2,Y1=[ ]0,0, Y2= 2 4
0 −0.7071 0 −0.7071
−1.0000 0 3 5.
In this case rank
•Z1TM˜ Z2TN˜
‚¶
=2 and therefore the shift index isκ=1. The DDAE (5.35) is
•0 0 0 0
‚ •x(t)˙
˙ y(t)
‚
=
•−1 0 0 −0.7071
‚ •x(t) y(t)
‚ +
•0 0 0 0
‚ •x(t−τ) y(t−τ)
‚ +
• t 0.7071e10t
‚ .
Remark 5.21. For linear time invariant DDAEs, the functionsM,P,G,Z2,T2,Z1,Y1, Y2 can be chosen to be constant matrices on the whole interval. Therefore, we only need to compute them at the initial point and can use them in the whole integration process.
Note that Algorithm5.2is successfully applied to the DDAE (1.2a) if and only if the following conditions hold:
i) Assumptions5.1,5.13are satisfied.
ii) The strangeness index is well-defined for the DAE (5.24) so that we can derive the regular, strangeness-free DDAE (5.20).
In more detail, we formulate the following hypothesis to guarantee that Algorithm5.2 can be performed and the regular, strangeness-free DDAE (5.20) can be achieved.
Hypothesis 5.22. Consider the DDAE (1.2a). Suppose that there exist integersκ, a,ˆ mˆ such that for each t∈I, the strangeness indexµ(t)ˆ is well-defined for the DAE(5.24)in a sufficiently small neighborhoodB(t)of t , and the functionsM,P, obtained by building the derivative arrays with k=µˆ(t)as in(5.26), satisfy the following properties:
1. There exists a smooth, full rank matrix function U of size (κ+1)(µ+1)n by mˆ such that its columns span the space corange (M(:, (2n+1) :end)) and UTP(:, (n+1) :end)=0. We set
M˜:=UTM(:, (n+1) : 2n), N˜:=UTM(:, 1 :n).
2. There exists a smooth matrix function Z2of sizem byˆ a onˆ B(t), such that Z2TM˜= 0and Z2TN˜has pointwise full row ranka.ˆ
3. There exists a smooth matrix function T2 of size n by dˆ:= n−a, and Tˆ 2 has pointwise maximal rank such that Z2TN˜T2=0.
4. For all s∈B(t), we have rank( ˜M(s)T2(s))=d so that there exists a smooth, point-ˆ wise maximal rank matrix function Z1 of size mˆ by dˆ satisfying rank(Z1T(s) ˜M(s)T2(s))=d .ˆ
As noted in Section3.3, in the case of advanced DDAEs, the orderµ(resp.,η) ofx(t) (resp.,x(t−τ)) in the strangeness-free formulation (5.14) can be arbitrarily high, and therefore it is complicated to determine system (5.14) pointwise. Furthermore, even if system (5.14) is already known, it will be hard and can be inefficient to implement numerical methods like Runge-Kutta or BDF methods to (5.14). Therefore, another approach is needed to determine whether it is possible to integrate an advanced DDAE.
This is the main purpose of the BVP method presented in the next section.