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.
5.4. Boundary Value Problem method 80 time interval (0,τ)
Ej(t) ˙xj(t)=Aj(t)xj(t)+Bj(t)xj−1(t)+fj(t), j =1, . . . ,`. (5.37) However, unlike Algorithm5.1, we do not solve the sequence of IVPs for the DAEs (5.37) step by step. Instead, we rewrite the IVP (1.2) as the differential-algebraic BVP for the following DAE of size`m×`n
2 6 6 6 4
E1(t)
E2(t) . ..
E`(t) 3 7 7 7 5 2 6 6 6 4
˙ x1(t)
˙ x2(t)
...
˙ x`(t)
3 7 7 7 5
= 2 6 6 6 4
A1(t)
B2(t) A2(t) . .. . ..
B`(t) A`(t) 3 7 7 7 5 2 6 6 6 4
x1(t) x2(t)
... x`(t)
3 7 7 7 5
+ 2 6 6 6 4
B1(t)x0(t)+f1(t) f2(t)
... f`(t)
3 7 7 7 5
, for allt∈(0,τ),
(5.38a)
together with the boundary condition 2
6 6 6 4
In 0 In
. .. ...
0 In 3 7 7 7 5
2 6 6 6 4
x1(0) x2(0)
... x`(0)
3 7 7 7 5
+ 2 6 6 6 4
0
−In 0 . .. ...
−In 0 3 7 7 7 5
2 6 6 6 4
x1(τ) x2(τ)
... x`(τ)
3 7 7 7 5
= 2 6 6 6 4
φ(0) 0
... 0
3 7 7 7 5
. (5.38b)
Let us introduce new notation and rewrite the BVP (5.38) in a compact form as
E(t) ˙X(t)=A(t)X(t)+B(t)x0(t)+F(t), for all t∈(0,τ), (5.39a)
C X(0)+DX(τ)=~v, (5.39b)
withX(t)=£
xT1(t) x2T(t) . . . x`T(t)⁄T
and~v:=£
φT(0) 0 . . . 0⁄T
.
Assuming that the strangeness indexµ=µ(E,A) is well-defined for the DAE (5.39a) and applying Algorithm2.1to (5.39a), one obtains the strangeness-free DAE
2 4
E1(t) 0 0
3
5X˙(t)= 2 4
A1(t) A2(t)
0 3
5X(t)+ 2 6 6 4
Pµ
j=0B1,j(t)x0(j)(t)+F1(t) Pµ
j=0B2,j(t)x0(j)(t)+F2(t) Pµ
j=0B3,j(t)x0(j)(t)+F3(t) 3 7 7 5
, (5.40)
where
•E1(t) A2(t)
‚
has pointwise full row rank.
By studying the solutionX to the BVP (5.39), we derive the solvability of the IVP (1.2) in the following theorem.
Theorem 5.23. Consider the IVP(1.2)on the bounded time intervalI=[t0,t0+`τ). As-sume that the strangeness indexµ=µ(E,A)is well-defined for the DAE(5.39a)with the strangeness-free formulation(5.40). Moreover, let Assumption5.1hold. Then we have:
i) The initial functionφ is consistent if and only ifφ isµ-times continuously
differen-tiable, and the following condition holds:
µ
X
j=0
B3,j(t)φ(j)(t−τ)+F3(t)=0 for all t∈(0,τ). (5.41)
ii) The BVP(5.39)has a unique solution only if the matrix function
•E1(t) A2(t)
‚
, as in(5.40), is pointwise nonsingular.
iii) If in addition E , A∈Cµ+1(I,Cm,n), then there exist pointwise nonsingular matrix-valued functions P∈C(I,C`n,`n), Q∈C1(I,C`n,`n)such that
P
•E1(t) 0
‚ Q=
•I 0 0 0
‚ , P
•A1(t) A2(t)
‚ Q−P
•E1(t) 0
‚ Q˙=
•A11(t) 0
0 I
‚
. (5.42)
Transforming the boundary condition analogously,
£C1 C2⁄
:=CQ(0), £
D1 D2⁄
:=DQ(τ),
then the BVP(5.39)(and consequently the IVP(1.2)) is uniquely solvable if and only if C1+D1is nonsingular.
Proof. i) The claim follows directly by observing thatx0(t)=φ(t−τ) for allt∈[0,τ].
ii) The second claim is obvious, since the number of undetermined components in the functionX is exactlyn−rank
•E1(t) A2(t)
‚¶
.
iii) Observe that the conditionE, A∈Cµ+1(I,Cn,n) implies thatE andA are continu-ously differentiable functions. Thus, there existsP andQ that satisfy (5.42). Partition-ingX =
•X1
X2
‚
correspondingly to (5.42), we see that the algebraic variableX2is fixed by the second block row equation of (5.40). Therefore, the BVP (5.39) is uniquely solv-able if and only if the initial conditionX1(0) of the differential variableX1is specified.
Note that due to the boundary condition (5.39b), one has (C1+D1)X1(0)=~v−C2X2(0)−D2X2(τ).
Therefore,X1(0) is uniquely specified if and only ifC1+D1is nonsingular.
Remark 5.24. Since the functionsP,Qin Theorem5.23cannot be computed numeri-cally, the condition ii) for the existence and uniqueness of the BVP (5.39) must be ver-ified in another way. We refer to [4,129] for one way of doing this by using the funda-mental matrix of the DAE (5.39a).
From the theoretical viewpoint, Theorem 5.23 completely analyzes the solvabil-ity analysis of the IVP (1.2). However, from the numerical viewpoint, it is only feasi-ble to compute the solution x of the IVP (1.2) if the differential-algebraic BVP (5.39) can be efficiently integrated. This naturally leads to the question of investigating the strangeness indexµof the DAE (5.39a), since high indexµwill lower the performance of numerical methods, see [75], Theorems 5.10, 5.12.
In the remaining part of this section, we discuss the strangeness indexµfor different
5.4. Boundary Value Problem method 82 types of the DDAE (1.2a). Certainly, the BVP method is expected to work for retarded or neutral systems. We consider one special case of (1.2a) in the next proposition.
Proposition 5.25. Consider the DDAE(1.2a)and assume that the following conditions hold.
i) The matrix function pair(E,A)is sufficiently smooth with a well-defined strangeness indexµand characteristic quantities dµ, aµ, vµ, uµthat satisfy vµ=uµ=0.
ii) Assumptions5.1,5.13hold.
Then, the DAE(5.39a)has a well-defined strangeness index which is exactlyµ.
Proof. Due to the condition i), we see that for eachj =1, . . . ,`, the DAE
Ej(t) ˙xj(t)=Aj(t)xj(t)+Bj(t)xj−1(t)+fj(t), (5.43) is regular and has a well-defined strangeness indexµ. This means that ifBj(t)xj−1(t)+ fj ∈Cµ+1([0,τ],Cn) then xj is continuously differentiable. This condition is further relaxed by Assumption5.13, which implies that only fj is required to beµ+1 times continuously differentiable, butxj−1only needs to be continuously differentiable.
Therefore, by simple induction, one sees that ifx0∈C1([0,τ],Cn) andfj∈Cµ+1([0,τ],Cn) for all j =1, . . . ,`, there exists uniquely`functionsx1, . . . ,x`∈C1([0,τ],Cn) that solve the DAE sequence (5.43). This implies that the DAE (5.39a) has a well-defined strangeness indexµ.
Remark 5.26. It is worth to note that if one removes the conditionvµ=uµ=0 in condi-tion i) of Proposicondi-tion5.25, then the DAE (5.39a) can have a strangeness indexµbigger thanµ, as demonstrated in the next example.
Example 5.27. Consider the Delay-DAE
•0 1 0 0
‚
| {z }
E
˙ x(t)=
•1 0 0 0
‚
| {z }
A
x(t)+
•0 0 0 1
‚
| {z }
B
x(t−τ)+
•f1(t) f2(t)
‚
(5.44)
on the time interval[0, 2τ]. Algorithm5.2applied to(5.44)results in the strangeness-free DDAE
•0 0 0 0
‚
˙ x=
•0 1 1 0
‚ x(t)+
• f2(t+τ) f1(t)+f˙2(t+τ)
‚ . The underlying DDE of (5.44)therefore is
•0 1 1 0
‚
x(t˙ )= −
• f˙2(t+τ) f˙1(t)+f¨2(t+τ)
‚ ,
is of retarded type, and hence the DDAE(5.44)is of retarded type. One also sees that the strangeness index of the pair(E,A)isµ(E,A)=0and the characteristic quantities are dµ=1, aµ=0, vµ=0, uµ=1.
The BVP approach, however, leads to the matrix function pair(E,A)given by
E = 2 6 6 6 4
0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 0 3 7 7 7 5
, A =
2 6 6 6 4
1 0 0 0 0 0 0 0 0 0 1 0 0 1 0 0 3 7 7 7 5 ,
and the strangeness index of the pair (E,A) isµ(E,A)=1, which is bigger than the strangenessµ(E,A). Therefore the BVP method can increase the strangeness index, even in the case that the considered DDAE is non-advanced.
Now let us turn to the case that the DDAE (1.2a) is of advanced type.
Example 5.28. Consider the following DDAE on the time interval(0,`τ)
•0 1 0 0
‚
| {z }
E
•x(t)˙
˙ y(t)
‚
=
•1 0 0 1
‚
| {z }
A
•x(t) y(t)
‚
−
•0 0 λ 0
‚
| {z }
B
•x(t−τ) y(t−τ)
‚ +
•f1(t) f2(t)
‚
, (5.45)
with a parameter λ∈R. By adding the derivative of the second equation to the first equation, we obtain 0=x(t)−λx(t˙ −τ)+g(t), where g(t) :=f˙2(t)+f1(t). Thus,(5.45) is of retarded type (resp., advanced type) if λ=0(resp., λ6=0). Now we consider the strangeness indexµof the pair(E,A)for different value ofλ.
i) Withλ=0, one sees that
E =diag(E,E, . . . ,E), A =diag(A,A, . . . ,A).
One can directly find that the Kronecker indexνof the pair(E,A)is equal to the Kro-necker index of the pair(E,A). It therefore implies that the two pairs(E,A)and(E,A) have the same strangeness index.
ii) Forλ=1the constraint0=x(t)−λx(t˙ −τ)+g(t)implies that x(t)=x(t˙ −τ)−g(t).
Thus, for any t∈((i−1)τ,iτ)we have x(t)=£
1 0⁄
φ(i)(t−iτ)−
i−1
X
j=0
λjg(j)(t−jτ).
As a consequence,
x`(t)=£ 1 0⁄
φ(`)(t−(`−1)τ)−
`−1
X
j=0
λjg(j)(t−jτ),
which yields that the strangeness index of the matrix function pair(E,A)is`.
Example5.28demonstrates that if the DDAE (1.2a) is of advanced type, the DAE (5.39a) may have an arbitrarily high index. This situation leads to very poor performance or even divergence of numerical methods whenever the desired integration time interval is long. In this case, the numerical error is rapidly amplified in time, as shown in Figure 5.1 below. Here the inhomogeneity f is chosen such that the system possesses the analytical solution
•x(t) y(t)
‚
=
• t et/10
‚
. The simulation is done by using the 3 stage Radau
5.4. Boundary Value Problem method 84 IIA method of [130] with constant step sizeh=0.01 for system (5.45) withλ=1.
0 2 4 6 8 10
10−20 10−10 100 1010
errory
0 2 4 6 8 10
10−10 10−4 102
108 errorx
Figure 5.1:Absolute error vs. length of the time intervalIfor (5.45).
The advantage of the BVP method is that even if the DDAE (1.2a) is of advanced type, the numerical solution to the IVP (1.2) can still be efficiently computed, provided that the strangeness indexµis sufficiently small. This fact is illustrated in the following example.
Example 5.29. Consider the IVP(1.2)for the DDAE
•0 1 0 0
‚
| {z }
E
•x(t)˙
˙ y(t)
‚
=
•1 0 0 1
‚
| {z }
A
•x(t) y(t)
‚
−
•0 0 0 1
‚
| {z }
B
•x(t−τ) y(t−τ)
‚ +
•f1(t) f2(t)
‚
, (5.46)
on the time interval(0,`τ)for some integer`∈N. The inhomogeneity is chosen such that the system possesses the analytical solution
•x(t) y(t)
‚
=
• t et/10
‚ .
Adding the derivative of the second equation of (5.46)to the first equation and elimi-natingx on both sides, one obtains˙
•0 0
‚
=
•x(t) y(t)
‚
−
•0 0 0 1
‚ •x(t−τ) y(t−τ)
‚
−
•0 1 0 0
‚ •x(t˙ −τ)
˙ y(t−τ)
‚ +
•f1(t)+f˙2(t) f2(t)
‚
, (5.47)
and therefore the DDAE(5.46)is of advanced type.
Furthermore, by direct computation, one can verify that the strangeness indexµof the pair(E,A)is exactly1. For illustration, we consider`=10which leads to the numerical solution and absolute error presented in Figure5.2. The simulation is done by using the three stage Radau IIA method of [130] with constant step size h=0.01for system(5.46).
So far in this chapter we have presented two different approaches to analyze the solvability analysis of general linear time varying DDAEs. We compare these two ap-proaches in Table5.3.
0 2 4 6 8 10 10−20
10−15 10−10 10−5
errorx errory
0 2 4 6 8 10
0 5 10
15 x
y
Figure 5.2:Numerical solution and absolute error of (5.46).
The generalized method of steps The BVP method Pros. Does not require any condition for
the time domain, which can be ei-ther bounded or unbounded.
Explicitly shows the necessary and sufficient consistency conditions for bothφandf.
Gives the strangeness-free DDAE pointwise.
Provides the exact smoothness re-quirements for bothφand f.
Explicitly shows the relation be-tweenx(t) andx(t−τ).
Suitable for the investigation of the long time behavior of the solution, for example the stability, controlla-bility.
Cons. Only shows necessary consistency conditions forφandf.
Only applicable for systems on bounded time domains.
In general, it can lead to extra smoothness requirements forφ.
Does not show the relation between x(t) andx(t−τ) explicitly.
Not suitable for the long time be-havior study of the solution.
Table 5.3:Theoretical comparison of the generalized method of steps and the BVP method.
5.4. Boundary Value Problem method 86
Numerical solutions of IVPs for DDAEs
This chapter is devoted to the numerical solution of the IVP (1.2). In view of Chapter5, the theoretical solvability of (1.2) can be completely analyzed via either the generalized method of steps (Algorithm5.1) or the BVP method. In this chapter, these approaches are implemented and examined consecutively in Sections6.1and6.2. Finally, some illustrative numerical experiments are presented.
Within this chapter, we assume that the (integration) time intervalIis bounded, i.e.,
`< ∞, and that the IVP (1.2) has a unique solutionxonI.
6.1 Application of the generalized method of steps for DDAEs
The integration strategy presented in this section is based on the theoretical results of Sections5.1-5.3under Assumption5.13, i.e., the DDAE (1.2a) is of either retarded and neutral type. The main idea is that we do not directly integrate the original DDAE (1.2a) but instead we integrate the strangeness-free formulation (5.35) which takes the form
•Eˆ1(t) 0
‚
˙ x(t)=
•Aˆ1(t) Aˆ2(t)
‚ x(t)+
•Bˆ1(t) Bˆ2(t)
‚
x(t−τ)+
•γˆ1(t) γˆ2(t)
‚
, (6.1)
where
•Eˆ1 Aˆ2
‚
is pointwise invertible.
Recall that in Section5.3, we have presented in Algorithm5.2an efficient approach to compute the strangeness-free DDAE (5.35) pointwise in a stable way. In parallel to the pointwise determination of (5.35), by checking the condition (5.30) we can also ver-ify whether Assumption5.13holds, i.e., the type of the DDAE (1.2a) is not advanced.
Therefore, the remaining work now is to integrate the system (6.1).
Due to the presence of the delay term, it is necessary to calculate a dense output rather than a discrete approximation to the solution for evaluating the solution at the lag points. This requirement is attainable by interpolation or even better, by continuous methods, see e.g. [13]. Due to this reason, the collocation methods seem to be good candidates. Adopted from the solver RADAR5 [60], we use the Radau scheme for the
87
6.1. Application of the generalized method of steps for DDAEs 88 numerical integration, which is given by nodes
0<δ1< · · · <δs=1, s∈N. (6.2) Since the DDAE (1.2a) is assumed to be of either retarded or neutral type and the delay τ is constant, we know all the discontinuity points that may occur in advance. We include all the discontinuity points ofx, ˙x, . . . ,x(s)into the mesh and denote a mesh by π : t0<t1< · · · <tN =tf. (6.3) It should be noted that the numberN in this chapter is the number of mesh points, which should not be confused with the nilpotent matrixNin the Kronecker-Weierstraß canonical form.
With the meshπas in (6.3), the collocation points therefore are
ti j=ti+hiδj, j =1, . . . ,s, (6.4) wherehi is the stepsize used at thei-th step. For the numerical approximation of the solution, we seek for the piecewise polynomial Xπof degrees, i.e., Xπ,i :=Xπ|[ti,ti+1]are polynomials of degrees, which are determined by the following set of equations
•Eˆ1(ti j) 0
‚
X˙π(ti j)=
•Aˆ1(ti j) Aˆ2(ti j)
‚
Xπ(ti j)+
•Bˆ1(ti j) Bˆ2(ti j)
‚
Xπ(ti j−τ)+
•γˆ1(ti j) γˆ2(ti j)
‚
, (6.5)
for alli =1, . . . ,N, j=1, . . . ,s.
If the delay factor
•Bˆ1(ti j) Bˆ2(ti j)
‚
Xπ(ti j−τ) does not occur, then (6.5) represents a linear system for the internal stages Xπ,i(ti j), j =1, . . . ,s for eachi =1, . . . ,`. Otherwise, in the case that
•Bˆ1(ti j) Bˆ2(ti j)
‚
Xπ(ti j−τ) is present in (6.5), then we still have to define the past function Xπ(ti j−τ) which is an approximation tox(ti j−τ). For this we choose
Xπ(ti j−τ)=
(φ(ti j−τ) ifti j−τ60,
Xπ,k(ti j−τ) for some 16k6N that satisfies tk<ti j−τ6tk+1. The continuous output polynomial Xπ,k at thek-th step is given by Lagrange interpo-lation of orders, i.e.,
Xπ,k(tk+θhk)=
s
X
j=0
Lj(θ)Xπ,k(tk+δjhk), (6.6)
whereLj(θ) is the Lagrange polynomial of degrees satisfyingLj(δk)=δk j withδk j
being the Kronecker delta symbol.
Remark 6.1. As noticed in [59,60], one can optionally replace the continuous output polynomial Xπ,kin (6.6) by another dense output polynomial given by
Xπ,k(tk+θhk)=
s
X
j=1
Lj(θ)Xπ,k(tk+δjhk).
The use of onlysinterpolation nodesδj, j=1, . . . ,sinstead ofs+1 nodesδj,j =0, . . . ,s is beneficial in the presence of a jump in the solution at the pointtk, i.e., Xπ,k(tk)6=
Xπ,k−1(tk).
The existence and uniqueness, and the convergence results for the numerical ap-proximation Xπare stated in the following theorem.
Theorem 6.2. Consider the IVP(1.2)and assume that it satisfies Assumptions5.1,5.13.
For N∈Nand s≥1, define the meshπas in(6.3)and for i =0, . . . ,N−1the collocation points ti j, j=1, . . . ,s as in(6.4). Then the following assertions hold.
i) For sufficiently small mesh widths h0, . . . ,hN−1there exists one and only one con-tinuous piecewise polynomialXπthat solves the DAE sequence(6.5)and it is con-sistent at all the mesh point ti.
ii) The convergence order of the collocation method with schemesδj as in(6.2)is s, i.e.,
kXe(t)−Xπ(t)k∞=sup
t∈IkXe(t)−Xπ(t)k =O(hs), whereXeis the exact solution x∈Cs+1(I,Cn)to the IVP(1.2).
Proof. For the proof see Theorem 4, [68] or Theorem 4.2, [60].
We demonstrate the application of the generalized method of steps in the following example.
Example 6.3. Let us revisit the IVP considered in Example5.20, which reads,
•1 0 0 0
‚ •x˙1(t)
˙ x2(t)
‚
=
•0 0 1 0
‚ •x1(t) x2(t)
‚ +
•0 1 0 0
‚ •x1(t−τ) x2(t−τ)
‚ +
•1−et10−τ
−t
‚
, for t∈[0,∞), φ(t)=
• t e10t
‚
, for t∈[−τ, 0].
(6.7) As shown in Example5.20, we see that the shift index of (6.7)isκ=1. Existing solvers, for example RADAR5 [60], fail to handle system(6.7)due to its noncausality. To compute the numerical solution of the IVP(6.7), we implemented the three stage Radau collocation method for the strangeness-free DDAE(6.1), which is pointwise computed automatically by Algorithm5.2. The numerical solution and the absolute error are presented in Figure 6.1. The constant stepsize is h=0.01.