Lie group integrators

0.25 MB·

by  Brynjulf Owren

Lie group integrators
BrynjulfOwren
6
1
0
2

n
a
J

8
1
Abstract In this survey we discuss a wide variety of aspects related to Lie group

integrators.Thesenumericalintegrationschemesfordifferentialequationsonman-
]
A ifolds have been studied in a general and systematic manner since the 1990s and
N theactivityhassincethenbranchedoutinseveraldifferentsubareas,focussingboth
. ontheoreticalandpracticalissues.Fromtwoalternativesetups,usingeitherframes
h
or Lie groupactions on a manifold,we here introducethe most importantclasses
t
a ofschemesusedtointegratenonlinearordinarydifferentialequationsonLiegroups
m
andmanifolds.Wedescribeanumberofdifferentapplicationswherethereisanatu-
[ ralactionbyaLiegrouponamanifoldsuchthatourintegratorscanbeimplemented.

Anissuewhichisnotwellunderstoodistheroleofisotropyandhowitaffectsthe
1
v behaviourofthenumericalmethods.TheordertheoryofnumericalLiegroupinte-
4 gratorshasbecomeanadvancedsubtopicinitsownright,andherewegiveabrief
6 introduction on a somewhat elementary level. Finally, we shall discuss Lie group
6
integrators having the property that they preserve a symplectic structure or a first
4
integral.
0
.
1
0
6 1 Introduction
1
:
v
[19].LeonhardEulerisundoubtedlyoneofthemostaccomplishedmathematicians
i
X ofalltimes,andthemodernthemeNumericalmethodsforLiegroupscanbetraced
r backtoEulerinmorethanonesense.Indeed,thesimplestandpossiblymostlyused
a
numerical approximation method for ordinary differential equations was first de-
scribed in Euler’swork InstitutionumCalculi Integralis(1768),VolumenPrimum,
Ch VII) and bears the name Euler’s method.And undoubtedly,the most used test
BrynjulfOwren
NorwegianUniversityofScienceandTechnology,AlfredGetzvei1,N-7491Trondheim,Norway,
e-mail:[email protected]
1
2 BrynjulfOwren
caseforLiegroupintegratorsistheEuler’sfreerigidbodysystem,whichwasde-
rivedinhisamazingtreatiseonMechanicsin1736.
Intheliteratureonstructuralmechanics,Liegroupintegratorshavebeenaround
fora longtime, butgeneral,systematic studiesofnumericalintegratorsfor differ-
entialequationsonLiegroupsandhomogeneousmanifoldsbeganasrecentasthe
1990s.SomenotableearlycontributionswerethoseofCrouchandGrossman[24]
andLewisandSimo[40].A seriesofpapersbyMunthe-Kaas[49, 50,51] caused
anincreasedactivityfromthelateninetieswhenalargenumberofpapersappeared
overashortperiodoftime.Manyoftheseearlyresultsweresummarisedinasurvey
paperbyIserlesetal.[34].TheworkonLiegroupintegratorshasbeeninspiredby
many subfields of mathematics. Notably, the study of order conditions and back-
ward error analysis uses results from algebraic combinatorics,Hopf algebras, and
hasmorerecentlybeenconnectedtopostLiealgebrasbyMunthe-Kaasandcoau-
thors,seee.g.[52].Inconnectionwiththesearchforinexpensivecoordinaterepre-
sentations of Lie groupsas well as their tangentmaps, the classical theory of Lie
algebrashas been put to use in many differentways. The theory of free Lie alge-
bras [63] has been used to find optimal truncationsof commutatorexpansionsfor
generalLiealgebras,seee.g.[53].Forcoordinatemapstakingadvantageproperties
ofaparticularLiealgebra,toolssuchasrootspacedecomposition[61]andgener-
alisedpolardecompositions[37]havebeenapplied.Alsothereisofcourseastrong
connectionbetween numericalmethodsfor Lie groupsand the area of Geometric
Mechanics.Thisconnectionisoftenusedinthesetuporformulationofdifferential
equationsinLiegroupsorhomogeneousmanifoldswhereitprovidesanaturalway
of choosing a group action, and in order to construct Lie group integratorswhich
aresymplecticorconserveaparticularfirstintegral.
Inthispaperwe shalldiscussseveralaspectsofLiegroupintegrators,we shall
however not attempt to be complete. Important subjects related to Lie group in-
tegrators not covered here include the case of linear differential equations in Lie
groupsandthemethodsbasedonMagnusexpansions,Ferexpansions,andZassen-
haussplittingschemes. These aremethodsthatcouldfit wellinto a surveyonLie
groupintegrators,butforinformationonthesetopicswereferthereadertoexcel-
lentexpositionssuchas[4,33].Anothertopicweleaveouthereisthatofstochastic
Liegroupintegrators,seee.g.[42].Weshallalsofocusonthemethodsandthethe-
ory behind them rather than particularapplications, of which there are many.The
interestedreadermaycheckoutthereferences[8,9,18,35,65,66].
In the next section we shall define a compact setup of notation for differential
equationsondifferentialmanifoldswith a Lie groupaction.Thenin Section3 we
discusssomeofthemostimportantclassesofLiegroupintegratorsandgiveafew
examplesofmethods.Section4brieflytreatsaselectionofgroupactionswhichare
interestinginapplications.TheninSection5weshalladdresstheissueofisotropy
inLiegroupintegrators,inparticularhowthefreedomofferedbytheisotropygroup
caneither beusedto reducethe computationalcostofthe integratororbe usedto
improvethe qualityofthe solution.We take ourownlookat ordertheoryandex-
pansionsintermsofageneralisedformofB-seriesinSection6.Finally,insections
Liegroupintegrators 3
7 and 8 we consider Lie groupintegratorswhich preserve a symplectic form or a
firstintegral.
2 The setup
Let M be a differentiable manifold of dimension d <¥ and let the set of smooth
vector fields on M be denoted X(M). Nearly everything we do in this paper is
concerned with the approximation of the h-flow of a vector field F ∈X(M) for
somesmallparameterhusuallycalledthestepsize.Inotherwords,weapproximate
thesolutiontothedifferentialequation
d
y˙= y= F| , F ∈X(M). (1)
dt y
Crouch and Grossman [24] used a set of smooth frame vector fields E1,...,En ,
n ≥donM,assuming
span(E1|x,E2|x,...,En |x)=TxM, foreachx∈M.
Itcanbeassumedthattheframevectorfieldsarelinearlyindependentasderivations
oftheringF(M)ofsmoothfunctionsonM,andwedenotetheirR-spanbyV.Any
smoothvectorfieldF ∈X(M)canberepresentedbyn functions f ∈F(M)
i
n
(cid:229)
F| = f(x)E| (2)
x i i x
i=1
wherethe f arenotnecessarilyunique.Wethenhaveanaturalaffineconnection
i
(cid:209) G=(cid:229) F(g)E
F i i
i
which is flat with constant torsion t ((cid:229) f E ,(cid:229) gE)=(cid:229) f g[E,E ]. For later,
j j j i i i i,j j i i j
we shall need the notion of a frozen vector field relative to the frame. The freeze
operatorFr:M×X(M)→V isdefinedas
(cid:229)
Fr(x,F):=F = f(x)E
x i i
i
We notethatthetorsioncanbe definedbyfreezingthevectorfieldsandthentake
theLie-Jacobibracket,i.e.
t (F,G)| =[Fr(x,F),Fr(x,G)]
x
Another setup is obtained by using a Lie group G acting transitively from the
leftonM [51]. TheLie algebraof G is denotedg. AnyvectorfieldF cannowbe
representedviaamap f :M→gandtheinfinitesimalactionr :g→X(M)
4 BrynjulfOwren
d
F| =r ◦f(x)| , r :g→X(M), r (x )| = exp(tx )·x (3)
x x x
dt
(cid:12)t=0
(cid:12)
Wenotethatthemap f isnotnecessarilyunique. (cid:12)
(cid:12)
3 Types ofschemes
There is now a large variety of numericalintegrationschemes available, typically
formulated with either of the setups of the previous section. In what follows, we
assumethatafinitedimensionalLiegroupGactstransitivelyonamanifoldMand
theLiealgebraofGisdenotedg.
3.1 Schemes ofMunthe-Kaas type
Using the secondsetup, a powerfulway of derivingnumericalintegratorsdevised
byMunthe-Kaas[51]isto
1. InaneighborhoodU⊂gof0introducealocaldiffeomorphismy :U→G,such
that
y (0)=1∈G, T y =Id
0 g
2. Observe that the map l (v)=y (v)·y is surjective on a neighborhoodof the
y0 0
initialvaluey ∈M.
0
3. ComputethepullbackofthevectorfieldF =r ◦f alongl
y0
4. ApplyonestepofastandardnumericalintegratortotheresultingproblemonU
5. MaptheobtainedapproximationbacktoMbyl
y0
Eventhoughtheideaisverysimple,thereareseveraldifficultiesthatneedtobe
resolved in order to obtain fast and accurate integrationschemes from this proce-
dure.
Observethatthederivativeofy canbetrivialisedbyrightmultiplicationofthe
Liegroup,suchthat
Tuy =TRy (u)◦dy u, dy u:g→g.
WiththisinmindwesetouttocharacterisethevectorfieldonU⊂g,thisisasimple
generalisationofaresultin[51]
Lemma1.LetMbeasmoothmanifoldwithaleftLiegroupactionL :G×M→M
andletgbetheLiealgebraofG.Lety :g→Gbeasmoothmap,y (0)=1.Fixa
pointm∈M,andsetL =L (·,m)sothat
m
r (x )| =T L (x )
m 1 m
Liegroupintegrators 5
SupposeF ∈X(M)isoftheform
F| =r (x (m))| , forsomex :M→g.
m m
Definel (u)=L (y (u),m).ThenthereisanopensetU ⊆gcontaining0suchthat
m
thevectorfieldh ∈X(U)definedas
h | =dy −1(x ◦l (u))
u u m
isl -relatedtoF.
m
The original proof in [51] where y =exp was adapted to general coordinate
maps in [61]. One step of a Lie group integrator is obtained just by applying a
classical integrator, such as a Runge–Kutta method, to the corresponding locally
definedvectorfieldh ong.ARunge–Kuttamethodwithcoefficients(A,b)applied
totheproblemy˙=h (y)inalinearspaceisofthefollowingform
s
(cid:229)
Y =y +h a k , r=1,...,s,
r 0 rj j
j=1
k = h | , r=1,...,s,
r Yr
s
(cid:229)
y =y +h b k
1 0 r r
r=1
Here y is the initial value, h is the step size, and y ≈y(t +h) is the approxi-
0 1 0
matesolutionattimet +h.Theparameters≥1iscalledthenumberofstagesof
0
the method.If the matrix A=(a ) is strictly lowertriangular,then the methodis
rj
calledexplicit.Theapplicationofsuchamethodtothetransformedvectorfieldof
Lemma1cabewrittenoutasfollows
s
u =h(cid:229) a k˜ , k =x ◦l (u ), k˜ =dy −1(k ), r=1,...,s, (4)
r ij j r y0 r r ur r
j=1
s
v=h(cid:229) b k˜ y =l (v). (5)
r r 1 y0
r=1
AmajoradvantageofthisapproachcomparedtoothertypesofLiegroupschemes
isitsinterpretationasasmoothchangeofvariableswhichcausestheconvergence
ordertobe(atleast)preservedfromtheunderlyingclassicalintegrator.Asweshall
see later, one generally needs to take into account additional order conditions to
accountforthefactthatthephasespaceisnotalinearspace.
3.1.1 Choosingtheexponentialmapascoordinates,y =exp
ThefirstpapersbyMunthe-KaasonLiegroupintegrators[49,50,51]allusedy =
expascoordinateson the Liegroup.In thiscase, thereare severaldifficultiesthat
6 BrynjulfOwren
need to be addressed in order to obtain efficient implementationsof the methods.
Oneisthe computationofthe exponentialmapitself. Formatrixgroups,thereare
a largenumberofalgorithmsthatcanbeapplied,see e.g.[47, 48].In[15, 16]the
authorsdevelopedapproximationstotheexponentialwhichexactlymapmatrixLie
subalgebrastotheircorrespondingLiesubgroups.Anotherissuetobedealtwithis
thedifferentialoftheexponentialmap
dexp :=T R ◦T exp
u exp(−u) exp(−u) u
Lemma2.(Tangentmapofexp:g→G).Letu∈g,v∈T g≃g.Then
u
d
T exp(v)= exp(u+sv)=TR ◦dexp v=dexp (v)·exp(u)
u ds exp(u) u u
(cid:12)s=0
(cid:12)
where (cid:12)
(cid:12)
1 ez−1
dexp (v)= exp(rad )(v)dr= (v)
u u z
Z0 (cid:12)z=adu
(cid:12)
(cid:12)
(cid:12)
Proof. Lety (t)=exp(t(u+sv))suchthat
s
d
T exp(v)= y (1)
u s
ds
(cid:12)s=0
(cid:12)
Butfornowwedifferentiatewithrespecttot(cid:12)toobtain
(cid:12)
d
y˙ := y (t)=(u+sv)y (t) (6)
s s s
dt
Wealsonotethaty (t)=exp(tu)+O(s)ass→0.From(6)wethenget
s
y˙−uy =svetu+O(s2)
s
andtheintegratingfactore−tu yields
d
e−tuy =se−tuvetu+O(s2)
s
dt
(cid:0) (cid:1)
Integrating,andusingthaty (0)=Id,weget
s
t
y (t)=etu+s eruve−ruetudr+O(s2)
s
0
Z
andso
d 1 1
y (1)= eruve−rudreu= eradu(v)dreu
s
ds(cid:12)s=0 Z0 Z0
(cid:12)
wherewehaveuse(cid:12)dthewell-knownidentityAd =exp(ad )inthelastequality.
(cid:12) exp(u) u
Formally,wecanwrite
Liegroupintegrators 7
1 1 ez−1
eradu(v)dr= erz| (v)dr= (v)
z=adu z
Z0 Z0 (cid:12)z=adu
(cid:12)
(cid:12)
Itisoftenusefultoconsiderthedexp-mapasaninfinites(cid:12)eriesofnestedcommu-
tators
1 1 1 1
dexp (v)=(I+ ad + ad2+···)(v)=v+ [u,v]+ [u,[u,v]]+···
u 2! u 3! u 2 6
In(4)itistheinverseofdexpwhichisneeded.Notethatthefunction
ez−1
f (z)=
1
z
isentire,thismeansthatitsreciprocal
z
ez−1
is analytic where f (z) 6=0. In particular this means that 1 has a converging
1 f1(z)
Taylor series about z=0 in the open disk |z|<2p . This series expansion can be
showntobe
¥
z =1− z +(cid:229) B2k z2k
ez−1 2 (2k)!
k=1
whereB aretheBernoullinumbers,thefirstfewofthemare:B = 1,B =− 1 ,
2k 2 6 4 30
B = 1 ,B =− 1 ,B = 5 .Themap
6 42 8 30 10 66
v=dexp−1(w) (wheneverw=dexp (v))
u u
isgivenpreciselyas
z 1 B
dexp−1(w)= (w)=w− [u,w]+ 2[u,[u,w]]+··· (7)
u ez−1 2 2!
(cid:12)z=adu
(cid:12)
(cid:12)
Weobservefrom(4)thatonen(cid:12)eedstocomputedexp−1(k )andthateachu =O(h).
ur r r
Thismeansthatonemayapproximatetheseriesin(7)byafinitesum,
dexpinv(u,w,m)=w−1[u,w]+(cid:229)m B2k ad2k(w).
2 (2k)! u
k=1
Onehasdexpinv(u,w,m)∈gforeverym≥0andfurthermore
dexp−1(k )−dexpinv(u ,k ,m)=O(h2m+1)
ur r r r
Aslongastheclassicalintegratorhasorder p≤2m+1,theresultingMunthe-Kaas
scheme will also have order p. Thereexists howevera clever way to substantially
reducethenumberofcommutatorsthatneedtobecomputedineachstep.Munthe–
8 BrynjulfOwren
KaasandOwren[53]realisedthatonecouldformlinearcombinationsofthestage
derivativesk˜ in(4)suchthat
r
r
Q = (cid:229) s k˜ =O(hqr)
r r,j j
j=1
forq aslargeaspossibleforeachr.ThenthesenewquantitiesQ wereeachgiven
r r
thegradeq andoneconsideredthegradedfreeLiealgebrabasedonthisset.The
r
resultwasasignificantreductioninthenumberofcommutatorsneeded.AlsoCasas
andOwren[13]providedawaytoorganisethecommutatorcalculationstoreduce
evenfurtherthecomputationalcost.HereisaRunge–KuttaMunthe–Kaasmethod
oforderfourwithfourstagesandminimalsetofcommutators
k =hf(y ),
1 0
k =hf(exp(1k )·y ),
2 2 1 0
k =hf(exp(1k −1[k ,k ])·y ),
3 2 2 8 1 2 0
k =hf(exp(k )·y ),
4 3 0
y =exp(1(k +2k +2k +k −1[k ,k ]))·y .
1 6 1 2 3 4 2 1 4 0
Forlaterreference,wealsogivetheLie-Eulermethod,afirstorderLiegroupinte-
gratorgeneralisingtheclassicalEulermethod
y =exp(hf(y ))·y (8)
1 0 0
3.1.2 Canonicalcoordinatesofthesecondkind
Theexponentialmapisgenerallyexpensivetocomputeexactly.FormatrixLieal-
gebrasg⊆gl(n,F)whereFiseitherRorC,standardsoftwareforcomputingexp
numericallyhasacomputationalcostofn3totheleadingorder,andtheconstantin
frontofn3 maybeaslargeas20−30.Another,yetcompletelygeneralalternative
totheexponentialfunctionisconstructedasfollows:Fixabasisforg,saye ,...,e
1 d
andconsiderthemap
y :v e +···+v e 7→exp(v e )·exp(v e )···exp(v e ) (9)
1 1 d d 1 1 2 2 d d
Although it might seem unnaturalto replace one exponentialby many, one needs
to keep in mind that if the basis can be chosen such that its exponential can be
computedexplicitly,itmaystillbeanefficientmethod.Forinstance,inthegeneral
linearmatrixLiealgebragl(mF)onemayusethebasistobee =eeT wheree is
ij i j i
theithcanonicalunitvectorinRn.Then
exp(a e )=1+a e , i6= j, exp(a e )=1+(ea −1)e
ij ij ii ii
Liegroupintegrators 9
So computing(9) takes approximatelynd operationswhich is much cheaper than
computingtheexponentialofageneralmatrix.
Thedifficultylies howeverin computingthe map dy −1 in an efficientmanner.
u
Amethodforthiswasdevelopedin[61].Themethodologyisslightlydifferentfor
solvableandsemisimpleLiealgebras.Wehereoutlinethemainidea,fordetailswe
refertotheoriginalpaper[61].Differentiate(9)toobtain
d
dy u(v)=v1e1+(cid:229) vi Adeu1e1◦···◦Adeui−1ei−1(ei)
i=2
Themainideaistofindanequivalentexpressionwhichisacompositionofcheaply
invertibleoperators.Forthis,weintroduceaprojectorontothespanofthelastd−k
basisvectorsasfollows
d d
(cid:229) (cid:229)
P : ve 7→ ve
k i i i i
i=1 i=k+1
where we let P and P equal the identity operator and zero operator on g re-
0 d
spectively. We may now define a modified version of the Ad-operator, for any
u=(cid:229) ue ∈g,let
i i
Adeukek =(Id−Pk)+AdeukekPk
Thisisalinearoperatorwhichactsastheidentityoperatoronbasisvectorse, i≤k,
i
c
andonbasisvectorsei, i≥kitcoincideswithAdeukek
Definition1. An ordered basis (e ,...,e ) is called and admissible ordered basis
1 d
(AOB)if,foreachu=(cid:229) u e ∈gandforeachi=1,...,d−1,wehave
j j
Adeu1e1◦···◦AdeuieiPi=Adeu1e1◦···◦AdeuieiPi (10)
Thisdefinitionisexactlywhatisneededtowcritedy uasaccompositionofoperators
Proposition1.Ifthebasis(e ,...,e )isanAOB,then
1 d
dy u=Adeu1e1◦···Adeuded
Anotherimportantsimplificationcacnbeobtainecdif anabeliansubalgebrahof di-
mensiond−d∗ canbeidentified.Inthiscasetheorderedbasiscanbechosensuch
thath=span(ed∗+1,...,ed). Then Adeuiei|h fori>d∗ is the identityoperatorand
thereforeAdeuiei istheidentityoperatoronallofg.Summarizing,wehavethefol-
lowingexpression
c dy u−1=Ad−eu1d∗ed∗ ◦···◦Ad−eu11e1
Choosing typically a basis consisting of nilpotent elements, the inversion of each
c c
Adeuiei canbedonecheaplybymakinguseoftheformula
c Adeuiei =1+(cid:229)K ukki!adek, adKei+1=0.
k=1
10 BrynjulfOwren
For choosing the basis one may, for semisimple Lie algebras, use a basis known
as the Chevalley basis. This arises from the root space decomposition of the Lie
algebra
g=h⊕ ga . (11)
a ∈F
a
HereF isthesetofroots,andga istheone-dimensionalsubspaceofgcorrespond-
ingtotheroota ∈h∗,seee.g.Humphreys[32].histhemaximaltoralsubalgebra
of g and it is abelian. In the previous notation, the number of roots is d∗ and the
dimension of h is d−d∗. The following result whose proof can be found in [61]
providesatoolfordeterminingwhetheranorderedChevalleybasisisanAOB.
Theorem1.Let{b ,...,b }, d =d−ℓ,bethesetofrootsF forasemisimpleLie
1 d∗ ∗
algebrag.SupposethataChevalleybasisisorderedas
(eb1,...,ebd∗,h1,...,hℓ)
whereebi ∈gbi,and(h1,...,hℓ)isabasisforh.SuchanorderedbasisisanAOBif
kb +b =b , m<i<s≤d ,k∈N ⇒ b +b 6∈F , m<n≤i−1. (12)
i s m ∗ m n
HereF =F ∪{0}.
Example1.Asanexample,weconsiderA =sl(ℓ+1,C),commonlyrealizedasthe
ℓ
setof(ℓ+1)×(ℓ+1)-matriceswithvanishingtrace.Themaximaltoralsubalgebra
isthenthesetofdiagonalmatricesinsl(ℓ+1,C).Thepositiverootsaredenoted
{b ,1≤i≤ j≤ℓ}.
i,j
Lettinge betheithcanonicalunitvectorinCℓ+1,therootspacecorrespondingto
i
b hasabasisvector
i,j
eieTj+1∈gbi,j, 1≤i≤ j≤ℓ
whereasthenegativerootsareassociatedtothebasisvectors
ej+1eiT ∈g−bi,j.
As a basis for h, one may choose the matriceseeT −e eT , 1≤i≤ℓ. The re-
i i i+1 i+1
mainingdifficultynowistochooseanorderingofthebasissothatanAOBresults.
Asindicatedearlier,thebasisforhmaybeorderedasthelastones,i.e.withindices
rangingfromd∗+1=ℓ2+ℓ+1tod=ℓ2+2ℓ.Withtheconventioneb ∈gb , b ∈F ,
h=span(e ,...,e ),let
h1 hℓ
B=(ebi1,j1,...,ebim,jm,e−bi1,j1,...,e−bim,jm,eh1,...,ehℓ),
wherei ≤i ≤···≤i andm=ℓ(ℓ+1)/2.OnecanthenprovebyusingTheorem1
1 2 m
thatBisanAOB.