Parareal for Higher Index Differential Algebraic Equations

This article proposes modifications of the Parareal algorithm for its application to higher index differential algebraic equations (DAEs). It is based on the idea of applying the algorithm to only the differential components of the equation and the computation of corresponding consistent initial conditions later on. For differential algebraic equations with a special structure as e.g. given in flux-charge modified nodal analysis, it is shown that the usage of the implicit Euler method as a time integrator suffices for the Parareal algorithm to converge. Both versions of the Parareal method are applied to numerical examples of nonlinear index 2 differential algebraic equations.


Introduction
The time-domain simulation of models from physics, finance or social sciences often leads to differential and algebraic equations systems. The systems of equations describing the transient behaviour of the required physical quantities can be both time dependent differential equations, such as e.g. in circuit models for the simulation of microchips and energy networks [28], or space and time dependent partial differential equations (PDEs), as is the case for the simulation of the electromagnetic behaviour of electric machines [25]. In the latter one, typically the method of lines is used, where first a spatial discretisation method is applied to the PDE to obtain an only time dependent system of differential equations which is then solved in time as an initial value problem (IVP).
The time domain simulation of large systems of equations e.g. obtained from fine meshes as well as fast dynamics of the excitations which require small time step sizes together with large time windows considerably increase simulation time. In these cases, parallelisation methods allow reducing computation time. When spatial parallelisation by means of domain decomposition methods is used up to saturation or whenever the time domain dynamics is the bottleneck of the simulation time, parallel-in-time methods [24,21,12,29] can be used. Parareal is such an algorithm which is based on the same idea as multiple shooting methods [21,12].
Many initial value problems arising from physical systems such as e.g. electric networks, spatial discretisation of some approximations to Maxwell's equations or constrained mechanical systems such as the pendulum are systems of differential algebraic equations (DAEs). These are systems that contain both ordinary differential equations (ODEs) as well as algebraic constraints. They convey analytical and numerical difficulties that do not arise when handling ODEs. This includes their potential large sensitivity towards small high frequent perturbations [3,20] as well as their non-trivial selection of appropriate initial conditions [20,7]. One way of classifying them according to the difficulties they pose is given by their index, a natural number ≥ 0. Especially when using less standard algorithms, such as e.g. Parareal, on DAEs with higher index, a correct handling of the equations is of utmost importance. This work focuses on the application of the Parareal algorithm to index 2 DAEs. Lowindex Problems have already been numerically solved by Parareal, e.g. [4]. The rather straight-forward index-1 case was already discussed in [27,11,28].
The paper is structured as follows: Section 2 introduces the concept of differential algebraic equations, the tractability index and the problems that arise for the choice of consistent initial conditions. A first result is given for the behaviour of the implicit Euler scheme on DAEs with a specific structure. In Section 3 the classic Parareal algorithm is presented and a modification of the algorithm is proposed for its application to index 2 tractable DAEs. The (possibly) index 2 DAE for circuit simulation arising from the modified nodal analysis is presented and Parareal is applied to two nonlinear index 2 DAEs in Section 4. The paper concludes in Section 5 with a summary.

