BIT36 :1(1996),1-13 . STABILIZEDMULTISTEPMETHODSFORINDEX 2EULER-LAGRANGEDAEst C .AREVALO 1 , C .FUHRER2 andG .SODERLIND 3 1 Dept .ofPureandAppliedMathematics,UniversidadSimonBolivar Apartado89000,Caracas1080-A,Venezuela, e -mail :camena@cesma .usb .ve 2 Inst .forRoboticsandSystemDynamics,GermanAerospaceResearchEst .(DLR) D-82230Wessling,Germany . e-mail:claus .fuehrer@dlr.de 3 GustafSoderlind,Dept .ofComputerScience,LundUniversity,Box118 S-22100Lund,Sweden . e-mail :gustaf@dna .lth .s e Abstract . Weconsidermultistepdiscretizations,stabilizedby0-blocking,forEuler-Lagrange DAEsofindex2 .Thuswemayuse"nonstiff"multistepmethodswithanappropriate stabilizingdifferencecorrectionappliedtotheLagrangianmultiplier term.Weshow thatorder p= k-1-1canbeachievedforthedifferentialvariableswithorder p=k for theLagrangianmultiplierfork-stepdifferencecorrectedBDFmethodsas wellasfor loworderk-stepAdams-Moultonmethods . Thisapproachisrelatedtotherecently proposed"half-explicit"Runge-Kuttamethods . AMSsubjectclassifications :65L05 . Keywords : differentialalgebraicequations(DAE),Euler-Lagrangeequations, mul- tistepmethods,,Q-blockedmethods,partitionedmethods,compoundmultistep meth- ods . 1 Introduction. Weshallconsidermultistepdiscretizationsofthe index2Euler-Lagrange DAEs where f : Rn-~ltn,g ; IR, and G(x)=8g/8x, with G(x)GT'(x) invertible. Equationsofthistypeoccurinthesi ulationofmultibodydyna When suchsystemsaresolvednumericallybymultistepmethods,itiscommontouse thebackwarddifferentiationmethods(BDF),alsointhecasewhenthe state- spaceformoftheproblemisnonstiff .Theapplicationofmethodssuitable for nonstiffproblemshasbeenstudiedrecently,seee .g .[3,4] . Inthecaseofmultistep methods,thisentailscombiningdifferentmultistepformulasfor differentparts ADecember1994.RevisedOctober1995 .
Transcript
BIT 36:1 (1996), 1-13 .
STABILIZED MULTISTEP METHODS FOR INDEX2 EULER-LAGRANGE DAEs t
C . AREVALO 1 , C. FUHRER2 and G. SODERLIND 3
1 Dept. of Pure and Applied Mathematics, Universidad Simon BolivarApartado 89000, Caracas 1080-A, Venezuela, e-mail: camena@cesma .usb .ve
2 Inst . for Robotics and System Dynamics, German Aerospace Research Est . (DLR)D-82230 Wessling, Germany . e-mail: claus [email protected]
3 Gustaf Soderlind, Dept . of Computer Science, Lund University, Box 118S-221 00 Lund, Sweden . e-mail: gustaf@dna .lth .s e
Abstract .
We consider multistep discretizations, stabilized by 0-blocking, for Euler-LagrangeDAEs of index 2. Thus we may use "nonstiff" multistep methods with an appropriatestabilizing difference correction applied to the Lagrangian multiplier term. We showthat order p = k -1 -1 can be achieved for the differential variables with order p = k forthe Lagrangian multiplier for k-step difference corrected BDF methods as well as forlow order k-step Adams-Moulton methods . This approach is related to the recentlyproposed "half-explicit" Runge-Kutta methods .
We shall consider multistep discretizations of the index 2 Euler-LagrangeDAEs
where f : Rn -~ ltn, g ; IR,
and G(x) = 8g/8x, with G(x)GT'(x)invertible.
Equations of this type occur in the si ulation of multibody dyna Whensuch systems are solved numerically by multistep methods, it is common to usethe backward differentiation methods (BDF), also in the case when the state-space form of the problem is nonstiff. The application of methods suitable fornonstiff problems has been studied recently, see e .g . [3, 4] . In the case of multistepmethods, this entails combining different multistep formulas for different parts
A December 1994. Revised October 1995 .
2
C. AREVALO, C . FUHRER AND G. SODERLIND
of the system. In [3] we studied how to apply these compound methods to theproblem with a linear algebraic constraint . This paper continues the study ofconvergence for such methods . Although the discretizations considered here areadditive linear multistep compounds they are still strongly related to some ofthe discretizations of the more general DAE systems considered in [1, 2] .
We introduce the generating polynomials p(() = Eo cxjcj and a(() = E0 i3j ;Jdefining a linear multistep method (p, a) . Let E denote the forward shift oper-ator corresponding to a stepsize h . For notational convenience, we let p and adenote the two difference operators p = E -k p(E) and a = E-ka(E) , respec-tively, i .e. for a sequence {xn } we have
k
k
Pxn = 1: ak-ixn-i ii=0
axn = > (3k-ixn-ii=0
With this notation, compound discretizations can be conveniently expressedwhile standard notions such as e .g. "the root condition for the polynomial p"still retain their classical meaning . We also want to modify a by adding abackward difference stabilizing term, i .e. we will consider operators of the forma+ CV k, where V = 1- E-1. We further assume that the method is normalizedso that al = 1. The coefficient c is selected with respect to the normalizationto make the compound discretization well defined .
The usual linear multistep discretization of (1.1) is
(1.2a)
h-1Pxn = a(f(x,,) - GT (xn)An)
(1 .2b)
0 = g(xn ) .
It is well-known that a must satisfy the strict root condition, [7, p . 498], forthis discretization to be convergent of order p when (p, a) is a p-th order stablelinear multistep method. This effectively limits the useful methods to a rathernarrow class, with the BDF methods as the most common choice . However,compound discretizations offer a wider range of feasible methods, [1, 2, 3] . Weshall consider the approach introduced in [3] . The strong method requirementsreferred to above are due to the way in which the multiplier appears in (1.1) .By using a compound method, it is possible to let the discretization (p, a) bebased on a "nonstiff method" if the multiplier term is discretized separately (0-blocking) with a difference operator providing sufficient stability . Although lessimportant, it is even possible to consider explicit methods provided that a(C)satisfies the strict root condition as a polynomial of degree k - 1, [2] . The lattertechnique is then a linear multistep analogue of the half-explicit Runge-Kuttaand extrapolation methods introduced in [4, 8, 9] .
The aim of the paper is to prove convergence results for multistep methodsstabilized by /3-blocking . We focus on Adams and difference corrected BDFmethods. At this stage we refrain from a full analysis of variable-step aspectsand finite-h stability properties . Instead we indicate that our methods are ofpractical value by solving a few test problems arising in applications . The com-putational results confirm the theoretical convergence orders .
STABILIZED MULTISTEP METHODS
3
The paper is built up as follows . In Section 2 we introduce /3-blocked multistepmethods (using backward difference stabilization) for index 2 Euler-LagrangeDAEs and prove our main result on the convergence order of the proposed meth-ods. We also formulate the criteria a compound method must satisfy in orderto be useful for (1.1). Then, in Section 3 we discuss the (nonstiff) difference-corrected BDF (DCBDF) methods [10] and show that they fulfill the necessaryrequirements if applied properly. In Section 4 we discuss using Adams methods,and why this approach is less successful . Finally, in Section 5, we give practicalcomputational results verifying theoretical order results as well as demonstrat-ing the feasibility of using nonstiff /3-blocked multistep methods for index 2Euler-Lagrange DAEs .
2 Convergence of difference corrected multistep methods .
Let us briefly consider a simple one-step compound discretization to introducethe idea of stabilizing a discretization by means of /3-blocking . Consider first(1 .1) discretized by the trapezoidal rule (this may be considered either as theone-step Adams-Moulton method or as the one-step difference corrected BDFmethod), for which o'TR = ( 1 + E-1 )/2,
(2.1a)
h -1Vx,, = o"TR(f(xf.) - GT(xn)A )(2.1b)
0 = g(xn,) .
Here QTR fails to satisfy the necessary stability requirements as stated in [7, p .498] . A minor modification, however, will stabilize the discretization. Thus weadd a stabilizing term to the right-hand side of the first equation :
(2.2a)
h-1Vx,,= oTR(f(xn) - GT(xn)A ,, ) - GT(x, )VA, /2
(2.2b)
0 = g(x,).
Since A (t) is continuous, the modification of the right-hand side introduces an0(h) deviation; the consistency of this discretization therefore drops by oneorder of h . Collecting terms, we obtain
The last term of (2 .3a) accounts for the "curvature" of the constraint manifold .The stabilizer could be introduced in several different ways, but the one proposedhere turns out to be particularly useful as will be seen below. The effect ofchanging the discretization here results in A being treated by the implicit Euler(BDF1) method. This recovers the stability of the overall discretization, i .e .(2.3) is a stable method for (1 .1) . This stabilization is in part bought at theprice of a reduced order of convergence, since the order of consistency of (2.2)is obviously lower than that of (2.1) . Here, however, as well as in the generalcase treated below, convergence order reduction only affects A . We refer to the
4
C. AREVALO, C . FUHRER AND G. SODERLIND
stabilization technique as "o-blocking" ; if one applies a method with a a notsatisfying the strict root condition to the given index 2 problem, the instabilitydue to the coefficients Oi of a can be blocked by employing a suitable stabilizingdifference correction to the A variables .
Now for the general case, consider a method (p, a) with a not satisfying thestrict root condition . Applied to the given index 2 problem, the method isnumerically unstable. This instability is associated only with the Lagrangianmultipliers . Thus it is confined to a subspace [3] defined by the structure of theEuler-Lagrange index 2 equation, and not by the individual problem, so that anaction taken to stabilize the method must be based on structure properties . Forlinear DAEs this was studied in [3] . Below we generalize to nonlinear systems andshow that it suffices to ask for strict stability only in the subspace of Lagrangianmultipliers . To this end we introduce an additional operator -r(0 = Eo -(l(3,acting on the algebraic variables only :
(2 .4a)
h-'pxn = o-(f (xn) - GT (xn)An) -I- GT (xn)rAn
(2 .4b)
0 = g(xn ) .
The stabilizer T is selected such that two objectives are met . First, we ask thata-T be strictly stable, and second, that r = O(hq) with q large or maximal. Asthe operator r affects consistency, we ask that the order drop be minimal . It isoptimal to have q = k, i .e. for some parameter c, we ask r = cVk . It depends on(p, a) if such an optimal choice can be made . For DCBDF methods, the obviouschoice of T produces no order reduction in x and a loss of one order in A. In thek step Adams-Moulton case, appropriate stabilizers achieving the correspondingresult are found for k < 3 . If further order drops are acceptable, then we canalso stabilize AM4.
The following theorem generalizes the corresponding result in [7] to the classof f3-blocked methods .
THEOREM 2 .1 . Consider the index 2 problem (1 .1) discretized by (2.4) withstarting values xi -x(ti) = O(hP) and Ai -,X(ti) = O(hk) for i = 0 : k - 1, wherep = k or p = k + 1 is the order of the pair (p, o - ) . Furthermore, let 'r be chosensuch that a -'r has roots only inside the unit circle and ry(t n ) = O(hk) for any
sufficiently smooth function . Then the discretization (2.4) is convergent, i .e. fornh = t constant,
xn - x(tn) = O(hP) and An - \(tn) = O(hk) .
PROOF. Let Si = fi - yi, and normalize coefficients such that bk = 1 . Theglobal errors of x and A at time to are denoted by ex and en . We have
where u = O(hP+1 ) . From this equation we subtract
Pxn - ha(f (xn) - GT (xn )An ) - hGT (xn)rA = 0
to obtain, after linearizing and assembling all higher order terms in s(l) and sn2) ,
where
C(h) _
En
hEn
STABILIZED MULTISTEP METHODS
5
(2.5a) pen = ho (ff rn - GTEn - GT A nEn) + hGT (-ren + 0(hk)) + hs(l) + u
(2.5b) GEn = s~2 >,
where Gx = 8G/ax,
Is(') I = cl((IEjI+IEnl) 2 +O(hk)IEnI),
Is$2) I = c2 1EX 1 2 .
Define Aj (h) = ajI - h/3j f,~ + h,3j (GxA) and let this matrix be evaluated at thepoints tn_k+j . Introduce the matrices
Q(h) = (GAk 1 (h)GT) -1 G,
P(h) = I - Ak1 (h)GTQ(h),
and note that P is an oblique projector onto the nullspace of Q, i .e . QP = 0 .Using PAk 1 GT = 0 (for all n and h), we may write (2 .5) in companion matrixform, [3], as
Now, using (2.5b) and A(h) = A(0) + O(h), 13(h) = 13(0) + O(h) yields
En
_ A(0) + O(h)
O(h)
En_1hEn
O(h)
13(0) + O(h) ( hEn-1
Making use of the strict root condition for Q-T, we obtain, with suitably chosennorms,
IIEnx II
_~
1 + 0(h)
0(h)
IIEn-111
M1
hIIEnII
O(h)
rb + 0(h)
hIIE2-1II
+
M2
'
where M1 = 0(hP+ 1)+0(hIs(1) I+Isn ) I) and M2 = O(hk+1 )+O(hls(1) l + IS(2) i)Let IS* I be the maximum modulus zero of o - r . Then r6 = I(* I if (* is simple,otherwise r6 = 1(* l + e for some e > 0. For BDF and /3-blocked DCBDF all rootsare 0 and rb = e. We can diagonalize the matrix in the difference inequality sothat
and
to get
Thus,
IIEnII
)hIIE,II
n
rb + O(h)) T - ( 0(h) 0(1) )
1
rb + O(h) )
( 00(1) ) 0(1)
O(IIEo II) + O(h 2 IIE0 II) )0 (hII Eo II) +O(hII Eo II)
O(hP) + klls(1)I +k2h-11s(2)I /
O(hk+1 ) + k3hIS(1) I + k41s(2) I
I ,n I
0(h P ) + kl ls(1) I +k2h-1 Is (2) l
lent < O(hk ) + k3 ls(1) I + k4h-1 Is (2)
Theorem 6 .2 in [7, p .492] can be extended for methods of the type (2 .4) . Byapplying this theorem we get
Is$2) 1 < CIs 2)1I, i = 1,2
with a constant C depending only on f and g and derivatives thereof .
From the definition of sn2) we obtain, by a simple induction over n, s~l) = O(hP)
and s~2) = O(hP+1) . Thus we conclude
STABILIZED MULTISTEP METHODS
7
,-I < O(hP), IEnl < O(hk ) .
0REMARK 2 .1 . If T and GT (xn ) are commuted in the discretization (2.4),
i .e .
(2.6a)
pxn = ha(f (xn) - GT (xn )An ) + TGT (xn)An
(2.6b)
0 = g(xn)
and all the initial values are of 0(hk), then the method is convergent of 0(hk)
for all variables, regardless of whether (p, a) has p = k or p = k + 1 .
REMARK 2 .2 . It is possible to consider an extended version of Theorem2.1, where 'ry(tn ) = O(hq), for a general q not related to k . In this case, ifstarting values satisfy x i - x(ti) = O(h`I"n(P,4+1)) and Ai - A(ti) = O(hmin(P , q ))
the method is convergent, and xn, - x(tn) = O(hi" in(P , q+l )) and An - .X(tn ) _O(hmin(P:q)) .
3 Difference corrected BDF methods .
The difference corrected BDF methods were introduced in [10] . Let (pk,1)denote the BDFk, where an ODE x = f (x) is discretized by h- ' pkxn = f (xn ) :
This discretization is consistent of order 0(hk) with an error whose principalterm is h-' Vk+1x(tn)/(k + 1) ;Z:~ Vkx(tn)/(k + 1) = vkf (x(tn))/(k + 1) . Theprincipal error term can therefore be recovered by a difference correction appliedto the right-hand side,
_?
xn=(1 -k+1 )f(xn)
We define vk = 1 - 7k 1(k + 1) . The method (pk, Qk) constitutes the DCBDFkmethod; it only differs from the BDFk by the difference correction -Vk/(k+1)to the QBDF - 1 operator. In linear multistep implementation for ODEs it isconvergent and of order p = k + 1 for k < 6, due to the structure of pk . In
8
C. AREVALO, C . FUHRER AND G . SODERLIND
their one-leg [5] form, the order of convergence drops to p = k for nonlinearODEs. The DCBDF methods have bounded stability regions, hence they areintended for nonstiff integration, cf. Fig. 3.1 . A variable step formulationof these methods is given in [10] . They have larger error constants than theAdams-Moulton methods of the corresponding order, cf. Table 3 .1, but theirclose relation to the BDF methods makes them especially suited for "type-insensitive" codes as well as for being combined with the BDFs in our context .Moreover, although the formulas are of the same order of consistency, the errorconstants of the DCBDFk are less than half those of the BDF(k + 1) .
Table 3 .1 : Error constants of Adams-Moulton (AM) vs . DCBDF methods
Figure 3 .1 : Stability regions for DCBDFk (largest circle : k = 2, smallest : k = 6 .)
The DCBDFs are also one-leg collocation methods, [5], hence for polynomialsP(t) and Q(t),
(3.1a) h-r (pkP)(tn) -P(aktn) = 0;
degP < k
(3.1b)
((YkQ)(t.) - Q(akt .n ) = 0 ;
deg Q < k -1 ; k>2.
STABILIZED MULTISTEP METHODS
9
The point vktn is termed the collocation point . Note that because Ck = 1 -Vk /(k + 1) we have Qktn = to for k > 2, i .e. the collocation point is the sameas for the BDFs. The property (3 .1b) gives us considerable freedom to manipu-late certain terms, in particular the multiplier term . Thus e .g. we have, undersufficient continuity conditions and for k > 2,
ok(GT (x(t))A(t)) GT (x(uktn))A(aktn) + O(h k )GT('kx( tn) + O(h k ))A(cTktn) + O(hk)GT(x(tn))A(rktn) + O(hk)GT (x(tn))Uk'\(tn) + O(hk)GT (x(tn))A(tn) + O(h k ) .
Therefore, instead of using the DCBDFk in (1 .2), which of course fails to satisfythe necessary stability requirements, we may replace the discretization by aDCBDF/BDF compound,
(3 .2a)
h-lpkxn = e7kf (xn) - GT(xn)An(3 .2b)
0 = g(xn),
without violating consistency by more than O(hk ) . The stability of the com-pound DCBDF/BDF follows from Theorem 2 .1 . While the DCBDFk is consis-tent of order p = k + 1 this discretization leads to an order reduction .
We may also view this as a matter of in which order we discretize the systemand impose the constraint, cf. Fig . 3 .2 . Thus, the standard (unstable) discretiza-tion is obtained by first constraining the system x = f (x) using d'Alembert'sprinciple to impose g(x) = 0, followed by discretizing the resulting DAE withDCBDFk. On the other hand, we may choose to first discretize x = f (x) withDCBDFk, and then impose the constraint on the discrete system. The latterstep is not uniquely defined, one needs a "discrete d'Alembert principle" to mo-tivate it thoroughly. Nevertheless, if the constraint is imposed using Lagrangianmultipliers, we obtain the stable discretization (3 .2) . Thus, the diagram in Fig .3.2 is commutative up to order O(hk), but only the lower branch yields a sta-ble discretization . The order of accuracy, however, drops and equals k in allvariables, cf. Remark 2.1 .
In order to stabilize the discretization without loss of accuracy except for inA, we use less obvious modifications than those indicated above and insteadconsider
Again, we only violate consistency by 0(hk), i .e . (3 .3) is O(hk ) away fromboth the pure DCBDFk discretization and the DCBDFk/BDFk compound (3 .2)above. Indeed, if G is constant (3 .3) reduces to (3 .2) . The purpose of the abovemodification is to avoid applying DCBDF to the multiplier term, which is in-stead treated essentially by the corresponding BDF method. For this particular
10
C. AREVALO, C. FUHRER AND G . SODERLIND
x=f
discretizeby DCBDF
h-1Px = of
apply d'Alembert
apply "d'Alembert"
x=f-GT A0=g
discretizeby DCBDF
h 1Px = v(f - GTA)0=g
h- 'px=af - GTA0=g
Figure 3 .2: Different DCBDFk discretizations . The diagram is commutative up toorder O(h k ) . The upper branch is unstable, the lower branch stable .
/3-blocking, Theorem 2 .1 yields a convergence order of p = k + 1 for x andp = k for A. These statements are also verified by our computational results, seeSection 5 .
4 /3-blocked Adams methods .
As the Adams-Moulton (AM) methods are the most accurate and popularamong multistep methods for nonstiff computations, it is desirable to construct/3-blocked AM methods for the given index 2 problem . The stabilizing T =
-V'/(k + 1) appeared naturally for the DCBDFk method and gave optimaldamping for error propagation in the A variables . Taking T = cVk will, however,only partly solve the problem for AM methods . Thus, the AMk method can onlybe stabilized for k < 3, and only for k = 1, when AMI-DCBDF1, is optimaldamping possible (at c = -1/2) . For k = 2, 3 the matrix 13(h) defined in Sec .2, which accounts for the propagation of the error in A, is no longer an O(h)perturbation of a nilpotent matrix like in the DCBDF case, and the spectralradius quickly becomes fairly large .As a result, the AMk methods can only be applied with order p = k + 1 for
x and order p = k for A when k < 3. For larger step numbers, one is forcedto use a rr with more than k steps or an order reduction will occur also in x .The computational results (cf. Sec. 5) confirm the orders predicted by Theorem2 .1. In addition, these results show that the stable 3-blocked AMk methodis more accurate than the corresponding DCBDFk . This is to be expectedbecause of the smaller error constant of the AMk, and observed computationalresults are again in perfect agreement with theory . We have not investigated the
STABILIZED MULTISTEP METHODS
4
1 ~
0.2
11
Figure 4 .1: Maximum root of v-r for AMk methods and 7 = cVk when the parameterc varies. Stabilization is only possible for k < 3 .
practical consequences of the weaker damping of the /3-blocked AM method,but a preliminary test was carried out in [3] . In real computations, the weakerdamping may to some extent reduce the efficiency of the AM methods . Theoptimally damped DCBDF methods may therefore be preferable in practice,and also provide a family of methods of orders up to p = 7 .
5 Computational results and conclusions .
Test problems arising in applications have been solved in order to verify thetheoretical results on convergence . We limit ourselves to assessing the practicalfeasibility of the stabilized methods by solving a sufficiently complex problem ;we simulate the motion of a truck in 2D . This model is described by 20 equationsin index 2 form . For a detailed description of the model and the correspondingequations of motion we refer to [3, 6] .
The methods were used for a sequence of constant stepsizes ranging from 2 -7to 2-15 with a fixed number of Newton iterations . The integration was startedat a steady state of the equations of motion . Errors were computed by com-paring the solutions to a reference solution calculated to a relative accuracy ofapproximately 10-13 . The mean value of the errors (corresponding to measur-ing the error in the L' norm over the range of integration) versus the stepsizeis given in Fig. 4 .2. In full agreement with theory, that diagram shows that thek = 3 step /-blocked DCBDF method and the Adams-Moulton method bothhave order p = k + 1 = 4 for the non-algebraic variables and order p = k = 3 forthe algebraic variables. The difference between these two methods reflects thesmaller error constant of Adams-Moulton. We note that for this test problem,
the /3-blocked methods are approximately 100 times as accurate as the standardBDF3, or equivalently 3-4 times faster for a given accuracy requirement . Forthe algebraic variable, however, differences between the methods are small .One could also compare the DCBDFk to the BDF(k + 1), since these methods
are of the same order . As pointed out earlier, the smaller error constant of theDCBDFk method then implies a 2-4 times higher accuracy than that of BDF(k+1) . The robustness of the DCBDFk methods needs further investigation, butin the nonstiff problems we have tried, no computational difficulties have beenobserved .
We also present effort/precision diagrams of /3-blocked DCBDFk for k = 1 : 5and AMk for k = 1 : 3. AM1 uses optimal damping and is therefore identical toDCBDF1 (or the trapezoidal rule) ; AM2 is /3-blocked with c = -0 .15, cf . Fig .4.1 .
Acknowledgments .
The first author would like to thank the CESMa center of the UniversidadSimon Bolivar for permitting her free use of its research facilities . The secondauthor would like to thank Prof. Ernst Hairer and Prof. Gerhard Wanner en-abling a month's visit to the University of Geneva, where part of this work wasdone. He would also like to acknowledge DAAD for travel support . The thirdauthor gratefully acknowledges partial support by the Swedish Research Councilfor Engineering Sciences TFR under contract no . 222/91-405 as well as travelsupport from the Swedish Institute SI in cooperation with DAAD .
12
10
8
6
2 3 4
12
10
8
6
STABILIZED MULTISTEP METHODS
-log 10(ReLError Position Variable)
DCBDF
Adams
2
10
8
6
4
2
3 4
2 3-log10(Stepsize)
-log 10(Rel .Error Algebraic Variable)
DCBDF
Adams
4
10
8
6
4
2
2 3 4
13
Figure 5 .1 : Effort/Precision diagram for ,0-blocked DCBDF and Adams-Moultonk-step methods
REFERENCES
1. C. Arevalo, Matching the structure of DAEs and multistep methods, Ph. D . thesis,Department of Computer Science, Lund University, Lund, Sweden, 1993 .
2. C. Arevalo and G. Soderlind, Convergence of multistep discretizations of DAEs,BIT-Numerical Mathematics, 35 (1995), pp . 143-168 .
3. C. Arevalo, C . Fuhrer and G . Soderlind, ,0-blocked multistep methods for Euler-Lagrange DAEs. I. Linear analysis, Technical report LU-CS-TR:94-122, Lund Uni-versity, Lund, 1994, submitted for publication .
4. V. Brasey and E. Hairer, Symmetrized half-explicit methods for constrained me-chanical systems, Appl. Num. Mathematics, 13 (1993), pp. 23-31 .
5. G. Dahlquist, On One-leg Multistep Methods, SIAM J. Numer. Anal., 20 (1983),pp. 1130-1138 .
6. E. Eich and C. Fuhrer, Numerical Methods in Multibody Dynamics, Teubner-Verlag,Stuttgart, to appear 1995 .
7. E. Hairer and G. Wanner, Solving Ordinary Differential Equations II. Stiff andDifferential-Algebraic Problems, Springer series in Computational Mathematics,14, Springer-Verlag, New York, 1991 .
9. A. Ostermann, A class of half-explicit Runge-Kutta methods for differential-algebraic systems of index 3, Appl. Numer. Math., 13 (1993), pp . 165-179 .
10. G. Soderlind, A multi-purpose system for the numerical integration of ODEs, Appl .Math. Comp., 31 (1989) , pp. 346-360 .
