Projected explicit and implicit Taylor series methods for DAEs

The recently developed new algorithm for computing consistent initial values and Taylor coefficients for DAEs using projector-based constrained optimization opens new possibilities to apply Taylor series integration methods. In this paper, we show how corresponding projected explicit and implicit Taylor series methods can be adapted to DAEs of arbitrary index. Owing to our formulation as a projected optimization problem constrained by the derivative array, no explicit description of the inherent dynamics is necessary, and various Taylor integration schemes can be defined in a general framework. In particular, we address higher-order Padé methods that stand out due to their stability. We further discuss several aspects of our prototype implemented in Python using Automatic Differentiation. The methods have been successfully tested on examples arising from multibody systems simulation and a higher-index DAE benchmark arising from servo-constraint problems.


Introduction
Higher-index DAEs do not only represent integration problems but also differentiation problems, as well (see, e.g., [4,22,23,28]). Therefore, it seems worthwhile to solve an associated ODE with classical integration schemes and the differentiation problems using Automatic Differentiation (AD). However, depending on the structure, both differentiations and integrations may be intertwined in a complex manner such that this plausible idea may be difficult to realize in general.
In this context, different approaches have been considered for DAEs in order to combine AD with ODE integrations schemes. By construction, the approaches are based on corresponding index definitions and lead, therefore, to quite different algorithms.
• In [26,27] and the related work, the structural index was used to determine the degree of freedom. • In [13], we used the tractabiliy matrix sequence to solve the inherent ODE for DAEs of index up to two. The generalization to higher-index DAEs seemed rather complicated. • In [17], we described briefly how an approach based on the differentiation index defined in [14,16] leads to an explicit Taylor series methods for DAEs. An analysis of the corresponding projected explicit ODE can be found in [18]. These methods can be considered as projected Taylor series methods.
In this paper, we analyze more general classes of the latter mentioned projected Taylor series methods. In particular, we discuss how projected implicit Taylor series methods can be defined for DAEs, generalizing the approach from [17]. Here we focus on the methods from [7,21].
The main idea in this context is that the computation of Taylor coefficients of a solution of an implicit ODE can be considered as the solution of a nonlinear system of equations. In this sense, we will see that a generalization for DAEs can be obtained by solving a nonlinear optimization problem [16]. The obtained solution corresponds to a projected method. There are several advantages of this approach: • We assume weak structural properties of the DAEs (1), such that ODEs and semiexplicit DAEs are simple special cases. Theoretically, we can consider DAEs of any index. • An explicit description of the inherent dynamics is not required for the algorithmic realization. Indeed, due to our formulation as an optimization problem, implicitly we consider the orthogonally projected explicit ODE introduced in [18]. • We can use higher-order integration schemes, also for stiff ODEs/DAEs.
The described methods were implemented in a prototype and first numerical tests for DAEs up to index 5 and integration schemes up to order 8 were successful.
The purpose of our implementation is the improvement of our diagnosis software [17]. Therefore, we do not focus on high-performance, but on information about the reliability of the numerical result, specially with regard to higher-index issues and the diagnosis of singularities like [19]. Due to their stability and order properties, the new higher-order methods presented here are a considerable improvement in comparison to [17], where only explicit Taylor series methods were considered.
The paper is organized as follows. In Section 2, we introduce DAEs and summarize some of our previous results that are crucial for the approach presented here.
The properties of Taylor coefficients in the derivative array of a DAE are discussed in Section 3. Using these properties, in Section 4, we define the general projected explicit/implicit Taylor series method that are, indeed, a generalization of the explicit method from [17].
Within this framework, in Section 5, we present four different types of projected methods: explicit methods, fully implicit method, two-halfstep (TH) schemes, and higher-order Padé (HOP) schemes.
The properties of the considered optimization problems that provide the projection are discussed in Section 6 and some practical considerations for the implementation are addressed in Section 7.
Our prototype implementation of the proposed projected methods for DAEs is tested on several well-known examples and benchmarks from the literature in Section 8. An outlook discussing directions for further investigations concludes this paper.
For completeness, in the Appendix, we summarize the stability functions and stability regions for the considered Taylor series methods, since they are essential for the development of HOP methods. We also summarize some linear algebra results for decoupling DAEs and provide the DAE formulation of the tested examples resulting from servo-constraint problems for multi-body systems.