Differential Algebraic Equations
We consider initial value problems consisting of quasilinear differential algebraic equations of the form ⊂ R is a time interval and n dof the number of degrees of freedom and initial condition x(t 0 ) = x 0 . Here, det A(x, t) can be zero, which yields a system of equations containing both differential equations as well as algebraic constraints. Differential algebraic equations are typically classified according to their index [3,17,20], which allows evaluating the analytical and numerical difficulties the system may convey. Higher index systems (index ≥ 2) require a special numerical handling. There are different types of index definitions, that essentially coincide for linear systems [23]. In this paper we introduce the projector-based tractability index [20] as it allows a separation of the degrees of freedom and equations into the purely differential and the algebraic ones.
Analogously to [7,5], we consider DAEs with a special structure and therefore take the following assumption.
Assumption 1 (Mass matrix) We assume the spaces ker A(x, t) and Im A(x, t) are independent of the degrees of freedom x and depend smoothly on t.
Remark 1 This assumption is mild, as many DAEs arising from physical systems fulfil those requirements. Later on, it is shown that e.g. RLC circuits described with modified nodal analysis have this property (see [9,7]) or even space discretised partial differential algebraic equations such as e.g. the eddy current problem [26]. Similar assumptions are taken e.g. in [19] and ensure the BDF method integrates index 2 problems well.
Let us consider quasilinear DAEs fulfilling Assumption 1 and their corresponding projectors Q(t) onto ker A(x, t) as well as its complementary P(t) = I − Q(t). We introduce the matrices , and the projectors Q 1 (y, x, t) onto ker A 1 (y, x, t) and P 1 (y, x, t) = I−Q 1 (y, x, t). Finally, the matrix is defined. As we only focus on index 2 systems, we present the tractability index definition accordingly. However, it can be generalised to systems with index > 2 (see [20]).
Definition 1 (Tractability index [20]) A quasilinear DAE (1) fulfilling Assumption 1 has tractability index For the rest of the paper, we will choose a specific projector Q 1 , which fulfils a property that is especially helpful in the setting of the implicit Euler method. This particular choice can be taken without the loss of generality (see [7]).
Assumption 2 (Canonical projector [7]) We define Q 1 to be the canonical projector, that is, for a projectorQ 1 (y, x, t) onto ker A 1 (y, x, t), we choose Note that here the projector fulfills the following property [7] Q 1 (y, x, t)Q(t) = 0 .

Consistent Initialisation
In contrast to initial value problems arising from ODEs, DAEs cannot be initialised with arbitrary initial conditions x 0 , as the algebraic constraints imposed by the system have to be fulfilled by the solution at every time point t ∈ I. In this context we introduce the concept of consistent initial conditions. Definition 2 (Consistent initial condition c.f. [3]) Let us consider an initial value problem consisting of the quasilinear DAE (1) on the time interval t ∈ I. Then, an initial condition x 0 at t 0 is called consistent, if there exists a solution x * : I → R n of (1) fulfilling x * (t 0 ) = x 0 .
Remark 2 Higher index DAEs (≥ 2) present hidden constraints that appear only after time differentiation of the original system and thus are not explicitly accessible This complicates the choice of appropriate initial conditions, as they have to fulfil both the explicit as well as the hidden constraints.
Following [7,20,22], we make usage of the projectors defined for the tractability index to separate the degrees of freedom of the DAE system and extract the purely differential components. Their ICs are freely choosable as they are not characterised by any type of algebraic constraint. Furthermore, when they are prescribed together with the quasilinear DAE (1), a uniquely solvable initial value problem arises [22]. They can be extracted with and, if they are fixed, the rest of the components, all of them algebraic, are uniquely determined by the values of x diff and time t [7]. Note that, for index 2 systems, two projectors are required to extract the differential components. Here a further reduction within the explicitly differentiated Px components is made by further applying the P 1 (y, x, t) projector.
With a second projector T(y, x, t) onto Im Q(t)Q 1 (y, x, t), the index 2 variables of the DAE can be extracted. For the differential index, the index 2 variables are the parts of x that require one time differentiation of the original system to be characterised by (hidden) algebraic constraints [3]. Its complementary projector U(y, x, t) = I−T(y, x, t) allows to extract the purely differential components together with the index 1 variables, that is, the ones prescribed by the explicit constraints of the system. Thus we can separate the degrees of freedom into three types (see [7])

Implicit Euler
Given a quasilinear DAE (1) defined on the time interval I with consistent initial condition x 0 at initial time t 0 , it can numerically be integrated with the implicit Euler method. For time steps t 0 , t 1 , . . . , t n with t n = t end , step size t i+1 − t i = h and approximated solutions x 0 , . . . , x i , the implicit Euler method performs for the (i + 1)th time step the approximation In the following section we consider DAEs with a simplified structure. For that we impose some additional requirements on the index 2 components.

