A new software package for linear differential-alg

ANEWSOFTWAREPACKAGEFORLINEARDIFFERENTIAL{ALGEBRAICEQUATIONSPETERKUNKEL,VOLKERMEHRMANNy,WERNERRATHy,ANDJORGWEICKERTyAbstract.WedescribethenewsoftwarepackageGELDAforthenumericalsolutionoflineardierential{algebraicequationswithvariablecoecients.Theimplementationisbasedonthenewdiscretizationschemeintroducedin[20].Itcandealwithsystemsofarbitraryindexandwithsystemsthatdonothaveuniquesolutionsorinconsistenciesintheinitialvaluesortheinhomogeneity.Thepackageincludesacomputationofallthelocalinvariantsofthesystem,aregularizationprocedureandanindexreductionschemeanditcanbecombinedwitheverysolutionmethodforstandardindex1systems.Nonuniquenessandinconsistenciesaretreatedinaleastsquaresense.Wegiveabriefsurveyofthetheoreticalanalysisoflineardierential{algebraicequationswithvariablecoecientsanddiscussthealgorithmsusedinGELDA.Furthermore,weincludeaseriesofnumericalexamplesaswellascomparisonswithresultsfromothercodes,asfarasthisispossible.Keywords.Dierential-algebraicequations,canonicalforms,backwarddierenceformulas,Runge-Kuttaformulas,leastsquareregularization,singularpencils,strangenessindex.AMSsubjectclassications.65L051.Introduction.WediscussthenewsoftwarepackageGELDAforthenumeri-calsolutionoflineardierential{algebraicequations(DAE’s)withvariablecoecientsE(t)_x(t)=A(t)x(t)+f(t);t2[t;t];(1)whereE;A2C([t;t];Cn;n);f2C([t;t];Cn)togetherwithaninitialconditionx(t0)=x0;t02[t;t];x02Cn:(2)HereC‘([t;t];Cn)denotesthesetof‘-timescontinuouslydierentiablefunctionsfromtheinterval[t;t]tothen-dimensionalcomplexvectorspaceCn.(Inthispaperwediscusscomplexfunctions.InthecaseofrealproblemsalltheresultsarevalidwhenCn,Cn;narereplacedbyRn,Rn;n,respectively.)Thetheoreticalanalysisofsuchsystems,regardingexistence,uniquenessofso-lutionsaswellasconsistencyofinitialconditionshasbeendiscussedindetailin[20,21,22].Wewillsurveyonlytherelevantpartofthisworkheretomaketheprocedurethatcomputestheinvariantsofthesystemtransparent.Themostim-portantinvariantintheanalysisoflinearDAE’sisthesocalledstrangenessindex,whichgeneralizesthedierentiationindex([4,5,11])forsystemswithundeterminedcomponentswhichoccur,forexample,inthesolutionoflinearquadraticoptimalcontrolproblemsanddierential{algebraicRiccatiequations,seee.g.[18,19,27].ItisknownthatmostofthestandardintegrationmethodsforgeneralDAE’srequirethesystemtohavedierentiationindexnothigherthanone,whichcorrespondstoavanishingstrangenessindex,see[21].IfthisconditionisnotvalidoriftheDAEhasundeterminedcomponents,thestandardmethodsasimplementedincodeslikeDASSLofPetzold[28]orLIMEXofDeuhard/Hairer/Zugck[7]oftenfail.FachbereichMathematik,CarlvonOssietzkyUniversitat,Postfach2503,D{26111Oldenburg,Fed.Rep.GermanyyFakultatfurMathematik,TUChemnitz{Zwickau,D{09107Chemnitz,Fed.Rep.Germany.ThisworkhasbeenpartiallysupportedbytheDeutscheForschungsgemeinschaftunderthegrantno.Me790/5-1\Dierentiell{algebraischeGleichungen.1Theimplementationofthenewsoftwarepackageisbasedontheconstructionofthediscretizationschemeintroducedin[20],whichrstdeterminesallthelocalin-variantsandthentransformsthesystemintoastrangeness-freeDAEwiththesamesolutionset.Wewillgiveabriefsurveyofthisscheme.Theresultingstrangeness-freesystemmaystillhavenonuniquenessinthesolutionsetorinconsistentinitialvaluesorinhomogeneities.Butthisinformationisnowavailabletotheuserandsystemswithsuchpropertiescanbetreatedinaleastsquaressense.InthecasethattheDAEisfoundtobeuniquelysolvable,wecancomputeaconsistentinitialvalueandapplythewell-knownintegrationschemesforDAE’s.InourpackagewehaveimplementedBDFmethods[3]andRunge-Kutta-schemes[14,15].2.Abriefsurveyofthebasicresults.Inthissectionwegiveabriefsurveyoftheresultsin[20,21,22].Themainresultsinthesepapersarethebasisfortheconstructionofthenewsoftwarepackage.Webeginwithtwoequivalencerelationsforpairsofmatrixvaluedfunctionsandmatrixpairs,whichplayacentralroleinthetheoryof(1).Definition1.Twopairsofmatrixfunctions(Ei(t);Ai(t))withEi;Ai2C([t;t];Cn;n);i=1;2arecalledequivalentifthereexistP2C([t;t];Cn;n)andQ2C1([t;t];Cn;n)withP(t);Q(t)nonsingularforallt2[t;t]suchthat(E2(t);A2(t))=P(t)(E1(t);A1(t))Q(t)_Q(t)0Q(t):(3)Forthedevelopmentofnumericalmethodsonealsoneedsacounterpartofthisequivalencethatcanbeobtainedlocallyataxedpoint.Ataxedpointt2[t;t]onecanchooseQ(t)and_Q(t)independently(see[12,20]).Thisleadstothedenitionoflocalequivalence.Definition2.Twopairsofmatrices(Ei;Ai);Ei;Ai2Cn;n;i=1;2arecalledlocallyequivalentiftherearematricesP;Q;B2Cn;nwithP;Qnonsingularsuchthat(E2;A2)=P(E1;A1)QB0Q:(4)NotethatthisisnotatransformationthatcanbeusedintheDAE(1),sincewecannottransformx(t)and_x(t)independently,butitismoreappropriatethantheusualequivalenceformatrixpencils,see[10].Forthelocalequivalence(4)thefollowingnormalformisprovedin[20]:Theorem2.1.LetE;A2Cn;n.Considerthefollowingmatricesandthespacesspannedbytheircolumns:(a)TbasisofkernelE(b)ZbasisofcorangeE=kernelE(c)T0basisofcokernelE=rangeE(d)Vbasisofcorange(ZAT):(5)Thenthequantities(withtheconventionrank;=0)(a)r=rankE(rank)(b)a=rank(ZAT)(algebraicpart)(c)s=rank(VZAT0)(str