DAEs: index, consistent values, and decoupling
In this article, we consider general DAEs: where the partial Jacobian f x is singular and ker f x (x , x, t) is constant. For our purposes, we define the constant orthogonal projector Q onto ker f x as well as the complementary orthogonal projector P := I − Q.
Recall that the singularity of f x means that (1) contains derivative-free equations, called explicit constraints, and that the differentiation of (1) may lead to further derivative-free equations, called hidden constraints. A consistent initial value x 0 has to fulfill all explicit and implicit constraints. The characterization of all these constraints motivated the following definition for the differentiation index, cf. [4].
Definition 1 [14] The differentiation index is the smallest integer μ such that: . . .
uniquely determines Qx as a function of (P x, t).
If μ is the differentiation index according to Definition 1, then the conventional differentiation index (see, e.g., [4]) results to be μ as well. According to this definition, in the following we will never prescribe initial values for Qx 0 , since we may compute Qx 0 evaluating a function at (P x 0 , t 0 ). Moreover, in the higher-index case, the Eqs. (2)-(4) contain explicit and hidden constraints that restrict the choice for P x 0 .
According to [16], for an initial guess α ∈ R n , consistent initial values x 0 can be computed solving the following constrained optimization problem: subject to all explicit and hidden constraints.
Equivalently, we can solve the system of equations: all explicit and hidden constraints (8) for a suitable orthogonal projector with rank = d ≤ rank P , where d is the degree of freedom of (1), cf. [15]. However, in particular for nonlinear DAEs, may be difficult to compute and therefore, in practice, (5) may be more convenient than (7), cf. Appendix 2. The approach (5)-(6) was implemented in InitDAE, a Python program to determine the index and consistent initial values for DAEs [12,17].
Furthermore, for linear DAEs of the form: in [18], it was shown that using the derivative array (2)-(4), the decoupling: can be obtained for suitable functions ϕ 1 , ϕ 2 , where (10) is an ODE in the invariant subspace im . Therefore, theoretically, we can set up (10) and solve it with an integration scheme for ODEs. Subsequently, (I − )x can be computed at each time point using (11). Notice that in doing so, the error in (I − )x depends only on the error made solving (10) and the properties of ϕ 2 from (11). Moreover, the values of (I − )x at previous time points do not influence (10). For nonlinear DAEs, analogous considerations can be undertaken considering the linearization along a solution. However, of course, the properties of f are decisive in practice. Since a detailed analysis for the nonlinear case goes far beyond the scope of this article, we focus on a general formulation of projected Taylor methods using (5)- (6), having in mind that at present, the theoretical basis has been developed for linear DAEs only. At least the numerical tests from Section 8 suggest the applicability for some classes of nonlinear DAEs.

Taylor series and DAEs
Since we wish to analyze one-step methods, we consider the computation of an approximation of the solution x(t) of the ODE/DAE (1) at time t j +1 , given an approximation of the solution at time-point t j . Consequently, in order to describe our method in terms of Taylor expansion coefficients, for k c ∈ N, we suppose that a suitable approximation: is given and that we look at adequate methods to compute: If we suppose that the ODE/DAE is described by (1), we require that: are fulfilled. For our purposes, given K ≥ 1, we further consider the order (K − 1) derivative array derivative array [4] containing (1) and K − 1 derivatives of (1): Therefore, we suppose that for the corresponding function: . . .
it holds: and In practice, the function r can be provided by automatic differentiation (AD) [17,31]. Using this notation, the index from Definition 1 is the smallest integer μ such that for K = μ, the derivative array r uniquely determines Qc 0 as a function of (P c 0 , t). For this purpose, the 1-fullness of the Jacobian of r is verified in [16,17], cf. Appendix 2. We emphasize that the main difference to the conventional differentiation index [4] is precisely that for 1-fullness the columns corresponding to c 0 (and not c 1 ) are considered. With this index definition in mind, we can define consistency for the Taylor coefficients.