Assumption 3 (Constant projector matrices [7])
We assume the mass matrix A and the space Im Q(t)Q 1 (y, x, t) are time and space independent and thus the projectors P, Q, T and U are constant. Without loss of generality (see [7]), we consider the projector T to fulfill TP = 0. Now we study the behaviour of the implicit Euler scheme applied on index 2 tractable DAEs with linear index 2 components and constant mass matrix, that is, systems written as Assumption 4 (Constant Q 1 ) We assume the index 2 tractable DAE in (3) has a constant projector Q * 1 onto Im Q 1 (Ux, t) and define P * Remark 3 We will show later nonlinear index 2 DAEs arising from application examples that fulfil the structural assumptions made in this section to demonstrate that they are not too restrictive.
Before discussing several special properties of the implicit Euler scheme, some important characteristics of the projectors are presented.
Proposition 1 (Projectors c.f. [7]) We consider a DAE as in (3) with its corresponding projectors fulfilling Assumptions 2, 3 and 4. Then it holds that where the matrix G 2 (Ux, t) is defined as (2).
Proof Property (i) follows from the feature of the canonical projector shown in Assumption 2 and (ii) from the definition of the projector T. In [7, Chapter 2.3] Property (iii) is shown and proven and for (iv) the equivalent expression B 2 T = G 2 (Ux, t)T is derived by applying the definitions of the projectors, Assumption 3 and (i). Finally, (v) follows from the definition of Q * 1 that implies Proposition 2 (Implicit Euler consistentialisation) We consider a DAE and the corresponding projectors fulfilling Assumptions 2, 3 and 4, with structure as in (3) and two initial conditions at t 0 , the first one, x 0 , being inconsistent and the second one, x 0 , being consistent Let x 2 be the solution obtained at time point t 2 after two implicit Euler steps of the IVP (3) with inconsistent IC (4) and x 2 the one obtained with consistent IC (5) both with the same time step sizes. If both initial conditions are such that Proof The proof proceeds similarly to the approach taken in [7, Chapter 2.5] for the computation of consistent initial conditions for index 2 DAEs. The superscript i is used to denote the solutions of the implicit Euler method for time step t i starting with inconsistent initial condition x 0 and subscript i for the solutions with consistent initial condition x i .
Application of the implicit Euler scheme yields for the first time step with step size h The equations are then multiplied by G 2 (Ux 1 , t 1 ) −1 and G 2 (Ux 1 , t 1 ) −1 , respectively, which leads to due to Proposition 1 (iii) and (iv).
Now we multiply both equations by T and U to split them into the parts defining Tx 1 (respectively Tx 1 ) and Ux 1 (respectively Ux 1 ). Thus we have Application of Proposition 1 (ii), (i) and (v) together with equality PP * 1 x 0 = PP * 1 x 0 to the equations defining the U components, and PT = 0 yields (7) and Here, (7) and (9) are to be solved to obtain the solutions for Ux 1 and Ux 1 , respectively. It can be seen that both equations are equivalent and their solution only depends on PP * 1 x 0 . Therefore they yield the same solution, i.e. Ux 1 = Ux 1 .
To obtain Tx 1 and Tx 1 , (6) and (8) have to be solved. Again both equations are equivalent, however, whereas the first one depends on Ux 1 and Ux 0 , the second one depends on Ux 1 and Ux 0 . Due to the previous step, Ux 1 = Ux 1 , but, Ux 0 is not necessarily equal to Ux 0 . Thus, if Ux 0 = Ux 0 , then Tx 1 = Tx 1 and only one implicit Euler step is required to obtain the same solution with both initial conditions x 0 and x 0 .
If Ux 0 = Ux 0 , then the analogous procedure is repeated to obtain the solution for t 2 . This time, the index 2 components Tx 2 and Tx 2 are again defined by the same equation and depend on Ux 1 and Ux 2 or Ux 1 and Ux 2 , respectively. Due to the previous step, Ux 1 = Ux 1 and, analogously as before, we obtain Ux 2 = Ux 2 and thus x 2 = x 2 .
Remark 4 Note that, as the implicit Euler scheme starting with a consistent initial condition yields solutions that are consistent for semiexplicit index 2 DAEs, x 2 is consistent [3]. Therefore, as x 2 = x 2 , x 2 is a consistent solution that is obtained after two implicit Euler steps with an inconsistent initial condition. In practice this consistency is only obtained up to a certain tolerance which depends e.g. on the accuracy of the Newton scheme.
In [7] a similar result is shown for the implicit Euler scheme. However, for schemes starting with only inconsistent index 2 variables, that is, schemes where Ux 0 = Ux 0 (in fact, only Px 0 = Px 0 is required). We have extended this result to consider DAEs, where also the PQ * 1 components might be inconsistent. Furthermore we have shown that the PP * 1 components of the initial condition are "remembered" by the time integration scheme, which is a key property for the Parareal algorithm considered next in this paper.

