GEOMS: A software package for the numerical integration of general model equations of multibody systems

für Angewandte Analysis und Stohastik

im Forshungsverbund Berlin e.V.

Preprint ISSN 0946 8633

GEOMS: A software pakage for the numerial

integration of general model equations

of multibody systems

Andreas Steinbreher

submitted: September25,2007

WeierstrassInstituteforApplied

AnalysisandStohastis

Mohrenstrasse39

10117Berlin

Germany

E-mail: steinbreherwias-berlin.de

No. 1259

Berlin2007

2000MathematisSubjetClassiation. 70E55,65L80.

Keywords andphrases. dierential-algebraiequations, equationsof motion,multibody sys-

tem,numerialintegration,simulation.

This work hasbeensupportedby DFGresearh grantno. Me790/13during the stay of the

authorat TUBerlin,FahgebietNumerisheMathematik. ThesoftwarepakageGEOMS ispart

ofthePhDthesis[38℄oftheauthordevelopedattheTUBerlin.

Weierstraÿ-Institut fürAngewandte Analysisund Stohastik (WIAS)

Mohrenstraÿe 39

10117Berlin

Germany

Fax: +49 302044975

E-Mail: preprintwias-berlin.de

World WideWeb: http://www.wias-berlin.de/

InthispaperwepresentthenewnumerialalgorithmGEOMSforthenumer-

ial integration ofthe mostgeneral formof theequations of motion of multi-

bodysystems,inluding nonholonomi onstraints and possibleredundanies

inthe onstraints, asthey mayappearinindustrialappliations. Besides the

numerialintegrationitoerssomeadditionalfeatureslikestabilizationofthe

model equations, use of dierent deomposition strategies, or heking and

orretionofthe initial valueswithrespettotheir onsisteny. Furthermore,

GEOMSpreserveshiddenonstraints and(possibly)existingsolutioninvariants

iftheyareprovided asequations.

We will also demonstrate the performane and theappliability of GEOMS for

twomehanial examplesof dierent degrees of omplexity.

1 Introdution

The multibody system (MBS) approah is frequently used in industrial simulation

pakages in robotis, vehile system dynamis, and biomehanis. A multibody

system model onsists of a nite number of rigid or elasti bodies and their in-

teronnetions like, e.g., joints, springs, dampers, and atuators. The equations of

motion may be generated in a systemati way by multibody formalisms that are

based onthe priniples of lassial mehanis[35℄.

The eient and robust numerial integration of these equations is a hallenging

problem in the development of simulation pakages, sine dynamial simulation is

frequently used and one of the most time onsuming analysis methods for MBS

models. The equations of motion with nonredundant onstraints form a nonlinear

system of dierential-algebrai equations (DAEs) of dierentiation index (d-index)

3,see [7,10,13,18℄. It iswellknown thatthe numerialtreatmentof DAEs ofhigh

index or higher index, i.e., d-index 2 or larger than 2, respetively, is nontrivial in

general. Eets arising in the numerial treatment are, for example, drift,instabil-

ities, onvergene problems, or inonsistenies. These diulties in the numerial

solutionofsuhhighindexproblemsaredisussedin[4,12,14,15,17,18,24,29,30℄.

However,theequationsofmotionareDAEswithaveryspeialstruturethatshould

be exploitedin the numerial solution[7,18℄.

In this report we will present the new software pakage GEOMS for the numerial

integration ofgeneral equationsof motionof multibodysystems indesriptor form.

Inontrasttostandardtextbookpresentationslike[18℄,wedonotrestritourselves

to lassial onstrained mehanial systems but onsider the more omplex model

equationsthatareatuallyusedinstate-of-the-artMBSsimulationpakages[7,34℄.

namialfore elements, ontat onditions, and (possibly) redundant holonomi as

wellasnonholonomionstraints. Furthermore,thepakage takesintoaountpos-

sibly existing information onerning solution invariants,e.g., energy onservation.

The ode is based on residual evaluations, i.e., the system need not be given om-

pletely in expliit form. It is suient that the right-hand side of the equations of

motionand the mass matrix are speied.

Although, the pakage GEOMS is able to treat also redundant onstraints, in this

paper we will restrit our onsiderations to regular equations of motion, i.e., the

onstraintsare assumed to benonredundant. Formore details on equationsof mo-

tion with redundant onstraints we referto[26, 38℄.

As baseof the integration methodGEOMS wewillpropose aremodelingof the equa-

tionsof motions. Theaimofthis remodelingistodetermineanequivalentformula-

tion, the so alled projeted-strangeness free form, whih has d-index 1but has the

same set of solutions as the original equations of motion. Beause of the redued

d-index, the numerial treatment of the projeted-strangeness free form by use of

impliitODE methodsis not aetedby instabilitiesarising fromthe higherindex.

Furthermore,all(hidden)onstraintsarepreservedsuhthatnodrift-oeetsarise

inthe numerialtreatment. The proposed remodelingan beseen asregularization

of the equations of motion. For more details on the regularization of equations of

motion we refer to [38℄. The integration method implemented in GEOMS ombines

animpliit Runge-Kutta-Methodof order 5 with this regularizationtehnique.

Thereportisorganizedasfollows. InSetion2weintroduetheequationsofmotion

whihwewanttotreatnumeriallyandwedisuss theremodelingtotheprojeted-

strangeness-free formwhih willbe used for the disretizationin GEOMS. In Setion

3 we introdue the ode GEOMS and we disuss its features and its appliability in

detail. InSetion4wedemonstratethepropertiesofthesoftwarepakage GEOMSby

two numerialexamples. Forthe usage and implementationof GEOMSthe manualis

presented in Appendix A.

2 The Equations of Motion and their Remodeling

Here and in the followingwe willuse the following notation.

Notation 2.1 Let

f

^{be a} dierentiable funtion

f : X → R ^m

,

X ⊂ R ⁿ

, and let

x

^{be a} dierentiable funtion

x : I → X

, where

I

is an open interval in

R

. The

i

^{th (total) derivative of}

x(t)

^{with respet to}

t

^{is denoted by}

x ⁽ⁱ⁾ (t) = d ⁱ x(t)/dt ⁱ

^for

i ∈ N ₀

. Note the onvention

x ⁽⁰⁾ (t) = x(t)

^,

x ⁽¹⁾ (t) = ˙ x(t)

^{, and}

x ⁽²⁾ (t) = ¨ x(t)

^{. The}

(partial) derivative of

f(x)

^{with respet to}

x

^{is denoted by}

f ,x (x) = _∂x ^∂ f (x)

^{. The}

same notation is used for dierentiable vetor and matrix funtions. The set of

l

^-times ontinuously dierentiable funtionsfrom

X

to

Y

isdenoted by

C ^l ( X , Y )

^.

⊳

Inthefollowingweinvestigateaspatialmultibodysystemwithholonomiaswellas

nonholonomi onstraints [19, 33℄. More preisely we onsider the following initial

valueproblemonthedomain

I = [t 0 , t f ]

^{onsisting ofthe equationsof motioninthe}

form

˙

p = Z (p)v,

^(1a)

M (p, t) ˙ v = f (p, v, r, w, s, λ, µ, t) − Z ^T (p)G ^T (p, s, t)λ − Z ^T (p)H ^T (p, s, t)µ,

^(1b)

˙

r = b(p, v, r, w, s, λ, µ, t),

⁽¹⁾

0 = d(p, v, r, w, s, λ, µ, t),

^(1d)

0 = c(p, s, t),

^(1e)

0 = H(p, s, t)Z (p)v + h(p, s, t) (= ˘ h(p, v, s, t)),

^(1f)

0 = g(p, s, t)

^{(1g )}

with the initialvalues

p(t 0 ) = p 0 ∈ R ^{n p} , v(t 0 ) = v 0 ∈ R ^{n v} , r(t 0 ) = r 0 ∈ R ^{n r} , w(t 0 ) = w 0 ∈ R ^{n w} ,

s(t 0 ) = s 0 ∈ R ^{n s} , λ(t 0 ) = λ 0 ∈ R ^{n λ} , µ(t 0 ) = µ 0 ∈ R ^{n µ} .

⁽²⁾

Here, the positionvetor

p

^{ontainsarbitraryposition oordinatesof the multibody}

system. The Euler-Lagrange formalism for modeling multibody systems yields the

equations of motion in seond order form. In order to transform the seond order

system to an equivalent rst order system we introdue a veloity vetor

v

^{and get}

the relation (1a) between the generalized veloities

p ˙

^{and the veloities}

v

^{with a}

matrix

Z (p)

^{, that determines the angular veloities. The equations (1a) are alled}

kinemati equations. The transformation matrix

Z (p)

^{ours only if there are ro-}

tations in three dimensional spae, it may be determined by Poisson's kinematial

equations[1, 7℄. In the twodimensional ase we have

Z(p) = I, p ˙ = v

^.

The equations (1b) are alled dynami equations of motion. They follow from the

equilibrium of fores and momenta and inlude the mass matrix

M (p)

^{, the vetor}

of the applied and gyrosopi fores

f(p, v, r, w, s, λ, µ, t)

^{, the onstraint matries}

G(p, s, t)

^and

H(p, s, t)

^{oftheholonomiand}nonholonomionstraints,respetively, whihontainthe inaessiblediretions ofmotionolumn-wise,the assoiatedon-

straint fores

G ^T (p, s, t)λ

^and

H ^T (p, s, t)µ

^{, and the Lagrange} multipliers

λ

^and

µ

^.

The holonomi onstraint matrix is dened as

G(p, s, t) = _dp ^d g(p, s(p, t), t)

^{. The}

massmatrix

M (p)

^ispositivesemi-denite,sine thekinetienergyisanonnegative quadratiform, and inludesthe inertia propertiesof the multibodysystem.

In a real multibody system, there are often dynami fore elements whih are de-

sribed by the vetor

r

^{and determinedby equations(1), see [7℄.}

Furthermore, not all onstraints of a multibody system are diretly desribed by

the position variables

p

^{or the veloity variables}

v

^{, but depend on ertain ontat}

points with oordinates

s

^{on the surfae of some bodies. The} relationship between these ontat point oordinates

s

^{and the position variables}

p

^{are given by (1e).}

Furthermore, the equations of motion are aeted by the

n λ

^{holonomi onstraints}

(1g) and

n µ

nonholonomi onstraints (1f). These onstraints are also alled the holonomi onstraints on positionlevel and the nonholonomi onstraints on velo-

ity level, respetively. Sometimes, fore laws and onstraints may be formulated

more onveniently using auxiliaryvariables

w

^{that are impliitly dened by the}

n w

Here,

n = n p + n v +n r + n w +n s + n λ +n µ

^{denotes thenumberofunknown variables.}

Furthermore,many motionsof mehanial systems have known solutioninvariants,

i.e., relations whih are satised along any motion of the mehanial system, like

the invariane of the total energy, momentum, or impulse. Let us denote the

m e

equationsdesribing suh solutioninvariantsby

0 = e(p, v, s, t).

⁽³⁾

In partiular, onservative multibody systems are energy onserving. In this ase

the total energy is onstant along every motion of the system. For more details on

solutioninvariants we referto[38℄.

The theoretial basis of the ode GEOMS is based onthe following assumptions.

Assumption 2.2 Consider the equations of motion (1). Then the matries

a) d ,w ,

^(4a)

b) c ,s ,

^(4b)

c) M

⁽⁴⁾

d)

GZM ⁻¹ G λ GZM ⁻¹ H µ

HZM ⁻¹ G λ HZM ⁻¹ H µ

(4d)

are assumed to be nonsingular with a bounded inverse for all

(p, v, r, w, s, λ, µ, t) ∈ M

, see (10), where

G λ = Z ^T G ^T − f ,λ + f ,w d ⁻¹ _,w d ,λ ,

⁽⁵⁾

H µ = Z ^T H ^T − f ,µ + f ,w d ⁻¹ _,w d ,µ .

⁽⁶⁾

Furthermore,it is assumed that

d ∈ C ¹ ( M , R ^{n w} ), c ∈ C ¹ ( M , R ^{n s} ), ˘ h ∈ C ² ( M , R ^{n µ} ), g ∈ C ³ ( M , R ^{n λ} ).

Remark 2.3 a) The nonsingularity of the mass matrix

M

^{is assumed only for}

reasons of simpliity. It is not neessary for the suessful numerial integration

with GEOMS.

b) Furthermore, note that in Assumption 2.2 redundant onstraints are exluded.

Redundant onstraints may result in a nonuniqueness of the Lagrange multipliers.

Nevertheless, GEOMS isable todeal withertaintypesof redundant onstraints. For

more detailson redundant onstraintssee [26, 38℄.

⊳

Usingthe equationsofmotion (1),the rst and seondderivativeswith respet to

t

of the holonomi onstraints(1g) are given by

0 = g ^I (p, v, s, t) = d

dt g(p, s, t)

^(7a)

= GZv + g ,t − g ,s c ⁻¹ _,s c ,t

^(7b)

0 = g ^II (p, v, r, w, s, λ, µ, t) = d ²

dt ² g(p, s, t)

^(8a)

= (g _,p ^I − g _,s ^I c ⁻ _,s ¹ c ,p )Zv + GZM ^{− 1} (f − Z ^T G ^T λ − Z ^T H ^T µ) + g _,t ^I − g _,s ^I c ⁻ _,s ¹ c ,t .

^(8b)

They are alled holonomi onstraints on veloity level (7a) and holonomi on-

straints on aeleration level (8a), respetively. The rst derivative with respet to

t

^{of the} nonholonomionstraints (1f)is given by

0 = h ^I (p, v, r, w, s, λ, µ, t) = d

dt (H(p, s, t)Z(p)v + h(p, s, t))

^(9a)

= (˘ h ,p − ˘ h ,s c ⁻ _,s ¹ c ,p )Zv + HZM ^{− 1} (f − Z ^T G ^T λ − Z ^T H ^T µ) + (˘ h ,t − h ˘ ,s c ⁻ _,s ¹ c ,t )

^(9b)

whih are alled nonholonomi onstraints on aeleration level (9a).

The holonomi onstraints on veloity level and on aeleration level in form (7b)

and (8b),respetively, aswellas thenonholonomionstraintsonaelerationlevel

in form (9b) turn out to be the hidden onstraints of the equations of motion, see

[38℄. Thehoieofvalues

(p, v, r, w, s, λ, µ, t) ∈ R ⁿ × I

isrestritedby allonstraints inludingthe hidden onstraints,i.e., (1d)-(1g), (7b), (8b), and (9b). Values whih

satisfyalloftheseonstraintsarealledonsistent andwegetthe setof onsisteny

M = { (p, v, r, w, s, λ, µ, t) ∈ R ⁿ × I : 0 = d(p, v, r, w, s, λ, µ, t),

⁽¹⁰⁾

0 = c(p, s, t),

0 = H(p, s, t)Z (p)v + h(p, s, t), 0 = g(p, s, t),

0 = h ^I (p, v, r, w, s, λ, µ, t), 0 = g ^I (p, v, s, t),

0 = g ^II (p, v, r, w, s, λ, µ, t) } .

Theorem 2.4 Letthe equationsof motion (1) satisfy Assumptions 2.2. Thenthere

exist matrix funtions

S p ∈ C ⁰ ( M , R ^{n fp ,n p} )

^and

S v ∈ C ⁰ ( M , R ⁿ fv ,n v )

^with

n f p = n p − n λ

^and

n f v = n v − n λ − n µ

^{suh that the matrix funtions}

S p (p, t) G(p, t)

and





S v (p, t)M(p, t) G(p, t)Z (p) H(p, t)Z (p)





^are nonsingular (11)

for all

(p, v, r, w, s, λ, µ, t) ∈ M

. Then the dierential-algebraisystem

S p (p, t) ˙ p = S p (p, t)Z(p)v,

^(12a)

S v (p, t)M(p, t) ˙ v = S v (p, t)f(p, v, r, w, s, λ, µ, t)

^(12b)

− S v (p, t)Z ^T (p)G ^T (p, s, t)λ − S v (p, t)Z ^T (p)H ^T (p, s, t)µ,

˙

r = b(p, v, r, w, s, λ, µ, t),

⁽¹²⁾

0 = d(p, v, r, w, s, λ, µ, t),

^(12d)

0 = c(p, s, t),

^(12e)

0 = H(p, s, t)Z(p)v + h(p, s, t),

^(12f)

0 = g(p, s, t),

^{(12g )}

0 = h ^I (p, v, r, w, s, λ, µ, t),

^(12h)

0 = g ^I (p, v, s, t),

⁽¹²ⁱ⁾

0 = g ^II (p, v, r, w, s, λ, µ, t)

^(12j)

has d-index 1 and the same set of solutions as the equations of motion (1).

Proof. The proof an befound in[38℄.

Remark 2.5 a)The matrix funtions

S p

^and

S v

^{are alled kinemati seletor and}

dynami seletor, respetively.

b)WewillalltheDAE (12)the projeted-strangeness-freeformulation of the equa-

tions of motion. In [21, 22, 23, 24℄ the strangeness-onept is introdued as tool

for the lassiation of general nonlinear DAEs inluding over- and underdeter-

minedDAEs. In partiular, so alled strangeness-free DAEs are introdued. Apart

from the over- or underdeterminedness strangeness-free DAEs behave like DAEs

with d-index 1 while nonstrangeness-free DAEs behave like DAEs with d-index 2

or larger. Strangeness-free DAEs do not ontain hidden onstraints. In partiular,

in [21, 22, 23, 24℄ it is pointed out that strangeness-free DAEs and, therefore, the

projeted-strangeness-free formulation of the equations of motion (12), are suited

and preferable for the numerial treatment using sti ODE solvers like impliit

Runge-Kutta-Methods orBDF methods.

) The algorithm GEOMS is based on a projeted-strangeness-free from (12) of the

equationsofmotionbutitisnotneessarythatthisformisprovidedbytheuser,i.e.,

theuserdoesnothavetoperformtheregularizationtotheprojeted-strangeness-free

form. It is suient, if the user provides the onstraints onveloity level (7b) and

on aeleration level (8b) and (9b) in addition to the original equations of motion

(1)and, if available, (3). By use of so alled order-n-formalisms for the evaluation

ofthe equationsofmotiontheonstraintsonveloityleveland onaelerationlevel

are omputed automatially,see [8, 34℄.

⊳

Withthese preparationswehave presented all the toolsto perform the onsisteny

preserving index redution of the equationsof motion (1)as follows.

The equationsof motion (1) are assumed tosatisfy Assumptions 2.2. Furthermore,

let

M ∈ C ⁰ ( M _p , R ^{n v ,n v} )

^and

Z ∈ C ⁰ ( M _p , R ^{n p ,n v} )

^{, where}

M _p = M ∩ ( R ^{n p} × I )

^{is the}

set of onsistent

(p, t)

^.

Thentheregularizationviaonsisteny preservingindexredutionisdonebyhoos-

ingaseletor

S p ∈ C ⁰ ( M _p , R ^{n fp ,n p} )

^andaseletor

S v ∈ C ⁰ ( M _p , R ⁿ fv ,n v )

^dependingon

(p, u)

^with

n f p = n p − n λ

^and

n f v = n v − n λ − n µ

^{, in the following way.}

1. Determinationof seletor

S p

(a) Determine

K p ∈ C ⁰ ( M _p , R ^{n p ,n fp} )

^dependingon

(p, t)

^{suhthattheolumns}

of

K p (p, t)

^span

ker(G(p, s(p, t), t))

^{for all}

(p, t) ∈ M _p

.

(b) Determinetheseletor

S p ∈ C ⁰ ( M _p , R ^{n fp ,n p} )

^dependingon

(p, t)

^suhthat

S p (p, t)K p (p, t)

^isnonsingular for all

(p, t) ∈ M _p

.

2. Determinationof seletor

S v

(a) Determine

K v ∈ C ⁰ ( M _p , R ^{n v ,n} fv )

^dependingon

(p, t)

^{suhthattheolumns}

of

K v (p, t)

^span

ker(

G(p, s(p, t), t)Z(p) H(p, s(p, t), t)Z(p)

)

for all

(p, t) ∈ M _p

.

(b) Determinetheseletor

S v ∈ C ⁰ ( M _p , R ⁿ fv ,n v )

^dependingon

(p, t)

^suhthat

S v (p, t)M (p, t)K v (p, t)

^isnonsingular for all

(p, t) ∈ M _p

.

3. Projeted strangeness-free form of the equations of motion

By appending the onstraints on veloity level (7b) and the onstraints on

aeleration level (8b) and (9b), the projeted-strangeness-free form of the

equationsof motion is given by (12).

Withthis algorithmweare abletodetermine aprojeted-strangeness-freeform(12)

of the equations of motion whih ontains all informationof the set of onsisteny

(10). The projeted-strangeness-free form (12) that is reated in this way is ana-

lytiallyequivalent to the originalequations of motion in the sense that both have

the same set of solutions. Therefore and beause of Remark 2.5, the projeted-

strangeness-free form(12) an be seen asa regularization tehnique. In partiular,

the semi-impliitformofthe projeted-strangeness-free form(12)isof greatadvan-

tage,sine allonstraintsare stated aspurely algebraiequations,and thereare no

redundanies among the algebraionstraints and the dierential equations.

Remark 2.7 Notethat Seletors

S p

^and

S v

^{satisfyingthe rankonditions (11) are}

not uniquelydetermined. Ratheritispossibletohoose theseletors inapieewise

onstant fashion. In priniple, the seletors may be kept onstant as long as the

Newton iteration matrix

N

(see Page 15) remains nonsingular. But the hoie of

tion. Therefore, with respet to the onditioning of the linear systems whih have

tobe solved duringthe Newtoniteration,the seletors should bereomputed early

enough and not just shortly before reahing a state, where the Newton iteration

matrix beomessingular. This fatis treatedinGEOMS by the reomputation ofthe

seletors if theolumn pivotingwith respet tothe algebraionstraintshangesor

onvergene problems of the Newton iteration our. This is demonstrated in two

simulationsenarios whih are depited in Tables 3 and 4.

Note that the pieewise onstant hoie of the seletors is of great advantage and

importanefor the numerialintegration, beause itoers the possibilityto redue

the amount of omputational work for the omputation of the seletors. In par-

tiular, this means, that the ondition number of the Newton iteration matrix

N

depends diretly onthe hoie of the seletors.

⊳

Example 2.8 The mathematial pendulum: Let us onsider a mathematial

pendulum, of length

L > 0

^{whih represents a point mass moving withoutfrition}

alonga vertial irle ofradius

L

^{undergravitydenoted by the gravity aeleration}

g

^{. Forthe desriptionofthe ongurationofthe pendulumwe hoose Cartesiano-}

ordinates

p =

x y T

denotingthe positionofthemass

m

^inthetwodimensional spae

R ²

. The equations of motion of rst order have the form

p ˙ 1

˙ p 2

= v 1

v 2

,

^(13a)

m 0 0 m

˙ v 1

˙ v 2

=

0

− mg

− 2p 1

2p 2

λ 1

,

^(13b)

0 =

p ² ₁ + p ² ₂ − L ²

.

⁽¹³⁾

The holonomi onstraintson veloity leveland onaelerationlevelare given by

0 =

2p 1 v 1 + 2p 2 v 2

,

^(13d)

0 =

2v ₁ ² + 2v ² ₂ − 2p 2 g − m ⁴ (p ² ₁ + p ² ₂ )λ 1

,

^(13e)

respetively. Following Algorithm 2.6 we have to onsider

G =

2p ₁ 2p ₂ .

^The

matrix funtion

K p

^{an be determinedas}

K p = − p 2

p 1

and, therefore,the seletor

S p

^{an be hosen as}

S p =

− p 2 p 1

suhthat

S p K p =

− p ₂ p ₁ − p 2

p 1

=

p ² ₂ + p ² ₁

= L ²

,

see the onstraints (13). Sine the mass matrix is given by

M = mI

^{, we an use}

S v = S p

^{and weget the} projeted-strangeness-free formulation

− p 2 p ˙ 1 + p 1 p ˙ 2 = − p 2 v 1 + p 1 v 2 ,

^(14a)

− mp 2 v ˙ 1 + mp 1 v ˙ 2 = − mgp 1 ,

^(14b)

0 = p ² ₁ + p ² ₂ − L ² ,

⁽¹⁴⁾

0 = 2p ₁ v ₁ + 2p ₂ v ₂ ,

^(14d)

0 = 2v ² ₁ + 2v ₂ ² − 2p 2 g − 4

m (p ² ₁ + p ² ₂ )λ 1 .

^(14e)

AsmentionedinRemark2.7,theseletors

S p

^and

S v

^{arenotuniquelydeterminedby}

the onditions (11) or the Algorithm2.6. In partiular, the seletors an behosen

tobe pieewise onstant.

Letusonsiderthisfatforthependulumwiththe initialstate

p 1 = 0

^and

p 2 = − L

^,

i.e., the pendulum is hanging downwards. In this position the seletors an be

determined as

S p (p, u) = S v (p, u) =

L 0

.

⁽¹⁵⁾

Keeping these seletors onstant, the leading matrix of the left-hand side of the

underlying ordinary dierential equations, (obtained by substituting the algebrai

equationsin (14) by their derivatives with respet to

t

^)is













L 0 0 0 0

0 0 mL 0 0

2p ₁ 2p ₂ 0 0 0

× × 2p 1 2p 2 0

× × × × m ⁴ (p ² ₁ + p ² ₂ )













.

⁽¹⁶⁾

Obviously, therank onditions (11)are fullled and the leadingmatrix (16) isnon-

singular,aslong as

p 2

^{doesnot beome zero. In partiular,this means that aslong}

as the pendulum does not reah one of the horizontal positions, i.e.,

p ₁ = ± L

^and

p 2 = 0

^{, the seletors an be hosen onstant asin(15). Otherwise, if the pendulum}

reahes orpasses the horizontalposition,the matrix (16) beomes singularand the

rst and third as well as the seond and fourth equations are redundant suh that

the solution is not uniquely dened. Furthermore, the ondition number of matrix

(16) goes toinnity as

p 2

^{goesto zero.}

For these reasons, in the neighborhood of the horizontal position of the pendulum

newseletorshavetobedetermined. SeealsotheExample4.1fornumerialresults.

⊳

3 GEOMS

The ode GEOMS is implemented in FORTRAN77 and furthermore, there exists a

MATLAB [20℄ interfae via MEX les for the diret usage of GEOMS in MATLAB.

GEOMS.

In GEOMS the 3-stage impliit Runge-Kutta Method Radau IIa of order 5, see [18℄,

asdisretization ofthe projeted-strangeness-free formulation(12) of the equations

of motion is implemented. Although, GEOMS bases on the presented stabilization

tehnique developed in [38℄ and presented in Theorem 2.4, i.e., GEOMS uses the

projeted-strangeness-free formulation(12) for the disretization, the user doesnot

have to provide the projeted-strangeness-free formulation. Instead the user has

to provide all neessary information, i.e., in partiular, the hidden onstraints in

additiontothe originalequations of motion(1) and ,if available,(3).

TheRunge-Kuttamatrix

A

^{,the weightvetor}

b

^{,and the nodevetor}

c

^{are given by}

the Buther tableau

c A

b ^T ⇔

4− √ 6 10

88−7 √ 6 360

296−169 √ 6 1800

−2+3 √ 6 225 4+ √

6 10

296+169 √ 6 1800

88+7 √ 6 360

−2−3 √ 6 225

1 ¹⁶⁻ ₃₆ ^{√ 6 16+} ₃₆ ^{√ 6 1} ₉

16 − √ 6 36

16+ √ 6 36

1 9

,

⁽¹⁷⁾

see [17, 18℄. The algorithmGEOMS is designed tohandle equations of motion of the

form(1)withpossibleredundantonstraintsaswellaswithpossibly knownsolution

invariants(3)whihmaybeprovidedasadditionalequations. Ifthemass matrix

M

is nonsingular and the onstraints are nonredundant then the equations of motion

havetosatisfyAssumption2.2. Ifthisisnottheasesomefurtherrankassumptions

have to be satised. Formore details see [38℄.

Here andinthe followingwewilluse thetypewriterstyleforobjetswhiharepart

of the soure odes of the implemented numerial algorithms. In partiular, this

involves names of subroutines like GEOMS, GEERREST, IVCOND, and variables like T,

X, NWTMAT, CALSEL.

In the following we will disuss the features of GEOMS in detail. For the use and

implementationof GEOMS see the manual inAppendix A.

The information of the equations of motion needed from the integration algorithm

has to beprovided inthe following form.

The vetor of unknown variables has to be inthe form

x ^T = X ^T =

w ^T λ ^T µ ^T r ^T v ^T s ^T p ^T

and the right-hand side in (1) and (3)of the hidden onstraints has to be speied

in a user-supplied subroutine with a name given by the user. The dierent parts

RDA =







































d(p, v, r, w, s, λ, µ, t) g ^II (p, v, r, w, s, λ, µ, t) h ^I (p, v, r, w, s, λ, µ, t)

g ^I (p, v, s, t)

H(p, s, t)Z(p)v + h(p, s, t) e(p, v, s, t)

c(p, s, t) g(p, s, t) b(p, v, r, w, s, λ, µ, t)

f(p, v, r, w, s, λ, µ, t) − Z ^T (p)G ^T (p, s, t)λ − Z ^T (p)H ^T (p, s, t)µ Z(p)v











































 (a)





 (b)

(c) (d) (e) (f)



 

 

 

 

 



 

 

 

 

 

 (A)





 (D)

(18)

In partiular, the right-hand side has to be ordered suh that the algebrai part,

i.e., the upper part (18A), ontains the algebrai onstraints ordered with respet

totheir dependenies, i.e., rst (18a) , the onstraints whih restritthe additional

variables

w

^{aswellastheLagrange}multipliers

λ

^and

µ

^{,seond(18b),theonstraints}

on veloity level and the information onerning solution invariants whih restrit

theveloities

v

^{,and third(18),the onstraintsonpositionlevel,whihrestritthe}

position

p

^{andthe ontatvariables}

s

^{. ThespeiedorderleadstoaJaobianofthe}

algebraipart withrespet to

x

^{whih has already blok upper triangularstruture}

that willbe exploitedin GEOMS.

The dierential part, i.e., the seond part (18D), ontains the right-hand side of

the dierential equations also ordered in the same way as the algebrai part. We

rst (18d) have the equations that desribe the behavior of the dynamial fore

elementsfollowedby (18e)the dynamialequationsof motionand,nally,(18f)the

kinematialequations of motion.

In some ases the onstraints of aeleration level (8) and (9), i.e.,

0 = g ^II

^and

0 = h ^I

^{, are not expliitly available or diult to evaluate. In this ase GEOMS is}

alsoappliable. Butone shouldnote thatonlyifallalgebraiinformation,inluding

0 = g ^I

^,

0 = g ^II

^{, and}

0 = h ^I

^{are provided,} instabilitiesand drift an be avoided by GEOMS. It ispreferable toprovideas muh informationaspossible. In the ase that

theonstraintsonaelerationlevelare missing,the providedinformationissimilar

toa DAE that behaves likea DAE with d-index 2.

This fat has to be ommuniated by the user to the ode GEOMS with help of the

optionIOPT(5)=FORM. If IOPT(5)=0 then the projeted-strangeness-free form(12)

of the equations of motion will be expeted as basis for the disretization. Thus,

the user has tospeify all informationof the hidden onstraints, i.e., up toaeler-

ation level. If IOPT(5)=1, then the disretization will be done without speifying

the onstraints on aeleration level

0 = g ^II

^and

0 = h ^I

^{. In the latter ase the}

used formulation of the equations of motion behaves like a system of d-index 2,

i.e., it is not strangeness-free. Beause of the fatthat the used formulation is not

strangeness-free, the suess of the numerialintegration depends highlysensitively

on the problem and on the onsisteny of the initial values, in partiular, on the

preserving invariantsolutions 4

preserving hidden onstraints 5

preserving nonholonomionstraints 2

taking into aount of redundanies in the on-

straints

19

IOPT( 2) LUN optionaloutput for integration information

IOPT( 3) NIT maximal numberof Newton iterations 17

IOPT( 4) STARTN startingvaluesfortheinternalstagesintheNew-

ton iteration

15

IOPT( 5) FORM inomplete regularization 11

IOPT( 6) NMAX maximal numberof integration steps 19

IOPT( 8) PRED step size ontrol 18

IOPT( 9) NWTMAT approximationoftheNewtonmatrixat

x 0

^orone

of the extrapolated stages possible

15

IOPT(10) NWTUPD update of the Newton matrix 17

IOPT(11) DECOMPC LU, QR, or SV deomposition for the algebrai

part

17

IOPT(12) DECOMPD LUorQRdeompositionforthedierentialpart 17

IOPT(13) SELCOMP seletorontrol 18

IOPT(14) AUTONOM exploitationof autonomous equationsof motion 19

IOPT(15) MASSTRCT exploitationof the struture of the mass matrix 19

IOPT(17) IVCNSST hek and orretion of the initial values with

respet toitsonsisteny

14

Table 1: Options and features of GEOMS

onsisteny of the Lagrange multipliers

λ

^and

µ

^.

An overview over the features of GEOMS is given in Table 1. Furthermore, in Table

2 the subroutines belongingto GEOMSand theirtask are listed.

The initial values are of great importane for the existene and the uniqueness

of the solution. For the existene of a ontinuous solution the onsisteny of the

initialvalues is neessary. In partiular, admissible initial values are restrited by

the (hidden) onstraints. On the otherhand onsistent initialvalues, in partiular,

onsistentinitialLagrangemultipliers,are notautomatially given by the modeling

proess and their determination by solving a system of nonlinear algebrai equa-

tions isdiultfor omplexmultibodysystems with alarge numberof onstraints.

Therefore, the algorithmGEOMSprovides the possibilitytodetermine onsistentini-

tialvalues.

In addition to the algebrai equations determining the set of onsisteny

M

, see (10), the user has to dene in a subroutine IVCOND additional onditions to deter-

mine onsistent initial values. Suh onditions oer the possibility to determine

some of the freely hoosable variables or to give further relations whih allows a

unique determinationof onsistent initialvalues.

GEBSUBST bakward substitution of the algebraipart

GECORE ore routine

GEDECCLU deomposition of the algebraipart with LU deomposition

GEDECCQR deomposition of the algebraipart with QRdeomposition

GEDECCSV deomposition of the algebraipart with SV deomposition

GEDECDLU LU deomposition of the dierentialpart

GEELIMFXQ eliminationinthedierentialpart aordingtoQRdeomposition

of the algebraipart

GEELIMFXS eliminationinthe dierentialpart aording toSV deomposition

of the algebraipart

GEELIMMIQ eliminationin the mass matrix and the identity of the kinemati-

al equations of motion aording to QR deomposition of the

algebraipart

GEELIMMIS elimination in the mass matrix and the identity of the kinemat-

ial equations of motion aording to SV deomposition of the

algebraipart

GEERREST error estimation, see Page 18

GEFXNUM numerialapproximation ofthe Jaobianof the right-handside of

the equationsof motion

GEGREPEQ piking relevant olumnsof the dierential part aording to QR

deomposition

GEGREPES piking relevant olumns of the dierential part aording to LU

and SV deomposition

GEINIVAL determinationof onsistent initialvalues, see Page 14

GEOMS main routine

GESOLDLU solving the dierentialpart by use of LU deomposition

GESOLDQR solving the dierentialpart by use of QRdeomposition

GETRFRHSC transformation of the right-hand side aording to the algebrai

part

User-supplied subroutines

EOM providesthe redued derivativearray RDA (18)

IVCOND providesadditionalinitialonditionsneeded forthe onsistentini-

tialization,see Page 12

JAC providesthe Jaobianof the redued derivative array

MAS providesthe mass matrix

SOLOUT outputofthenumerialsolutionandadditionalinformationduring

integration

Table 2: Subroutines of GEOMS

Example 3.1 The mathematial pendulum: In Example 2.8 we have intro-

dued the mathematial pendulum. The position variables

p

^{are restrited to the}

irlewithradius

L

^{,i.e., theonstraintonpositionlevelisgivenby}

0 = p ² ₁ + p ² ₂ − L ²

^.

If one of the positionvariablesis given, the other is uniquely determined up tothe

By dening additionalonditions via the subroutine IVCOND the user an forethe

pendulum into a deviation of

π/4

^{by setting}

p 1 = L/ √

2

^{or by}

0 = p 1 + p 2

^{, for}

instane. Furthermore, a ertain angular veloity

ω

^{an be presribed by}

0 = p v ₁ ² + v ₂ ² /L − ω

^.

⊳

The determination of onsistent initial values is done in the subroutine GEINIVAL

andis basedontheolletionofallalgebraionstraints(1d)-(1g)and (7), (8),and

(9)together with the onditions dened inthe subroutine IVCOND.

Theuserhastodeideifthegiven initialvaluesareassumed tobeonsistentornot.

By setting IOPT(17)=IVCNSST=1, the initial values are assumed to be onsistent

andnohek ofonsisteny ororretionof theinitialvaluesisdone duringtherun

ofGEOMS. Notethat nononsistentinitialvaluesould leadto onvergene problems

in the integration proess whih leads to an abortof the run of GEOMS. Otherwise,

by settingIOPT(17)=0, theinitialvaluesareonsidered tobepossiblyinonsistent.

Thus, onsisteny willbeheked andthe initialvalueswillbe orretedduringthe

run of GEOMS,if neessary. Ifthe user does not providesuiently many additional

onditions, onlythe onsisteny is heked. Ifthe initialvalues are onsistent,then

the integration will be ontinued, otherwise the run of GEOMS will be stopped. If

the user providesmore additionalonditions than neessary, then the orretion (if

neessary) is done regarding the overdetermined nonlinear system. Ifall onditions

together are nonontraditory, then onsistent initial values will be determined.

Otherwise, the Newton iteration used in this proess will diverge and the run of

GEOMS willbe stopped.

ThesolutionofthenonlinearsystemofequationsisobtainedviaasimpliedNewton

method with the possibilityof a ertainnumberof updates of the iteration matrix,

asdesribed atPage17. Thestoppingriterionisthesameasthatforthesimplied

Newton methodduring the integration proess desribed at Page 16.

Remark 3.2 Note the fat that the onditions provided to IVCOND by the user

dominate the given initialguess, i.e., if the given initialguessis onsistent but does

not satisfy the (possibly wrong) onditions provided by IVCOND, then the initial

guess will be orreted in suh a way that both, the onstraints (1d)-(1g) and (7),

(8), and (9)and the initialonditionsprovided toIVCOND are satised.

Inaseofaninitialguesswhihisonsistenttotheonstraints,theoptionIOPT(17)

an be set to one toavoid suha orretion. Otherwise, the onditions provided to

IVCOND should be adapted.

⊳

Ifthere is only interest inthe omputation of onsistent initialvalues, the user has

toset T=TENDand IOPT(17)=0. Thenthe ode GEOMS determinesonsistent initial

values, willall the user-suppliedsubroutine SOLOUT, and nally willreturn to the

allingsubroutine.

In the following we willdisuss the approahwhih is used in the algorithmGEOMS

forthe numerialintegrationof the equationsof motion(1)and, ifavailable,(3)by

use ofthe three stageRunge-Kuttamethodof type RadauIIaof order 5. Let

s = 3

denotethe numberof stages.

As mentioned above, the ode GEOMS ombines the disretization method with the

regularization tehnique presented in Theorem 2.4. Therefore, the algorithm uses

the projeted-strangeness-free form (12) as basis for the disretization. For more

detailsonthedisretizationwereferto[38℄. This disretizationleads toanonlinear

stage equation for the determination of the three stages

X ki ∈ R ⁿ

,

i = 1, 2, 3

^on

the urrent integration interval

[t k , t _k+1 ]

^with

t _k+1 = t k + h k

^{. Here}

h k

^{denotes the}

urrent step size. The stages

X ki ∈ R ⁿ

,

i = 1, 2, 3

approximate the solution at the points

t ki = t k + c i h k

^{The nonlinear stage equation has to be solved by use of a}

(simplied) Newtonmethod.

A good hoie of starting values

X _ki ⁰

^,

i = 1, 2, 3

^{is very important for the on-}

vergene of the Newton iteration. In the ode GEOMS two dierent possibilities for

the determinationof starting values for the integration step from

t k

^to

t k+1

^{are im-}

plemented. The user has to dene in advane whih of both shall be used during

the integration proess.

BysettingIOPT(4)=STARTN=1thestartingvaluesfortheinternalstagesarehosen

by

X _ki ⁰ = x k

^,

i = 1, 2, 3

^{, where}

x k

^{denotes the already known value whih approx-}

imates the solution at the point

t k

^{. This}

x k

^{orresponds either to the initialvalue}

inthe rstintegrationstep, i.e.,

k = 0

^{,oritorresponds tothe value determinedat}

the end of the preedingintegration step.

On the other hand setting IOPT(4)=0 (whih is the default) the starting values

X _ki ⁰

^,

i = 1, 2, 3

^{forthe Newtoniterationareobtained byevaluatingthe}interpolation polynomial

q(t)

^ofdegree

s

^{overthealreadypassed}integrationinterval

[t k −1 , t k ]

^with

t k −1 = t k − h k −1

^andwith

q(t k −1 ) = x k −1 , q(t k −1 +c i h k −1 ) = X k −1i , i = 1, 2, 3.

^Inthis

wayweobtainthestartingvaluesfortheNewtoniterationas

X _ki ⁰ = q(t k +c i h k ), i = 1, 2, 3,

^where

x k −1

^{denotes the numerial solution at the point}

t k −1

^{. In partiular,}

this means that the new starting values in the integration step from

t k

^to

t _k+1

^are

obtained by extrapolation to the points

t k + c i h k

^,

i = 1, 2, 3

^{based on the internal}

stages of the earlierintegration step from

t k −1

^to

t k

^{. Of ourse, this is not possible}

inthe rst step. Formore details see [18℄.

In GEOMS asimpliedNewton methodis implemented. Formore detailson Newton

methodswerefer to[6℄. In partiular, this meansthat a onstant Newtoniteration

matrix

N

is used during the whole or several parts of the Newton iteration inside the urrent integration step

[t k , t k+1 ]

^{. We use the simpliedNewton method,sine}

a onstant Newton iteration matrix redues the amount of omputationbeauseof

thesaved evaluationofJaobiansandsaved deompositionsoftheNewtoniteration

matrix inevery exept the rst Newton iterationstep. But the partiular hoie of

theNewtoniterationmatrixinuenestheonvergene ofthe Newtoniteration. For

thisreason,the odeGEOMS oersthe possibilitytohoose between several referene

points

(X ^∗ , t ^∗ )

^forthedeterminationoftheNewtonmatrix. Thehoiehastobede- termined by the user by settingthe optionIOPT(9)=NWTMAT. The range ofpossible

hoies isrelatedtothe stages duringthe integration step. As disussed previously,

thedeterminationofthe internalstages. In aseof IOPT(4)=0the initialvaluesare

obtained by extrapolation of the solution omputed so far in the points

t k + c i h k

^,

i = 1, 2, 3

^{. This oers the} possibility to approximate the Newton iteration matrix at four dierent referene points

(X ^∗ , t ^∗ ) = (X _ki ⁰ , t k + c i h k )

^for

i = 0, ..., 3

^{, where}

c 0 = 0

^and

c i

^,

i = 1, 2, 3

^{orrespond tothenode vetor ofthe} Runge-Kuttamethod,

see Table 17. Furthermore,

X _ki ⁰

^{orresponds tothe} extrapolated startingvalues for the internal stages at the times

t k + c i h k

^,

i = 0, ..., 3

^{, and, in partiular,}

X _k0 ⁰ = x k

orresponds to the initial state of the urrent integration interval. Note that this

possibility is only given if the initial values for the Newton iteration are extrapo-

lated, i.e., if IOPT(4)=0. In the ase of initial values hosen suh that

X _ki ⁰ = x k

for all

i = 1, 2, 3

^this possibility is not given and the Newton iteration matrix will be approximated at the initial point

(x k , t k )

^{with the initial state of the urrent}

integration step

[t k , t _k+1 ]

^.

Several numerial experiments have shown that the onvergene of the Newton it-

eration an be improved by use of extrapolated initial values, i.e., IOPT(4)=0 in

onnetion with an approximation of the Newton iteration matrix at the seond

internalstage,i.e.,

(X ^∗ , U ^∗ ) = (X _k2 ⁰ , t k + c 2 h k )

^{withIOPT(9)=2. But, ifthe Newton}

iterationdetets onvergene problems,and the integration step has tobe repeated

with a smaller step size, then the Newton iteration matrix has to be reomputed

suh that the overall omputationtime may inrease if the number of times aon-

vergene problems is deteted is large. This number is reeted in the ounter

NCRJCT=IWORK(11) whih orresponds to the number of step rejetions aused by

onvergene test failures.

Itshould benoted thatthe hoie of dierent Newtoniterationmatrieswithinthe

Newton iterationis not availableinthe ode RADAU5. Furthermore,the ode GEOMS

oersthe possibilityof a ertainnumberof updates of the Newtoniterationmatrix

duringthe Newtoniteration inside of one integration step, see the following.

The onvergene rate of the simplied Newton method is investigated in detail in

[6℄, see also [17, 27℄. One important question in the use of an iterative method for

solvingnonlinear systemsinside anintegrationproess iswhentostop theiteration

suh that the obtained auray of the omputed solution of the nonlinear system

is within the presribed tolerane without performing too many Newton iteration

steps.

Theonvergene estimation andthe stoppingriterionimplementedinGEOMS isde-

sribed in [18℄ and adopted from the ode RADAU5 [17, 18℄. The estimation of the

onvergene isbased onthe weighted root square norm

|| · || ^sc

^{whihis dened for a}

ζ ∈ R ⁿ

by

|| ζ || ^sc = v u u t 1 n

n

X

i=1

ζ i

sc i

2

(19)

with

sc i = ATOL (i) + max( | x ki | , | x k+1i | ) RTOL (i)

^{, see [18℄. This norm allows to pre-}

sribe thatsome solutionomponentshave tobemore preiselyapproximated than

and relative tolerane, respetively. For more details on the error estimation and

the stoppingriterionof the Newton iterationwe referto [18, 38℄.

In the ase of a very slow onvergene or, in partiular, in the ase of divergene,

the number of Newton iteration steps has to be restrited by a maximal number

k max = NIT = IOPT ( 3 )

^{. Thus, the Newton iteration will stop} unsuessfully if a)

the stopping riterion is not satised within the maximal number

k max

^{of allowed}

Newton iterationsteps, orif b) the iterationdiverges.

In ase a) the user has to deide whether the whole integration step has to be re-

jeted beause of onvergene failures and toberepeated with aredued step size,

or if the Newton iteration should be ontinued with an updated Newton iteration

matrix. InGEOMS this deision ismade by dening the maximal number of updates

in the option IOPT(10)=NWTUPD. However, several numerial results suggest that

the number of allowed updates should not exeed 1.

It shouldbe notedthat the possibility of anupdate of the Newtoniterationmatrix

withinthe Newton iterationisnot available inthe ode RADAU5.

Duringthe Newtoniterationalinearsystem has tobesolved ineahstep. Thishas

to be done in an eient but stable way. The ode GEOMS oers the possibility to

deompose the dierential part and the algebrai part via dierent deomposition

methods. The user has to speify inthe option IOPT(11)=DECOMPC if the algebrai

part, i.e., the Jaobian of the onstraints, should be deomposed by use of the LU

deompositionwithfullpivoting(IOPT(11)=1),byaQRdeompositionwithpivot-

ing(IOPT(11)=2), orbyaSVdeomposition(IOPT(11)=3). Heuristiallyseen,the

LUdeompositionwith(partial)pivotingisagoodompromiseonerningeieny

and stability. Therefore, it isthe default inGEOMS, although,the SV deomposition

oersexellentstability properties but is more expensive.

Furthermore, with the option IOPT(12)=DECOMPD, the user an speify how to de-

ompose the dierentialpart. By setting IOPT(12)=0 the LU deomposition with

partialpivoting isused and by setting IOPT(12)=1 the QRdeomposition is used.

Remark 3.3 a)Theseparate deompositionimplemented inGEOMS has theadvan-

tage that the deompositionof the algebraipart an bedone independentlyof the

stepsize

h

^{. Onlythe}deompositionofthedierentialparthastobedoneseparately dependingon

h

^{. Inpartiular,iftheNewtoniterationhasonvergene problemsand}

thealgorithminterruptstheNewtonfortoreduethestepsize,thentheinformation

withrespettothe algebraipartmaybereyledwhihsavesomputationalwork.

b)Forthe linear algebraomputations like QRdeompositions and SVdeomposi-

tions we use BLAS 1

(Basi Linear Algebra Subprograms) [25℄ and LAPACK 2

(Linear

Algebra PACKage) [2℄ subroutines.

⊳

Forstrangeness-freedierential-algebraisystems insemi-impliitformlikethepro-

jeted-strangeness-freeform(12) thesaling ofthe algebraionstraintswith

1/h

^is

1

BLAS - http://www.netlib.org/blas/

2

LAPACK - http://www.netlib.org/lapak/

reommended in [32℄,where

h

^{is the urrent step size. Sine the numerial integra-}

tion of the equationsof motion in GEOMS isbased on the projeted-strangeness-free

formulationof the equationsof motion, the onstraintsare saledby

1/h

^.

The step size ontrol of the integration proess is a very sensitive topi in the im-

plementation of numerial algorithms for the integration of ODEs as well as for

DAEs. An overview overseveral step size ontrolstrategies isgiven in[37℄, see also

[4, 5, 11, 18℄. The ode GEOMS works with two dierent step size ontrol strategies

asusedinthe odeRADAU5, butadapted tothestruture ofthe equationsofmotion

(1). The basis for a step size ontrol mehanism is a loal error estimation. For

more details we refer to [18℄. The error estimation is implemented in the subrou-

tine GEERREST. For the hoie of a new step size for the next integration step or a

repeated integration step two possibilities are implemented inGEOMS whih have to

beseleted by use ofthe option IOPT(8)=PRED.WithIOPT(8)=2 the lassial step

sizeontroller developed in[11℄ isusedandwith IOPT(8)=1thepreditivestep size

ontroller,developed by Gustafssonin[16℄, isused. Thepreditivestep sizeontrol

is not possible in the rst step, so, the lassial step size ontroller will be used

instead. The preditive step size ontroller needs slightly more work and storage

than the lassial step size ontroller but is more exible in adaptating the step

size. Byuse of the preditive step size ontroller a faster redution of the step size

without step rejetions is possible than by use of the lassial step size ontroller.

This leads to a possible redution of the overall amount of omputation by use of

thepreditivestepsize ontroller. Experiments suggestthatthepreditivestepsize

ontroller seems to produe safer results for simple problems. On the other hand,

thehoieofthe lassialontrolleroftenproduesslightlyfaster runs,see also[18℄.

The preditive step size ontroller willbe used inGEOMS by default.

Sine the ode GEOMS is based on the ombination of disretization and regulariza-

tion tothe projeted-strangeness-freeformulation ofthe equationsof motionwhih

is inuened by the hoie of the seletors

S p

^and

S v

^{, see Theorem 2.4, an eient}

omputationofthese seletors isalsoimportantand willbedisussed inthe follow-

ing.

In general, it is not neessary to reompute seletors in every integration step, see

Remark 2.7.

If the LU deomposition is used for the dierential part then it is possible to de-

ide whether the determination of the seletors is done in eah integration step

(IOPT(13)=SELCOMP=1)ortheseletorsarekeptonstantforthoseintegrationsteps

wherethe pivoting inthe algebraipart doesnot hange (IOPT(13)=0). The latter

ase is the default.

The ode GEOMS oers the possibility to integrate the equations of motion of form

(1)withpossiblyredundantonstraints. As disussed inthe literature[26℄, see also

Remark 2.3, the solution may not be unique in this ase, but under ertain ondi-

tions the nonuniqueness is only restrited tothe Lagrangemultipliers

λ

^,

µ

^{, and}

w

^.

Formore details see [38℄.

is the detetion of the degree of redundany, i.e., the determination of the rank of

the Jaobianassoiated withthe onstraints. The reliablenumerialdetermination

of the rank of a matrix isa deliate taskand the SV deomposition isa ommonly

used tool for doing this. Therefore, the numerial integration of equations of mo-

tion with redundant onstraints is only allowed via the SV deomposition for the

onstraints,i.e., IOPT(11)=DECOMPC=3.

The rank of the onstraints will be determined in every integration step. If it is

deteted in the rst step that the onstraints are redundant, a reliable numerial

integration requires the use of the SVdeomposition atleast for the deomposition

of the onstraints. Furthermore, if a possibly hange of the rank from one step to

anotheris deteted, then the integration possibly has reahed a singularpoint and

willbe stopped with anerror message.

If the equations of motion have solution invariants (3), then it is often desirable

to preserve these solution invariants expliitly. GEOMS is able to preserve solution

invariantsif they are provided by the user as equations(3)inthe RDA(18). See the

Example 4.1.

Theusermayrestritthemaximalnumberofallowedintegrationstepsbysettingthe

option IOPT(6)=NMAX. The default value of NMAX is 100000. Furthermore, the user

may fore the ode to exploit some speial strutures of the problem. If the prob-

lemisautonomoustheamountofomputationalworkfor thenumerialintegration

may be redued. By settingIOPT(14)=AUTONOM=1 the user tellsthe ode that the

problemisautonomous and theode GEOMSexploits this inthe integrationproess.

The default is IOPT(14)=0, i.e., the problem is not autonomous. In partiular, if

the massmatrix isonstant and/ordiagonalalarge amountof omputationalwork

an be saved. Therefore, the user an speify by use of IOPT(15)=MASSTRKT if the

mass matrix is diagonaland onstant IOPT(15)=4, full and onstant IOPT(15)=3,

diagonal and time and/or state dependent IOPT(15)=2, or full and time and/or

state dependent IOPT(15)=1 (default).

4 Numerial experiments

In the following we will demonstrate the appliability and the performane of the

new solver GEOMS. The integration with GEOMS willbe performedfor three dierent

formulations of the regularized equations of motion. First, the numerial results

obtainedwithGEOMS usingtheprojeted-strangeness-freeform(12)ofthe equations

of motion will be abbreviated by GEOMS(psfEoM). Seond, the numerial results

obtained with GEOMS without providing the onstraints on aeleration level, see

optionIOPT(5) onPage 11,will beabbreviated by GEOMS(pEoM1). Furthermore,

ifthesolutionoftheonsideredexamplesatisessomesolutioninvariantswewilluse

the projeted-strangeness-free form ofthe equationsof motion with expliitforing

ofthesolutioninvariantsinadditiontothe twoformulationsabove,see Page4. The

The numerialintegrations are done onanAMD AthlonXP 1800+, 1533 MHz.

Letusnotethatwewillabstainfromtheuseofphysialunitslikemetersorseonds.

Example 4.1 The mathematial pendulum: In Example 2.8 we introdued

theequationsofmotionofthemathematialpendulumandwedidregularizethemto

the projeted-strangeness-freeform(14) whihisused forthe numerialintegration

via GEOMS.

Forthenumerialsimulationsof themovementweused themass

m = 1

^{,the length}

L = 1

^{, and the} gravitational aeleration

g = 13.75

^{. Let us note that we did}

modify the gravitational aeleration to approximately

g = 13.75

^{suh that the}

exat solution has a periodof

2

^{whih allows the omparison of the auray every}

period.

0 100 200 300 400 500 600 700 800 900 1000

−20

−15

−10

−5 0

5 x 10 ⁻⁴ Total energy

simulation time

total energy

prescribed tolerance 1e−07

PSfrag replaements

GEOMS(psfEoM) GEOMS(psfEoM1) GEOMS(psfEoM+I)

RADAU5(EoM) RADAU5(EoM2) RADAU5(EoM1)

RADAU5(GGL)

ODASSL(oEoM)

DASSL(EoM1) DASSL(EoM1)

MEXAX HEDOP5

Figure1: MathematialPendulum: Conservation of thetotal energy by thenumer-

ialsolutionsfor presribed RTOL=ATOL=

10 ⁻⁷

^{onthe time domain}

I = [0, 1000]

Themathematialpendulummodeledasin13representsamehanialsystemwhih

onserves the total energy. This total energy isgiven by

E(p, v) = 1

2 m(v ² ₁ + v ₂ ² ) + mgp 2

⁽²⁰⁾

and is onserved suh that

0 = E(p(t), v(t)) − E 0 = e(p, v)

^with

E 0 = E(p 0 , v 0 )

⁽²¹⁾

for

t ∈ I

and every solution of the equationsof motion (13).

Letusonsiderthe holonomionstraints(13) andtheirderivatives, whihrestrit

tion of the total energy (21). Wehave

0 = p ² ₁ + p ² ₂ − L ² ,

^(22a)

0 = 2p 1 v 1 + 2p 2 v 2 ,

^(22b)

0 = 2v ₁ ² + 2v ² ₂ − 2p ₂ g − 4

m (p ² ₁ + p ² ₂ )λ ₁ ,

⁽²²⁾

0 = 1

2 m(v ² ₁ + v ₂ ² ) + mgp ₂ − E ₀ .

^(22d)

The onstraints (22) are nonredundant for all

p

^,

v

^{, and}

λ

^{satisfying (22). In par-}

tiular, in addition to the holonomi onstraints and their derivatives the energy

onservation restrits the solution as well. The dimension of the solution manifold

withthe energy onservationis thereforesmallerthanwithoutthe energy onserva-

tion.

Foromparison,inFigure1thetotalenergyinthenumerialsolutionisdepited. In

additiontoGEOMS,thenumerialsolutionisomputedwithRADAU5[17,18℄fordier-

entformulations,i.e., (EoM)theequationsof motion(1)ofd-index3,(EoM2)thed-

index2formulation(usingthe onstraintsonveloitylevelinsteadofthe holonomi

onstraints),(EoM1)thed-index1formulation(usingtheonstraintsonaeleration

level instead of the holonomi onstraints), and (GGL) the Gear-Gupta-Leimkuhler

formulation, see [13℄. Furthermore, the solution is omputed with ODASSL [9, 10℄,

DASSL [4, 31℄, MEXAX [28℄, and HEDOP5 [3℄. Expeting GEOMS(psfEoM+I) the nu-

merially omputed total energy is far from being onstant. This an be expeted

beausethe energyonservationisontainedasanequationinthe usedformulation

and is therefore expliitely fored during the numerialintegration. However, even

the other numerial results obtained with GEOMS satisfy the onservation of total

energy very aurately. The preserving of the total energy yields a stabilization of

the solution.

In the Figures 2 and 3 the eieny is depited, i.e., the relation between the

obtained auray and the onsumed omputation time of the dierent used for-

mulations. Obviously, the integration with use of GEOMS based on the projeted-

strangeness-free formulation(14) plus solutioninvariantsGEOMS(psfEoM+I) oers

thebestperformaneforthisexample. NotethattheapproximationoftheLagrange

multipliersby GEOMS(psfEOM+I)is muh better than of the other results.

Averyimportantfatforthenumerialintegrationandthestabilityofthenumerial

algorithmsregarding the integration of DAEs is the satisfation of the onstraints,

inluding the hidden onstraints. In Figure 4 the residual of the onstraints of

positionlevel,ofveloitylevel,andofaelerationleveldependingonthesimulation

time isdepited. As one an see, GEOMS satises allonstraintswell.

Above we disussed the strategy for the determination of appropriate seletors,

onerningIOPT(13), see Page 18. Furthermore,the projeted-strangeness-free for-

mulation(14)ofthependulumhasbeen developedinExample2.8andthehoieof

theseletors

S p

^and

S v

^{has beenonsidered. Wehavestatedabovethat,inpriniple}

the seletors may be kept onstant as long as the deviation of the pendulum does