On Improving the Convergence of Radau IIA Methods Applied to Index 2 DAEs.
This paper presents a simple new technique to improve the order behavior of Runge--Kutta methods when applied to index 2 differential-algebraic equations. It is then shown how this can be incorporated into a more efficient version of the code {\sc radau5} developed by E. Hairer and G. Wanner.
Introduction
In recent years, differential algebraic equations (DAEs) have been studied by various authors (see
[HW91], [HLR89], [BCP91]), and their importance acknowledged by the development of specific solvers
such as DASSL from [Pet86] or RADAU5 from [HW91]. An especially important class of DAEs arising
in practice are semi-explicit systems of the form
where g y f z is assumed to be of bounded inverse in a neighbourhood of the solution of (S). Here, we
are interested in obtaining a numerical approximation to (S) accurate for both the differential and the
algebraic components. Although some of the ideas presented in this paper also apply to more general
Runge-Kutta methods, we will focus on Radau IIA methods, given that they were used to build the
code RADAU5. Their construction as well as some of their properties are briefly recalled in Section 1.1.
When applying a s-stage Radau IIA method to (S), the orders of convergence are respectively
for the y-component and s for the z-component (see [HLR89]). In some situations, where getting an
accurate value of z may be important (in mechanics for instance), one is led to use a different approach.
Generally speaking, the order reduction phenomenon may be overcome by the following techniques:
(i) a first possibility consists in applying the Radau IIA method to the index one formulation (S) of
Since Radau IIA methods applied to index one DAEs exhibit full order of convergence for y and
z (see Theorem 3-1 [HLR89]), the order of convergence is now 2s \Gamma 1 also for z. However, solving
(S) can be considerably more costly: as a matter of fact, this requires to evaluate the Jacobian
of the function F (y; at each step (or whenever the convergence rate of
the Newton iteration gets too small), instead of the function F (y; Another
drawback is that the numerical solution is not forced any longer to lie on the constraint manifold
(ii) a second idea consists in computing the z-component by solving the additional equation g y (y)f(y; z),
i.e. in projecting the numerical solution on the so-called "hidden constraint". The corresponding
numerical scheme now reads
The order of convergence for the y-component is still 2s \Gamma 1 and can be shown to be now
also for the z-component by using the Implicit Function Theorem. However, this technique is
once again computationally more demanding than the original one: solving the new implicit
part of (S n ) requires an accurate evaluation of g y (y) at each step, and not only, as previously
mentioned, whenever the convergence rate of the Newton iteration becomes too small. It can be
nevertheless noted that in a parallel environment, this additional cost would be shadowed by the
use of a second processor.
In this paper, we present a third approach which does not require an analytical form of g y and whose
computational cost is basically equal to what it is for the standard formulation. It is based on the
PI n-976
4 A. Aubry , P. Chartier
observation that the errors in the z-component are essentially of a local nature, at least up to the
order of convergence of the y-component. As a consequence, making z more accurate is a matter of
recovering the significant terms that appear in the so-called "B-series of the error". This is made
possible by considering the composition of the basic method with itself over several steps. It can be
noted that similar ideas were used by R.P.K. Chan to deal with the order reduction of Gauss methods
when applied to certain stiff problems (see [BC93]).
The new order conditions will be determined in Section 2. While they can be derived in a straightforward
manner from the work of E. Hairer, C. Lubich and M. Roche ([HLR89]), they are still relatively
unknown since they are not satisfied by most classical Runge-Kutta methods (see also [BC93] or
[CC95]). It will then be shown that some of those conditions are actually redundant and may be
omitted. This is a crucial aspect of the method, since it allows the construction of formulas with a
manageable level of complexity.
In Section 3, the implied modifications to the code RADAU5 are listed. Finally, numerical results are
presented that illustrate the advantages of this new technique.
1.1 Basic properties of Radau IIA methods
Radau IIA methods can be defined by the s
A particular method R will be characterized by the triple (A; b; c) where (a ij ) i;j=1;\Delta\Delta\Delta;s is a s \Theta s matrix,
s-dimensional vector and s-dimensional vector. In the sequel,
we will furthermore use the notations
1.1.1 Construction of Radau IIA methods
Their coefficients are uniquely determined by the following conditions:
1. are the ordered zeros of the Radau right polynomial
dx
x
2.
3. The coefficients a ij of the matrix A satisfy C(s):
As the c i are all distinct, A is non-singular.
1.1.2 Some useful properties
We will refer here to the additional simplifying assumptions D(-) introduced, as B(p) and C(j) of
previous subsection, by J.C. Butcher (see [HW91] page 75).
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 5
1. Due to the conditions C(s), B(s) and c methods are stiffly accurate (i.e.
vectors b and (a are solution of the same Cramer
system).
2. They are collocation methods (Theorem 7-8 [HNW91], page 212).
3.
4. D(s \Gamma 1) is satisfied (Lemma 5-4 [HW91], page 78).
Optimal convergence results have been obtained for those methods by J.C. Butcher on the one hand
(Theorem 5-3 of [HW91]) and E. Hairer, C. Lubich and M. Roche on the other hand (Theorems 3-1,
4-4 and 4-6 of [HLR89]). They are collected in Table 1.

Table

1: Optimal global error estimates for the s-stage method Radau IIA
1.2 Increasing the order of convergence of the z-component
When applying a s-stage Radau IIA method c) to the system (S), we obtain
z
are the coefficients of the matrix A \Gamma1 . In (4), z n vanishes because us
now replace the vector (b T A by an adjustable vector (4). By doing so,
we define a new method Rw . It is easily seen that the order of convergence of the y-component
remains unchanged. As for the z-component, the lack of accumulation makes the errors purely local.
The convergence behaviour of the z-component is thus entirely determined by the following order
conditions from Theorem 8-6 and 8-8 of [HW91].
Proposition 1 Let ffi y and ffi z be the local errors respectively for the y and z components of a Runge-Kutta
method. Then we have,
where DAT2 y , DAT2 z are sets of trees, \Phi a vectorial function and fl; ae scalar ones associated with
the trees 1 .
Nevertheless, w has not enough components to allow for an order of convergence greater than s
is the optimal vector). Hence, to get sufficient freedom, we need to consider the
composition R oe of R over oe steps. As variable steps h i , are considered, we also have to
introduce the ratios r
characterized by the triple (A; B; C) where
ffl A is the blockmatrix blocks of the form
1 See [HW91] for a definition of these notions.
PI n-976
6 A. Aubry , P. Chartier
ffl B is the blockvector (B i ) i=1;\Delta\Delta\Delta;oe with blocks of the form B
ffl C is the blockvector (C i ) i=1;\Delta\Delta\Delta;oe with blocks of the form C
e.
Replacing
by an adjustable vector w offers oe \Theta s degrees of freedom, i.e hopefully enough
for oe ? 1 to increase the order of convergence of the z-component. It is our aim now to show how to
construct w and how to implement the new method R oe;w .
2 Construction of the vector w
In this section, the effective construction of w is described. It should be emphasized that its components
depend on r forcing one evaluation per step. However, these additional computations
become negligible as soon as the dimension of the system (S) is large enough.
2.1 Order conditions
The conditions for order k are enumerated below together with the associated trees.
required with
z -
are required with
s
z -
z -
Let U s be AC
s+1 C s+1 . If (s) is satisfied, then (C s+1 ) is equivalent to
are required with 2
z -
z -
s
z -
z -
s
z -
Let U s+1 be AC
equivalent to
2 The dot stands for the componentwise product.
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 7
are required with
z -
z -
z -
z -
z -
s
z -
z -
s
z -
s
z -
z -
s
z -
s
z -
z -
Let U s+2 be AC
s+3 C s+3 . If (s + 7) is satisfied, then (C s+3 ) is equivalent to
These conditions are obtained by Proposition 1 and by using simplifying assumptions (it is important
to note that the composite method R oe satisfies B(2s \Gamma 1), C(s) and D(s \Gamma 1)).
2.2 Preliminary calculus
In order to later simplify the equations for w, we now state some basic results.
eb T , then F
Proof: By definition, we have eb T A \Gamma1 . As the method c) is stiffly accurate,
ve T
s where ve T
eb
ve T
ve T
ve T
ve T
be the s \Theta s blocks of the matrix A \Gamma1 , then
PI n-976
8 A. Aubry , P. Chartier
Proof: By definition, AA
Lemma 2 is then easily proved by induction on i. 2
oe;n
be the oe \Theta s dimensional vector AC
n+1 C n+1 and u n be the
s dimensional vector
is given by
k=s
r k+1
Proof: Let n be less than or equal to 2s \Gamma 2. For all
S i;n can be expanded as follows
r k+1
r k+1
e:
Owing to C(s), we have
k=s
r k+1
e:
Similarly, T i;n can be expanded as
l
r l
implies that
and as n - 2s \Gamma 2, we obtain
r l+1
e
e
e
On improving the convergence of Radau IIA methods applied to index 2 DAEs 9
Lemma 4 For all n less than
Proof: This follows at once from B(2s \Gamma 1). 2
Lemma 5 Let
oe;n
be the vector A \Gamma1 U n , then for all integer n -
by
k=s
Proof: By definition, X
A \Gamma1 U i;n . Using
Lemma 3, we obtain
k=s
r k+1
k=s
and we use Lemma 4 to complete the proof. 2
Lemma 6 For all n less than 2s
Proof: This follows straightforwardly from the order conditions for the trees
z -
z -
1. Let
oe;n
be the vector C:A
k=s
r k+1
2. Let
oe;n
be the vector A \Gamma1 (C:U n );
k=s
r k+1
Proof: By definition, Y and the first part of the result is obtained by applying
Lemma 5. Now, let
oe;n
denote the vector C:U n . From Lemma 3, we can write T i;n as
k=s
r k+2
We have furthermore Z
so that Lemma 2 leads to
k=s
s n\Gammak
s n+1\Gammak
k=s
r k+1
PI n-976
A. Aubry , P. Chartier
The result then becomes a consequence of lemmas 4 and 6. 2
Lemma 8 For all n less than 2s
Proof: This follows from the order conditions for the trees
z -
z -
z -
Proof: The results follow from the order conditions for the trees [[-; [
z -
z -
z -
z -
z -
2.3 Results for the 2-stage method
In order to get a third-order method, w has to satisfy the following linear system (S L;2
Taking leads to a system with 4 equations and 4 unknowns. For convenience, we recall below
the coefficients of the 2-stage Radau IIA method,
R =312
The matrix M corresponding to (S L;2 ) is then of the form
and we have
Hence, for all r 1 2 (0; 1), M is non-singular and (S L;2 ) has the following unique solution
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 11
2.4 Results for the 3-stage method
In order to get a fifth-order method, w must satisfy the following linear system (S L;3
Taking now leads to a system with ten equations but only nine unknowns. However, we will show
in this section, that one of these equations is identically satisfied. Collecting all results of Section 2.2,
we get
r 4
r 4
r 4
3
A
r 4
r 4
r 4
r 4
r 4
r 4
r 4
r 4
r 4
It is found that
r 4
r 4
r 4
r 4
r 4
r 4
and the system (S L;3 ) is equivalent to
Theorem are linearly dependent.
Proof: It is enough to show that v 1 , v 2 and u 3 are linearly dependent. Since c) is of order
5 for the differential component, we have
so that b T u implies the result.2
the vectors A \Gamma1 U 3 , U 3 and V 1 (for example) become linearly
dependent. To prove this, it is sufficient to show that A \Gamma1 u 3 , u 3 and v 1 are dependent. As b T A \Gamma1 u
PI n-976
A. Aubry , P. Chartier
0, we can conclude as in Theorem 1. In this case, w depends on one parameter
that can be chosen so as to minimize the quantity
where DAT2 z 5g. This seems a natural goal to achieve, since this is an
attempt to minimize the local error. For convenience, Table 2 collects the trees of DAT2 z (5) and the
values of the associated functions ff, fl and \Phi. To compute the ff's, we refer to [Hig93].
tree u ff(u) fl(u) \Phi(u)

Table

2: Trees of DAT2 z (5) and their associated functions
2.5 Sketch of the case
To achieve order 7 for the algebraic component, we have to solve the system (S L;4 ) composed of (C 4 ),
Section 2.1), that is to say 24 equations. Comparing the number of equations
and the number of unknown, we could consequently think of taking In fact, it can be shown
that is sufficient, since 5 equations are identically satisfied (see Section 6.1). However, it does
not seem reasonable any more to consider a practical implementation of the corresponding method,
owing to the complexity of the formulas for variable stepsize.
Remark 2 If the stepsize is constant (i.e. r equations are identically satisfied
(see Section 6.1), and oe can be chosen equal to 4.
3 Modifications to the code RADAU5
The 3-stage Radau IIA method has been implemented by E. Hairer and G. Wanner in order to solve
problems of the form MY and DAEs of index less than or equal to three can be solved
by this code, called RADAU5. A precise description is given in Section IV-8 [HW91] and we will adopt
the notations used there. Implementing our method requires slight modifications to the subroutines
"radcor" and "solout" which are actually replaced respectively by "radcorz" and "soloutz".
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 13
3.1 Modifications to radcor
Only the computation of the algebraic component (z) is modified (if an index 2 DAE is solved). Once
the n th step has been accepted, two cases are considered:
1. less than three steps have been computed. Then, we keep the internal stages z n;1 , z n;2 , z n;3 and
the step size h n . The value of (y
2. three or more steps have been computed. Then, we keep the internal stages z n;1 , z n;2 , z n;3 and
the step size h n . r 1 , r 2 are computed by the formulas
and r
and w by the subroutine "vectw2". For (y n+1 ; z n+1 ) we put
z
Remark 3 For continuous outputs, we need also to keep the internal stages over two steps for the
differential components (see section 3.2.2).
3.2 New subroutines
3.2.1 Vectw2
In order to compute the formal expression of w and to create the associated fortran subroutine, the
manipulation package Maple was used. A call to vectw2 uses the format vectw2(icas,vw,r1,r2),
where the inputs are one of the five cases described below (icas) and the parameters r 1 , r 2 (r1,r2).
The output is the vector w (vw). Five cases are considered:
1. 1. In this case, we have seen that w depends on one parameter. It is optimized
as explained in Remark 1.
2. r
3. r
4. r
5. r 1
This allows us to reduce the cost of computation and to eliminate computational problems: had we
used the general expression of w (case 5), divisions by zero would have occured in the cases 1 to 4.
3.2.2 Continuous outputs
In the code Radau5, the subroutine "solout" provides the user with approximations at equidistant
output-points. The corresponding interpolation formulas are implemented in the subroutine "contr5".
"contr5(I,x)" gives an approximation U I (X) to the I th component of the solution Y at the point x
should lie in the interval [x n ; x n+1 ]). U is the collocation polynomial: it is of degree 3 and defined by
3:
PI n-976
14 A. Aubry , P. Chartier
For index 2 DAEs, are polynomials of degree 3 which satisfy
3:
By Theorem 7-8 ([HW91]), we have
As our aim is to increase the order of convergence for the algebraic component, it seems natural to
search for an approximation P
Approximation of z(x)
Let x be of the form x n\Gamma2 We define q as follows
where the vector satisfies the linear system
Using the notations of section 2.4, we have
Proposition 2 For all ' 2 (0; 1], S z
L;3 (') possesses a solution and
Proof: From Theorem 8-5 and 8-6 in [HW91], we have
According to the analysis of section 2.1, this leads to the system S z
L;3 (') which, by Theorem 1, possesses
a solution. 2
The subroutine "vectwz" computes w('). As for the vector w, five cases are considered. When h
depends on one parameter (case 1) which is not adjusted as in Remark 1. In
this case, the value defined by continuity for w(') is choosen. A call to "vectwz" uses the format
vectwz(icas,vwz,r1,r2,t), where the inputs are one of the five cases described before (icas) and
the parameters r 1 , r 2 , ' (respectively r1,r2,t). The output is the vector w(') (vwz).
Approximation of y(x)
Let x be of the form x
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 15
where the vector satisfies the linear system
C) is the Runge-Kutta method R 2 obtained by the composition of the 3-stage Radau IIA
method c) over two steps,
re
and E is the vector (z -
Proposition 3
Proof: Let us introduce the vector and the following polynomial
u:
From Theorem 8-5 and 8-6 in [HW91], we have
Using the points of the internal stages, it follows
so that (P ) is equivalent to the system S y
The subroutine "vectwy" computes B(j). A call to "vectwy" uses the format vectwy(vwy,r,t),
where the inputs are the parameters r and j (respectively r,t) and the output the vector B(j) (vwy).
The subroutine "soloutz" provides the user with approximations (p(x i
out
out
out
out
at equidistant output-points
out ) i=1;\Delta\Delta\Delta;N . For the differential component, x i
out is of the form x
choosen so as to satisfy x
out - x n+1 and for the algebraic one, x i
out is
of the form x n\Gamma2 is choosen in satisfy x
out - x n . A call to
"soloutz" uses the format
where the inputs are the number of accepted steps (nr), x n (xold), x n+1 (x), (y n+1 ; z n+1 ) (y), the
system's dimension (neq), one of the five cases described before for the computation of w(') (icas),
the parameters r1, r2 (r1,r2), the stepsize h n\Gamma2
variable (last) to indicate if the last computational step has been reached.
PI n-976
A. Aubry , P. Chartier
significant digits
modified radau5

Figure

1: Precision versus computing time - algebraic components - test problem
4 Numerical experiments
4.1 Test problem
We consider the index two problem:
where \Psi is the following infinitely smooth function:
with consistant initial values. The exact solution is
In both codes, we set
(work(4) is the parameter - in the stopping criterion for Newton's method. work(4), work(5) are
the parameters c 1 , c 2 in the stepsize control, see [HW91], page 130-134).
In figure 1, we plot the cpu time against the number of significant digits (\Gammalog 10 (absolute error)) of
the algebraic components, for both codes. For this, we use continuous outputs : outputs are required
at we compute the global error and then take the maximum over all values. In figure 2,
we plot the cpu time against the number of significant figures of the differential components, for both
codes.
In the following problems, only the modified parameters will be shown.
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 17261014
significant digits
modified radau5

Figure

2: Precision versus computing time - differential components - test problem
4.2 Pendulum
The simplest constrained mechanical system is the pendulum, whose equations of motions are described
in [HW91], page 483-485. We have applied the code Radau5 and the modified code to the GGL (Gear,
Gupta and Leimkuhler) formulation
with consistant initial values
simplicity, we took
In figure 3 (respectively 4), we plot the cpu time against the number of significant digits of the
algebraic (resp. differential) components, for both codes.
4.3 Multibody mechanism
A seven body mechanism is described in [HW91], page 531-545. We have applied the code Radau5
and the modified code to the index 2 formulation
with consistant initial values.
In figure 5 (respectively 6), we plot the cpu time against the number of significant digits of the algebraic
differential) components, for both codes. Here, outputs are required at t
4.4 Discharge pressure control
This simplified model of a dynamic simulation problem in petrochimical engineering is described
in [HLR89], page 116-118. We have applied the code Radau5 and the modified code to the following
PI n-976
A. Aubry , P. Chartier246810
significant digits
modified radau5

Figure

3: Precision versus computing time - algebraic components -
significant digits
modified radau5

Figure

4: Precision versus computing time - differential components - pendulum
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 19
significant digits
modified radau5

Figure

5: Precision versus computing time - algebraic components - seven body mechanism
significant digits
modified radau5

Figure

versus computing time - differential components - seven body mechanism
PI n-976
A. Aubry , P. Chartier13579
significant digits
modified radau5

Figure

7: Precision versus computing time - algebraic components - discharge pressure control26100.01
significant digits
modified radau5

Figure

8: Precision versus computing time - differential components - discharge pressure control
with consistant initial values. Here, initial step size h is equal to 10 \Gamma5 .
In figure 7 (respectively 8), we plot the cpu time against the number of significant digits of the
algebraic (resp. differential) components, for both codes. Here, outputs are required at t
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 21
5 Conclusion
A new simple technique to overcome the order reduction phenomenon, appearing for the algebraic
component when Radau IIA methods are applied to index two DAEs, is proposed.
Increasing the order of convergence of z is made possible by considering the composition of the basic
method with itself over oe steps. As z n+1 is defined in the basic method as a linear combination of the
internal stages, a good choice of oe should provide enough freedom for the order conditions associated
with the composite method to be satisfied. We have determined these order conditions which derive
straightforwardly from the work of E. Hairer, C. Lubich and M. Roche (section 2.1). Then we have
shown that some of those conditions are redundant and might be omitted for s-stage Radau IIA methods
with s - 4 (section 2.3 to 2.5).
It could be interesting to generalize this simplifications to any Radau IIA methods. A general question
will be to determine how many compositions of a s-stage Radau IIA method have to be considered to
obtain an order of convergence equal to 2s \Gamma 1 for the algebraic component. However, it does not seem
reasonable any more to consider a practical implementation of the corresponding method for s - 4,
owing to the complexity of the formulas for variable stepsize.
The formulas for incorporated in the code Radau5 developed by E. Hairer and G.
Wanner. Slight modifications were required (section 3). Only the computation of the algebraic component
was modified and a new procedure in order to have continuous outputs was created where we
have used our technique for both components (algebraic and differential) to compute approximations
of order five at equidistant output-points.
According to our numerical experiments, results for the differential components are disappointing.
However, the use of our technique in the code Radau5 leads to an increase of the accuracy for the
algebraic components (when tolerances are sufficiently small).
6

Appendix


6.1 Construction of the vector w in the case
In this section, we explained the calculus of linear algebra used to show that composed five times the
s-stage Radau IIA method is sufficient in the case expand the following vectors of (C 6 ) (We
use Lemma
A
A
A
and the following vectors of (C 7 ) (We use Lemma
PI n-976
22 A. Aubry , P. Chartier
A
r 6
r 6
A
A
A
A
r 6
r 6
r 6
r 6
Let introduce the following vectors
(idem for the vectors t
r 6
r 6
r 6
r 6
r 6
r 6
r 6
and the system (S L;4 ) is equivalent to
Irisa
On improving the convergence of Radau IIA methods applied to index 2 DAEs 23
Proposition 4
1. are linearly dependent.
2. are linearly dependent.
3. are linearly dependent.
4. are linearly dependent.
5. are linearly dependent.
Proof: Because of the expression of the vectors it is sufficient to show that
are linearly dependent (idem for the part 2 to 5 of the proposition). The method
c) is of local order 8 for the differential component. Thus, order conditions associated with
the trees of DAT2 y (7) are satisfied. In particular,
Finally, we obtain
but b 6= 0, hence the proposition is shown. 2
Remark 4 If the step size is constant (i.e. r
1. are linearly dependent.
2. are linearly dependent.
3. are linearly dependent.
4. are linearly dependent.
Hence, nine equations are identically satisfied and p can be choosen equal to four.
PI n-976
A. Aubry , P. Chartier



--R

On smoothing and order reduction effects for implicit Runge-Kutta formulae
Numerical solution of initial value problems in differential-algebraic equations
A Composition Law for Runge-Kutta Methods Applied to Index-2 Differential-Algebraic Equations
Coefficients of the Taylor expansion for the solution of differential-algebraic systems
The Numerical Solution of Differential Algebraic Systems by Runge-Kutta Methods
Solving Ordinary Differential Equations (Vol
Stiff Problems and Differential Algebraic Problems (Vol.
A description of DASSL: A differential/Algebraic System Solver.
--TR

--CTR
Frank Cameron , Mikko Palmroth , Robert Pich, Quasi stage order conditions for SDIRK methods, Applied Numerical Mathematics, v.42 n.1, p.61-75, August 2002
Frank Cameron, A Matlab package for automatically generating Runge-Kutta trees, order conditions, and truncation error coefficients, ACM Transactions on Mathematical Software (TOMS), v.32 n.2, p.274-298, June 2006
rooted trees;runge-kutta methods;composition;differential-algebraic systems of index 2;Radau IIA methods;simplifying assumptions