Parareal
In the following section we will introduce the parallel-in-time method Parareal and study its application to differential algebraic equations.
Let us consider the initial value problem of the quasilinear DAE (1). To apply the Parareal algorithm, first the time interval I is partitioned into N smaller time windows I n = (T n−1 , T n ] of size ∆T = (t end −t 0 )/N , with T 0 = t 0 and T N = t end . In each iteration k Parareal solves in parallel the N initial value problems with X k 0 = x 0 and n = 1, . . . , N . This, however requires initial conditions for each subwindow I n , which are a priori unknown. Therefore, the algorithm has to start with incorrect initial conditions at the interface points T n . This yields jumps in the solution across subwindows, which Parareal tries to iteratively eliminate by updating the initial conditions. Thus, in addition to the parallel computation of the N initial value problems described previously, Parareal performs an update formula for the initial conditions which for the (k + 1)th iteration reads [21,14] X k+1 n = F(T n , T n−1 , X k n−1 ) + G(T n , T n−1 , X k+1 n−1 ) − G(T n , T n−1 , X k n−1 ) , (11) for n = 1, . . . , N . Here, F and G are solution operators of the initial value problem called the fine and coarse propagator, respectively. The fine propagator F(T n , T n−1 , X k n−1 ) returns the solution of (10) at time step T n . It solves the problem in a very accurate way and is thus computationally expensive to apply. However, as for the (k + 1)th iteration this solution only requires the initial condition computed at the previous Parareal iteration X k n−1 , it is performed in parallel. Therefore, simulation time is still reduced with respect to a sequential computation.
The second operator, G(T n , T n−1 , ) gives the solution of the initial value problem at time point T n with initial condition at T n−1 . However, X k+1 n−1 is a solution of the current Parareal iteration and thus can not be computed in parallel. Therefore, the coarse solver has to be applied sequentially. To ensure computation time is reduced as much as possible, this operator has to be cheap to compute and as a consequence less accurate. Here, for example, a larger time step size can be employed or a reduced system can be solved.
Combining the parallel computation of the expensive and accurate fine propagator with the sequential computation of the cheap coarse propagator yields each Parareal iteration to require less computation time than the sequential simulation approach. This establishes a link to other parallelisation methods such as multiple shooting methods (see [12] for an historical overiew) or multigrid approaches [10].

Parareal for DAEs
The application of the Parareal algorithm for nonlinear ordinary differential equations and its convergence is already studied [13]. However, for the case of higher index DAEs, its convergence has not been studied yet and its applicability is not ensured. Here, especially for nonlinear DAEs, it can happen that the update formula (11) yields an inconsistent initial condition for the fine solver, which may lead to divergence of the algorithm, slower convergence or an incorrect solution. Nevertheless, the algorithm has been applied to DAEs previously. For example in [4] the algorithm is applied to a system of DAEs without a special handling. In [27] it is applied to index 1 DAEs with a special structure and in [11] with a modified Parareal algorithm. Both approaches are special cases of the theory that is given in this paper.
A modification of the Parareal algorithm to be applied to quasilinear index 2 DAEs is presented here. The method is similar to the one given in [18,19] for multiple shooting methods. There it is ensured that the algorithm works by means of extending the Jacobian that is computed to update the initial conditions to avoid it being singular. This is achieved by including an equation for the calculation of consistent initial conditions. The idea behind the method here is to only apply the Parareal algorithm on the purely differential components of the DAE and then compute the rest of the degrees of freedom accordingly. For that, the update formula (11) is restricted tô whereX k n := F(T n , T n−1 , X k n−1 ) andX k n := G(T n , T n−1 , X k n−1 ). The resulting solutionX k+1 n is then used to compute the consistent initial condition X k+1 n such that This can either be done analytically for simple DAEs or with numerical techniques as proposed e.g. in [7,8]. The algorithm of [8] is implemented in a Python package called InitDAE 1 , which is able to numerically compute consistent initial conditions for DAEs. Note that is approach can also be performed on DAEs with index higher than two. In that case, the projectors to extract the purely differential components in (12) and (12) have to be adapted accordingly (see [20]).