Definition 2
For K ≥ μ and 0 ≤ k c ≤ K − μ, the Taylor coefficients up to k c are consistent if they are in the set: Note that T 0 0 corresponds to the set of consistent initial values and that if sufficient smoothness of f is given, we can suppose that for all c 0 ∈ T 0 0 in regularity regions, there is a unique solution of the initial value problem. For a discussion of regularity regions and singularities within a projector-based analysis, we refer to [14,23].
For sufficiently smooth regular linear DAEs (9), Theorem 1 from [18] implies that for any c 0 ∈ T 0 0 , there is a unique solution fulfilling Indeed, Theorem 1 from [18] provides a general description of the inherent dynamics in terms of the associated orthogonally projected explicit ODE (10).
Let us focus on the relation between k c , K, the DAE-index μ and the computation of consistent initial values and Taylor coefficients at a particular t 0 considering: • For initial value problems x(t 0 ) = c 0 for ODEs g(x (t), x(t), t) = 0 with regular g x , if we consider (17), then we can compute K consistent coefficients: at t 0 , since c 0 is given. In the above notation, the maximal value for k c is k c = K. • If we consider an uniquely solvable nonlinear time-dependent equation g(x(t), t) = 0 with regular g x and a corresponding system of equations (17), then, at t 0 , we cannot prescribe c 0 and compute therefore K consistent coefficients: For the coefficients (c K ) 0 , no equations are given, since in this particular case, they do not appear in (17). Note that in principle, g(x, t) = 0 can be considered an index-1 DAE. In this sense, it fits into the case below. • For DAEs (1), if we consider (17) and fix the free initial conditions of c 0 , then in general, we may compute K + 1 − μ consistent coefficients: cf. [17]. In this case, we have at most k c = K − μ. In general, the coefficients c K+1−μ , . . . , c K cannot be computed considering (17). Another crucial aspect is that not all components of c 0 can be prescribed, since all the constraints have to be satisfied.
Note that according to (5)-(6), for an arbitrary initial guess α that, in general, may be not consistent, the optimization problem: provides consistent initial values (18). Moreover, in terms of (7)- (8), this minimization problem is equivalent to the system of equations: where describes an appropriate orthogonal projector, and the rank of coincides with the number of degrees of freedom of the DAE [15,18]. Note further that, in general, the coefficients [(c K+1−μ ) 0 , . . . , (c K ) 0 ] are not uniquely determined neither by (19)- (20) nor by (21)- (22). In our implementation from [17], the minimum norm with the general solution: For this clearly structured example, the projector-based approach will lead to: This means that Qx corresponds to x 2 and according to the notation introduced in [18], the EOPE-ODE (the essential orthogonally projected ODE describing the dynamics) that corresponds to (10) will be formulated in terms of x 1 . Let us have a closer look to the derivative array with respect to the index determination and the computation of consistent initial values, both related to the computation of , see Appendix 2.
According to (12), for n = 5 and K = 4, we consider: such that the equations on the left, that are formulated for functions as described in (13), correspond to the equations on the right, that can be formulated at some t = t j for the scalar numbers (14)): The colors visualize the chains of calculation by with each entry of (I − )c 0 is uniquely determined by the equations r = 0, cf. (15). In particular, the red equations permit the computation of Qc 0 , since: Since no representation of c 20 is possible with less differentiations, the index is μ = 4. Moreover, the violet, blue, and green expressions provide values for components of P c 0 , in particular c 50 , c 40 , and c 30 , respectively. This means that we cannot prescribe initial values for (P − )c 0 .
Indeed, the EOPE-ODE reads: Summarizing, we see that to compute (c 0 ) 0 , we have to prescribe a value for x 1 (t 0 ) and consider at least derivatives of order up to three (with K = μ = 4) in the derivative array, i.e., r( We emphasize that the gray items must vanish in order to satisfy r = 0, but for K = μ − 1, they do not uniquely determine all coefficients of (c ) 0 for any > 0. This means that the value for k c from Definition 2 is k c = 0 = K − μ.
Since with the approach (19)-(20) the projector is not computed explicitly and at least we consider nonlinear under-determined systems of equations, they are solved in a minimum-norm sense. Therefore, the used solver obtains values for all higher derivatives, whereas we cut off (c K−μ+1 ) 0 , . . . , (c K ) 0 , since only (c 0 ) 0 , . . . , (c K−μ ) 0 are consistent in the sense of Definition 2.
In Table 1, we present the results of the computation of consistent initial values with InitDAE [12] that solves (19)- (20). For the considered initial value We can appreciate that for K = 5, only the Taylor coefficient (c 0 ) 0 and (c 1 ) 0 are consistent, i.e., k c = 1 = K − μ. Increasing K, correspondingly more consistent Taylor coefficients could be computed.
The numerical solution delivered by the methods defined in the following corresponds to: • the numerical solution obtained by Taylor series methods applied to the projected explicit ODE (10) for x, and • corresponding values for the components (I − )x that result from (11).
Therefore, the stability and order properties of the integration methods defined below can be transferred from ODEs to DAEs. Due to the formulation as an optimization problem, the inherent dynamics of the DAE that can be expressed in terms of x is not considered explicitly, but implicitly. Table 1 Numerical solution of the initialization problem for system (23)

Recall that:
• P = = I holds for ODEs and that, therefore, for ODEs, the approach (19)- (20) means to compute the Taylor coefficients if c 0 is prescribed. • We assumed that ker f x (c 1 , c 0 , t) and therefore also P do not depend on (c 1 , c 0 , t). Therefore, the Taylor coefficients of P x(t) at t j correspond to: With these two properties in mind, we can present a very general formulation for implicit and explicit Taylor series methods for ODEs and DAEs defining suitable objective functions instead of (19).
In a first step, we focus on consistency.

Lemma 1 Consider:
Then for any solution: of the minimization problem: are consistent at t j +1 for all k i ≤ K − μ .
Proof According to Definition 2, the coefficients (29) are consistent, since the constraints (28) are fulfilled. Recall further that, under suitable assumptions, the solvability of (27)-(28) follows from the Definition 1 of the index μ and has been discussed in [16].  (28) for the ODE (10) that is invariant in the subspace im , and the subsequent computation of the remaining components according to (11).
Proof On the one hand, according to (21)- (22), the solution of (27)-(28) for the original DAE delivers the same result (29) as: On the other hand, the derivative array (31) for the DAE contains the derivative array of the ODE (10) and the nonlinear equations (11).
With the notation of the objective function, different Taylor integration methods can be described with corresponding weights ω e e , ω i i . This provides us a very flexible way to implement schemes with different order and stability properties.

Explicit Taylor series method for DAEs
In terms of the above notation, the explicit Taylor series method for ODEs corresponds to k e ≥ 1, ω e e = 1 for 0 ≤ e ≤ k e , k i = 0, ω i 0 = 1. Recall that the approach from [17] for DAEs consists of the following steps: • Initialization: Solve the optimization problem: for an initial guess α. • For time-points t j +1 , j ≥ 0, h j = t j +1 − t j : Solve the optimization problems: successively for k e ≤ K − μ, where k e is the order of the series in t j .
This method is called explicit, since (32) is an explicit equation for (c 0 ) j +1 that does not involve any value (c ) j +1 for ≥ 1. In contrast to explicit ODEs, where Taylor coefficients may be obtained by function evaluation, with this approach for DAEs in general, we have to solve a nonlinear system of equations. Therefore, it seems reasonable to consider also implicit Taylor approximations in the integration scheme.

Fully implicit Taylor series methods for DAEs
The implicit counterpart of the explicit Taylor series method for ODEs corresponds to k i ≥ 1, ω i i = 1 for 0 ≤ i ≤ k i , k e = 0, ω e 0 = 1. Our generalization for DAEs consists, therefore, of the following steps.
• Initialization as in Section 5.1, Eqs. (32)-(33). • For time-points t j +1 , j ≥ 0, h j = t j +1 − t j : Solve the optimization problems: Obviously, if, instead of (34) and (36), more general conditions of the type (27) are considered, then the dimension of the system that has to be solved remains equal. Therefore, it seems natural to consider more general schemes with better convergence and stability properties than the explicit and the fully implicit Taylor series methods.

Two-halfstep explicit/implicit schemes
One straightforward combination of the explicit and implicit integration schemes is to approximate x(t j + σ h j ) = x(t j +1 − (1 − σ )h j ) for 0 ≤ σ ≤ 1 as follows: and equalize the expressions from both right-hand sides. The properties of the methods (38)-(39) are described in [21]. The choice σ = 1 2 , which can be interpreted as a generalization of the trapezoidal rule, turns out to be convenient. For k e = k i , σ = 1 2 , it coincides with the one tested in [2,9].
Remark 1 Note that another closely related implicit/explicit scheme is described in the literature, see, e.g., [29]. There, the first step is implicit and the second one explicit, in contrast to the approach from above. According to the extensive analysis from [29], σ = 1 2 is convenient also in that case. However, for 0 < σ < 1, these methods are less suitable for our DAE-scheme since the Taylor coefficients would be considered at t j + σ h j , whereas the constraints have to be fulfilled at t j and t j +1 .
In the notation from Section 4, choosing σ = 1 2 in (38)-(39) means to consider: (c e ) j For shortness, we denote these two-halfstep methods by (k e , k i )-TH.
Recall that, in general, the stability function R(z) (cf. Appendix 1) of a (k e , k i )-TH method is not a Padé approximation of the exponential function. Consequently, the maximally achievable order of the integration method for fixed k e and k i is not given for these particular weights in general. Therefore, further higher-order schemes for stiff ODEs have been developed, namely the HOP-methods described in [7].

Higher-order Padé methods
According to [7], HOP may be interpreted as Hermite-Obrechkoff-Padé or higherorder Padé. The corresponding integration schemes may be considered as implicit Taylor series methods based on Hermite quadratures.
In our notation, a (k e , k i )-HOP scheme means choosing: These coefficients correspond to the (k e , k i )-Padé approximation of the exponential function such that the stability function R(z) is precisely this approximation, see Appendix 1. Indeed, (k e , k i )-HOP methods have the following properties, cf. [7]: • the order of consistency is k e + k i , • the order of the local error is k e + k i + 1, Note that also in this case, the trapezoidal rule corresponds to k e = k i = 1 and the implicit Euler method to k e = 0, k i = 1. In this sense, the methods with k e = k i could be viewed as a generalization of the trapezoidal rule and those with k e = k i − 1 as a generalization of the implicit Euler method, cf. [7]. In Section 8, we numerically verify the outstanding properties of these methods.

Properties of the minimization problems
In [16], we analyzed the properties of the minimization problem obtained when computing consistent initial values. That analysis can directly be applied to the explicit Taylor series method, cf. [17]. To appreciate the properties for implicit methods (i.e., k i > 0), we define, for k ≥ max {k e , k i }, the matrices: T e := P ω e 0 P ω e 1 h j P ω e 2 h 2 j . . . P ω e k h k j ∈ R n×n·(k+1) assuming ω i = 0 for l i > k i and ω e = 0 for e > k e , and the vectors: With this notation, we write: ω e e (c e ) j h j Therefore, as in [16], for α := T e X j , X := X j +1 , we consider the objective function: Observe that the matrix: is, by construction, positive semi-definite. However, it is not an orthogonal projector in general. Therefore, Theorem 1 and Corollary 1 of [16] cannot be applied directly. Hence, the solvability of the optimization problem is more difficult than for explicit Taylor methods. More precisely, we want to emphasize that, for: the nullspaces: ker P i G T G 0 and ker P G T G 0 may be different. However, since P i depends on h j , it is reasonable to assume that a suitable stepsize h j can be found such that the optimization problem becomes solvable in the sense discussed in [16].

Dimension of the nonlinear systems solved in each step
For a given K ∈ N, the Lagrange approach for solving (27)-(28) leads to a nonlinear system of equations of dimension 2n · (K + 1), cf. [16]. Thereby, consistent coefficients: . . , (c K−μ ) j +1 are obtained. In contrast, the coefficients c K−μ+1 , . . . , c K will not be consistent in general and the introduced Lagrange-multipliers are not even of interest.
However, increasing K by one means solving a nonlinear system containing 2n additional variables and equations.

Setting k e and k i in a simple implementation
Dealing with automatic differentiation (AD), the number K ∈ N has to be prescribed a priori in order to consider (K + 1) Taylor coefficients. Since 0 ≤ k e , k i ≤ K − μ must be given in general, for the (k e , k i ) TH and HOP methods, we set: k i := K − μ and k e := k i , by default. We further tested k i := K − μ, k e := k i − 1 for HOP methods. So far, we considered schemes with constant order and stepsize only.
• The Jacobian of the constraints (28) is described in [17], since it is also used for the computation of consistent initial values. • To describe the Jacobian of the objective function (27), which is a gradient, we define: and realize that:

Order validation
To visualize the order of the methods, we integrate Example 1 in the interval [0, 1] with different stepsizes. The results can be found in Fig. 1. On the left-hand side, we show the results for the index-4 DAE. On the right-hand side, we report the results obtained for the corresponding ODE described in (25). Summarizing, we observe that: • For k i , k e ≤ 1, the methods coincide with the explicit and implicit Euler methods or the trapezoidal rule. Therefore, the graphs coincide up to effects resulting from rounding errors. • The similarity of the overall behavior for the DAE and the ODE is remarkable. • As expected, the HOP methods are considerably more accurate due to the higher order. • For small h and large k e , k i , scaling and rounding errors impede more accurate results in dependence of the tolerance ftol from the module minimize from SciPy.

Validation of known results
We report numerical results obtained by the methods (3,3)-HOP and (4,4)-HOP for the following examples from the literature: • Mass-on-car from [30], see Appendix, Section C.1, • Extended mass-on-car from [25], see Appendix, Section C.2, • Pendulum index 3, which can be found in almost all introductions to DAEs, in the reduced to first-order form, with the positive y axis pointing upwards and the parameters m = 1.0, l = 1.0, and g = 1.0. (x 1 , x 2 ) are the coordinates, (x 3 , x 4 ) the corresponding velocities and x 5 the Lagrange parameter. In our computation, the system starts from rest at 45 degrees to the vertical. • Car axis index 3 formulation with all parameters as given in [24]. In order to avoid a disadvantageous scaling of the Taylor coefficients, we changed the independent variable t to τ = 10 t. This is advantageous, since the time-dependent input-function is y b = 0.1 sin(10 t). For large K, the corresponding higher Taylor coefficient lead to considerable scaling differences that are avoided by the substitution with τ . For the details of our reformulation, we refer to [17].
For all examples (see also Table 2), we used ftol for the tolerance of the module minimize from SciPy. To estimate the error, we considered the difference between the results obtained by (3,3)-HOP and (4,4)-HOP.
All tests confirmed the applicability of the method. The solution graphs look identical with those given in the literature [24,25,30]. The graphs of the estimated errors of the (3,3)-HOP methods in Fig. 2 confirm the order expectations.
Since it is obvious that our implementation is not competitive with respect to runtime (see Table 3), we have not made a systematic comparison with other solvers here.

A challenging index 5 DAE
We consider now the index-5 example from [26] resulting from a model of two pendula, where the Lagrangian multiplier λ 1 of the first one controls the length of the second one: Note that in this formulation from [26] the positive y axis is pointing downwards. The DAE has index 5 and four degrees of freedom. For the numerical tests, we use the gravity constant g = 1, the length of the first pendulum L = 1, the parameter c = 0.1 and the interval [0, 80]. In [26], this example was integrated with constant stepsize h = 0.05 and order k = 7 as well as h = 0.025 and order k = 8. For the component x 2 , the two solutions were very close until about t = 30 and clearly diverging from there to about t = 50 and totally unrelated from t = 55 on.
Our implementation leads to a very good results in the sense that two solutions have a small difference up to much larger t. We compare the solutions for K = 9 (with HOP method k i = k e = 4 and order 8) and K = 8 (with HOP method k i = k e = 3 and order 6) for the (numerically) consistent initial value:   Table 3  Table 3 Overview of the computations carried out for Fig. 2  The quality of our results is visualized for x 2 in Fig. 3. For slightly perturbed initial values for x 2 and y 2 , the obtained solutions behave analogously. For arbitrarily perturbed initial values, convergence difficulties during the computation of consistent initial values may appear.

Andrews squeezing mechanism
Finally, we want to report here the behavior we obtained for a well-known index-3 benchmark problem with an extreme scaling. According to [20,24], the problem is of the form: ⎛ for q ∈ R 7 , λ ∈ R 6 . To our surprise, using for α the initial value given in the literature, which leads there to an dynamic behavior, our computation of consistent initial values delivers a stationary solution, such that all our integration methods provide the same constant solution for: and all other components are (numerically) equal to zero. Therefore, we explain here why this happened.
First of all, we want to mention that the indicated condition number, introduced in [14], corresponding to the DAE at α is about 10 11 , such that a clear hint to the scaling difficulties is given. In contrast, at the given stationary solution, the condition number is about 10 6 .
To simplify further considerations, we notice that the last four equations of g(p) = 0 are used to compute , δ, , , such that we can neglect them and consider only g 1,2 and f 1,2,3 to determine β, , γ, λ 1,2 . For a constant solution v = w = 0, the relevant equations are therefore: At the equilibrium point corresponding to (41)-(42), the constant drive torque mom, the spring force f 3 , and the Lagrangian forces are equalized. Therefore, it only remains to explain why the approach (19)-(20) delivers a stationary solution. Due to the high condition, the Lagrange multipliers λ are numerically difficult to compute. In fact, for other numerical computations, the accuracy for λ is not controlled [20, p. 536ff]. In contrast, if no (numerical) full row rank of the JacobianG is given, then (19)-(20) computes a minimal norm solution [15], that in this case minimizes the values for λ, leading to the stationary solution. To our knowledge, this stationary solution was not reported before in the DAE literature.
We plan further investigations on this unexpected behavior. Indeed, for some perturbed initial values, we obtained a solution that converge towards the stationary solution (41)-(42). Moreover, with scaled equations and very different initial values, a nonconstant solution that behaves like the one described in [20,24] has been obtained.

Summary and future work
In this article, we presented a projection approach that permits the extension of explicit/implicit Taylor integrations schemes from ODEs to DAEs. As a result, we obtained higher-order methods that can directly be applied also to higher-index DAEs. The methods are relatively easy to implement using InitDAE and convenient since, thanks to the formulation as an optimization problem, the inherent dynamics of the DAE are considered indirectly. We analyzed in detail explicit, fully implicit, two-halfstep (TH) and higher-order-Padé (HOP) methods. Particularly HOP methods present excellent stability and order properties.
The results obtained by a prototype in Python that is based on InitDAE [12] outperform our expectations, in particular with regard to the accuracy for higher-index DAEs, cf. Section 8.2.2. Until now, our focus was on the extension from ODEs to DAEs in order to use higher-order and A-stable methods with InitDAE for our diagnosis purposes during the integration to monitor singularities [17]. With this promising first results, we think that further developments of these projected methods are worthwhile.
In fact, at present, our implementation is not competitive by far. One reason is that setting up the nonlinear equations (27)- (28) and the corresponding Jacobians with AlgoPy, cf. [31], is still very costly. If equations (27)- (28) and the corresponding Jacobians are supplied in a more efficient way, competitive solvers might be achieved. At present, we do not even consider the sparsity of matrices. Furthermore, an improvement seems likely if we take advantage of specific structural properties, e.g., solving subsystems step-by-step, cf. [10,11]. Another reason for our high computational costs is that the package minimize from SciPy often performs more iterations than we expected (often more than 30), although a good initial guess computed with the explicit Taylor series method is given in general. This behavior has to be inspected in more detail. For linear systems, a direct implementation considering the projector from (21) (or, more precisely, a corresponding basis) should deliver an efficient algorithm. This could be of interest, e.g., for the applications from [25,30]. Last but not least, competitive solvers require adaptive order and stepsize strategies-a broad field for future work.
Although these algorithms open new possibilities to integrate higher-index DAEs, we want to emphasize that, in practice, a high index is often due to modelling assumptions that should be considered very carefully. The dependencies on higher derivatives should always be well founded.
Funding Open Access funding enabled and organized by Projekt DEAL

A.1 Stability functions
The general definition (27)-(28) allows for a straightforward description of the stability function. Applied to ODEs (and therefore P = I ), the stability function R : C → C results if we consider the test-ODE: and describe the numerical method for constant h = h j in terms of: For ODEs, the methods described in Section 4 imply: ω e e (c e ) j h j e and, for the test-equation (43), we obtain from: i.e., for z = h j λ ∈ C :

A.2 Stability regions
The corresponding stability regions can thus be characterized by: For the methods discussed in this article, we obtain: • Explicit Taylor: z e e ! .

Fig. 5
Colored stability regions S of two-halfstep schemes considering R ke,ki for all combinations of k e , k i = 0, . . . , 6, where k i corresponds to the rows and k e to the columns. The symmetry with respect to the main diagonal of the images results from (44). We can realize that, for k i = k e = 5, 6, they are not A-stable. This contradicts the statement in [21] Moreover, analogously as before, we have: .
Therefore, in Fig. 6, we obtain again symmetric stability regions that are in accordance with the stability properties reported in Section 5.4.
Since we obtained (44) and (46) for TH and HOP methods with k = k e = k i , it holds for these methods that: 1 = R k,k (it)R k,k (−it) = R k,k (it)R k,k (it) = R k,k (it) for all t ∈ R. Colored stability regions S of HOP-methods considering R ke,ki for all combinations of k e , k i = 0, . . . , 6, where k i corresponds to the rows and k e to the columns. The symmetry with respect to the main diagonal of the images results from (46). For k e = k i , we observe S = C − is 1-full with respect to s 1 , if s 1 is uniquely determined for any consistent b, cf. also [5]. Since we focus on properties of the matrix A , according to [16], we prefer the following equivalent formulation: For our analysis of DAEs, we focus on 1-full matrices A of the following particular form: for an orthogonal projector P . For DAEs, G 1 , G 2 result from the derivative array and P is the orthogonal projector introduced in Section 2. Let W 2 denote an arbitrary matrix fulfilling ker W 2 = im G 2 . Then for N := W 2 G 1, the 1-fullness of (48) implies: For linear DAEs, this projector allows for the formulation of the orthogonally projected explicit ODE (10) in the invariant subspace im , cf. Theorem 1 in [18]. In practice, we can avoid the explicit computation of the orthogonal projector considering the constrained optimization problem (5)-(6) instead of (7)-(8).

C.1 Example: mass-on-car
The DAE resulting from the spring-mass system mounted on a car from [30] corresponds to: y d (t) = y 0 + p 9 t t max · (y f − y 0 ) for 0 ≤ t ≤ t max y f for t ≥ t max for y 0 = 0.5, y f = 2.5, t max = 6.0 and the polynomial p 9 described in Appendix C.3. For 0 < α < π 2 , the DAE-index is 3 and the projector from (21) reads: i.e., it depends on α only and is independent of the other parameters.