Implicit Euler as propagator
Note that, as the update (12) in the Parareal algorithm is performed sequentially, the computation of the initial conditions X k+1 n out of the obtained solution after the updateX k+1 n is also done in a sequential manner. If this operation is computationally expensive, then it can considerably increase the simulation time of the Parareal algorithm. Furthermore, one of the advantages of the Parareal method is that, as it is not intrusive, i.e., it can even be applied to black box simulators (as long as you can prescribe initial values). In such cases, obtaining the explicit matrices of the DAE system that is solved might not be possible. These two inconveniences can be overcome for DAEs with a specific structure by means of using the implicit Euler method as a time integrator. . Implicit Euler yields a consistent solution after at most two time steps (see Proposition 2). Furthermore, its solution corresponds to the one obtained with a consistent initial condition where the PP * 1 coincide, that is, Thus, if on the fine level at least two Euler steps are performed, the solution at the end of the interval is equivalent to first computing the consistent initial condition X k+1 n and then starting the simulation with it. On the coarse level, already one implicit Euler step suffices, as only the PP * 1 components are relevant for the update and those are remembered by the time integration method also in the first time step, as shown in the proof of Proposition 2.
Therefore, for DAEs fulfilling the requirements of Proposition 2 it is not necessary to use the modified Parareal update (12) and the subsequent computation of a consistent initial condition. Here, the usage of the implicit Euler as coarse propagator and performing at least two implicit Euler steps on the fine level are sufficient for the Parareal algorithm to converge.
In [27] this property is exploited for an index 1 DAE solved with Parareal and the implicit Euler scheme.
The approach taken in [11], however, is a special case of the algorithm proposed in Section 3.1 for index 1 DAEs. There, all the components are added in the Parareal update formula (11), which is equivalent to (12) if the projector P is constant and the algebraic variables are made consistent afterwards. Both cases are covered within the formalised results of this paper, and the theory is expanded to index 2 DAEs.

Numerical Examples
In the following, two nonlinear index 2 differential algebraic equations are solved with two versions of the Parareal algorithm, the classic one and the modified Parareal algorithm for DAEs from Section 3.1.

Nonlinear Index 2 DAE
To test the proposed algorithm, we first consider the following index 2 DAE as a toy example with degrees of freedom x = (x 0 , x 1 , x 2 ) and nonlinear function Here, the projector matrices P, P 1 (x, t) are Note that the index 2 component x 2 appears nonlinearly in the DAE (14). Therefore, this example does not fulfil the requirements for the consistentialisation of the implicit Euler scheme and Proposition 3 does not necessarily hold. This implies that starting the Euler scheme with an inconsistent initial value x 0 does not necessarily yield the same solution after two time steps than starting with the consistent initial condition x 0 with the same PP 1 (x 0 , t) components.
To exemplify this behaviour a counterexample is presented. Let us consider the inconsistent initial condition x 0 = (0, −1, 0) at t 0 = 0. Its corresponding consistent initial condition x 0 with is x 0 = (0, 0, 0.3π) . After two implicit Euler steps with e.g. time step h = 1/3 and starting with the inconsistent initial condition x 0 , 045 sin(20π/3) + 3)/3 = 0 is obtained. The same scheme for initial value x 0 , however, yields the correct solution x 0,2 = 0 if initialised with the corresponding consistent value x 0 and thus x 2 0 = x 0,2 . Note that this example is an artificially created DAE without dynamics, due to the definition of the nonlinear function g(x) and the fact that x 2 is always ≤ 1.
To study the proposed modification of the Parareal algorithm, we apply the variants:

-PR Euler
In the first algorithm no special handling is implemented, that is, classic Parareal with implicit Euler as time integrator is applied with a small time step size δt for the fine solver and a larger one ∆T for the coarse propagator.

-PR Init
For the second simulation, in the update formula we only consider the PP 1 (x, t) components as in (12). Afterwards, the corresponding consistent initial condition with the same PP 1 (x, t) components as in (13) is computed.
For this particular example, the Parareal update (12), with PP 1 as shown in (16), is Once the updated value (x 0 ) k+1 n is obtained, the consistent initial condition is computed analytically with (x 2 ) k+1 n = 0.3π cos(20πT n ) (x 1 ) k+1 n = 0.015 sin(20πT n ) , and using (17) with projector matrices (16) the differential variable is obtained  For both simulations N = 21 processors are chosen and the simulation time window I = [0 1). The time step size of the fine implicit Euler propagator is set to δt = 10 −5 and the coarse solver is chosen to perform one time step per window and thus has time step size ∆T = 1/N . The Parareal algorithm is iterated until the l 2 norm of the difference between the PP 1 components of the solution at the end of interval I n−1 and the initial condition of I n for all T n is below a relative tolerance of 5 · 10 −4 and an absolute tolerance of 10 −10 (see error norm of [16]).
Whereas the first algorithm (PR Euler) requires 3 iterations to reach the required tolerance, the second algorithm (PR Init) converges immediately after the 1st Parareal iteration. The obtained solutions for the purely differential component x 0 and the index two variable x 2 are depicted in Figure 1. Here it can be seen that the proposed algorithm for DAEs 'PR Init' obtains the correct solution for the index two variable x 2 already at the first iteration. Therefore, the algorithm converges immediately, whereas the classic Parareal algorithm has jumps on x 2 that yield an incorrect solution also for the differential component x 0 . After 3 Parareal iterations the classic Parareal algorithm manages to reduce the jumps on both x 2 and x 0 but this is not covered by theory.
Even in this simple case without dynamics, where the Parareal algorithm should converge immediately, the classic algorithm requires 3 iterations due to the inconsistency of the index 2 variable, which introduces large errors in the nonlinearity term affecting thus the solution of the entire system.

Circuit with Modified Nodal Analysis
To exemplify the theoretical results presented in Sections 2 and 3.1.1 as well as demonstrate that the assumptions that are taken are realistic for real-life applications, we further present an example arising from the description of an electric network.

Modified Nodal Analysis
We consider electric networks containing capacitors (C), inductors (L), resistors (R) and voltage (V) and current (I) sources. In modified nodal analysis, networks are described by means of incidence matrices A , ∈ {C, L, R, V, I} that characterise the branch-to-node relation of the underlying graph for the corresponding elements. Applying Kirchhoff's current law and the lumped parameter models of the different components, we obtain the following system of DAEs [9,15] Here the system of DAEs is given in the flux-charge formalism. In this formulation, the additional degrees of freedom are e : I → R ne , the vector of node potentials, i : I → R n , the vector of currents through branches containing the element , q : I → R nC , the vector of charges in capacitances and φ : I → R nL , the vector of fluxes in inductances. Finally, q C (·), φ L (·), g R (·), i s (·) and v s (·) are (nonlinear) functions describing the lumped parameter relations for the different elements. The vector of node potentials and the incidence matrices allow extracting the voltages across the branches containing a given element with the relation v = A e. The tractability index of this system has already been analysed (see e.g. [9]) and is in the worst case 2. This result is given by only topological properties of the underlying graph.
Note that, unlike in the classic MNA formulation, the flux-charge system (18) yields a system of DAEs with a constant mass matrix. This is achieved thanks to the introduction of the degrees of freedom q and φ.
The index analysis in [9] shows that the possible index 2 components of the system are currents through voltage sources i V and voltages across inductances A L e. These two degrees of freedom appear linearly in the original system (18) and thus flux-charge MNA has linear index 2 components [2]. Therefore, system (18) is, in the worst case, an index 2 tractable DAE with linear index 2 components and constant mass matrix as described in (3). Finally, in [9] it is also shown that Assumption 4 is fulfilled, as Im Q 1 (Ux, t) is constant. This allows the application of Proposition 2 to the system of DAEs obtained from flux-charge MNA and thus the implicit Euler scheme returns a consistent solution after at most two time steps even if an inconsistent initial condition is given. Fig. 2: Index 2 circuit with nonlinear inductance L 2 (i L,2 ) as described in [6].
For the nonlinear inductance the model of [6] is used with nominal inductance L nom = 10 −3 H, deep saturation inductance L deepsat = 8 · 10 −4 , smoothness factor σ = 5 · 10 −2 and current I * L = 90 A. As it has been explained in Section 4.2.1, the system of DAEs arising from flux-charge MNA fulfills the requirements to apply Proposition 2. Thus, following the idea of Section 3.1.1, we do not require any special handling and can apply classic Parareal to that system of equations without any drawback.
To exemplify this, we apply again Parareal twice. Once the classic version 'PR Euler' and the second version designed for DAEs 'PR Init'. Unlike in the first example, the consistent initial conditions are not computed manually, but numerically with the Python package InitDAE. For the 'PR Init' algorithm, the projectors PP 1 * are used for the update formula. In both simulations N = 15 processors are chosen and the simulation time window is set to I = [0, 0.2). For the fine solution implicit Euler with a time step size of δt = 10 −5 is used and the coarse implicit Euler solver performs one time step per window. The initial condition x 0 is computed by starting an implicit Euler scheme with an inconsistent value x −2 = 0 at time step t 0 − 2δt and performing two steps until t 0 . Exploiting Proposition 2, the obtained solution x 0 for the rest of the simulation is a consistent initial condition. The error is computed as in the previous example and the relative tolerance is set to 10 −4 , whereas the absolute tolerance is chosen to be 10 −8 .  This time both algorithms require 4 Parareal iterations to reach the required tolerance. The result is not surprising, as both algorithms supposed to be doing the same on the differential, and thus transient part, of the system and the convergence depends only on the problem.
In Figure 3, the purely differential component φ L,1 and the index 2 component v L,2 of the fine solutions on the last two windows I n , n = 14, 15 are depicted for the first and the last Parareal iterations. It can be seen that in both cases the purely differential component has a jump at the first iteration and becomes continuous at the last one, which is a typical behaviour of the Parareal algorithm. As in the previous example, the index 2 component is immediately smooth and converged at the first iteration for the 'PR Init' algorithm, whereas 'PR Euler' starts with an inconsistent solution at the first step in the first iteration. Here the behaviour of the implicit Euler scheme when starting with inconsistent initial conditions can be observed: the solution jumps to the correct value due to Property 2. This, however, does not negatively affect the convergence of 'PR Euler', as the purely differential components are handled equally in both algorithms.
Remark 5 Note that, even for a DAE fulfilling the requirements of Proposition 2, it can happen that 'PR Init' may reach the required tolerance for the differential components, while the algebraic ones have not reached the required accuracy yet. Let us consider the solution at I n . For example, for index 1 components it suffices to take the result of the last time step of the previous window I n−1 instead of the initial value at the beginning of I n to ensure a consistent value at time T n−1 . The index 2 components would in addition require to either ignore the first (possibly) inconsistent time step of I n at T n−1 + δt (with δt being the time step size of the fine propagator), or perform one extra time step of the solution of I n−1 to arrive at the time T n−1 + δt.

Conclusions
This article has presented a modification of the Parareal algorithm for its application to quasilinear index 2 tractable DAEs. Its extension to higher index systems requires extra care, however it follows analogously from the projector-based decoupling of differential algebraic equations. For a large class of DAEs i.e. linear index 2 components and constant mass matrix as given in flux-charge formulated modified nodal analysis, a new property of the implicit Euler scheme is proven. This property allows the usage of the classic Parareal algorithm, as long as the implicit Euler scheme is used as the time integrator of both the first two time steps of the fine as well as for the coarse propagator.
The theoretical results are backed up by numerical simulations of two DAEs, one toy example with nonlinear index 2 components and the other one arising from a flux-charge modified nodal analysis formulated circuit. As theoretically expected, the modified Parareal algorithm speeds up the convergence when applied to a DAE with nonlinear index 2 components.