Role of the particle method DFPM for solving linear equations

S. Edvardsson, M. Neuman, P. Edström

∗

Department of natural s ien es, MidSweden University, SE-871 88, Härnösand, Sweden

M. Gulliksson and H. Olin

Department of natural s ien es, Mid Sweden University, SE-851 70, Sundsvall, Sweden

The present work is on erned with exploiting a se ond orderdynami al system approa hfor

solvingequations. Theparti ularfo ushereisforrealvaluedlinearequationsystemshavingawide

eigenvalue spe trum,often en ounteredinappli ations as dis retization be omes dense. The Dy-

nami alFun tionalParti leMethod(DFPM)isdevelopedasaninterdis iplinary methodbetween

physi sandnumeri almathemati s.Themethodisthis ontextoptimizedtotakeadvantageof riti-

allydampedos illatorsresultingingoodperforman e. Convergen eisrea hedwhenalleigenvalues

are eitherpositiveor negative. The ase with omplexeigenvalues isalsostudied. Theparti ular

stru tureofthe matrixturnsoutto beunimportantfor onvergen e. Furthermore,DFPMis not

limitedwithrespe ttoaspe tralradiusasis ommonforiterativemethods.Arstorderdynami al

systemisalsostudiedand ompared. Itsperforman eisnot ompetitivein omparisonwithDFPM.

Theperforman eofDFPMs alesequallytothewellknown onjugategradient method,butla ks

someof itslimitationsimposedbymatrixsymmetry. Several physi altest examplesare provided

and omparedwithvariousexistentnumeri alapproa hes.

I. INTRODUCTION

Parti le methodsin omputational physi s

Thenumeri altreatmentthroughtheappli ation ofparti lemethods hasbeenaverya tiveresear hareaduring

many years but these methods have still not been exploited to their full potential [1, 2℄. Parti le methods are

on eptuallyattra tivesin etheyaresprungfromtheveryfundamental oreofphysi s. Thebasi assumptionisthat

allmatter onsistsofintera tingparti lesthatobey ertainlawsofmotion(su hasNewton's se ondlaw). Therst

parti lemethodsweredevelopedsoonafterthebirthoffast omputerssome40yearsago. Afewexamplesofparti le

treatmentsofdynami alobje tsin elestialme hani saregivenin[3,4℄. Thelong-termevolutionofdynami alobje ts

intheSolarsystem(orasmallgalaxy) anreadilybe al ulated[5℄. Thedis reteelementmethod(DEM)isamongthe

earlymethodswhi happearedinthe1970s[6℄. TodayDEMisoftenappliedtotreatengineeringproblemsingranular

anddis ontinuousmaterials. Anotherrelativelyearlyexampleisthemole ulardynami ssimulationte hnique(MDS)

that was developed during the 1980s. Today, millions of mole ules oratoms are more or less non-problemati to

simulateduetothesimpli ityofparallelizingtheN-bodyalgorithms[7℄. Smoothed-parti lehydrodynami s(SPH)is

yet anotherparti le te hniquethat dealswith simulations ofuids. SPH wasoriginally developed by Lu y [8℄ and

GingoldandMonaghan[9℄todealwith3Dproblemsinastrophysi s. SPHhasalsobeenadoptedtosolidme hani s

to study impa t fra tures in solids. The method is in this ontext abbreviated SPAM (smooth parti le applied

me hani s)[10℄. TherearealsoMDS-likeparti lemethods(HybridLatti eParti leModeling),whereLennard-Jones

potentialsorsimilar intera tion potentialsare utilized forquasi-parti les[11℄. Inthelimit whenthe quasi-parti les

approa h the atomi s ale the parti le method be omes an ordinary MDS method. The Car-Parrinello method is

anextremelywellknownparti le method fordynami al omputationof ele troni quantum statesduring mole ular

dynami ssimulation[12℄. Alsoindensityfun tionaltheory,parti lemethodshavebeenappliedto omputeele troni

states, see agoodreview by Payneet al. [13℄. Parti le methods haveoften beenfound attra tivefor runningfast

omputer graphi s of deformable obje ts[14, 15℄. A quite omprehensivereview onvarious parti le approa hesis

providedbyLiandLiu[2℄.

Inare entpaperbyEdvardssonetal. [16℄,itwasarguedthattheparti le on eptismoregeneralthanitsstandard

usage in physi s or me hani s. Inspired by the onne tionbetweenparti les in physi s and theirusual translation

into dierential equations,the opposite possibility to translatemathemati al problemsinto aquite generalparti le

s hemewasre ognized. Thedynami alfun tionalparti lemethod(DFPM)wasdeveloped. Edvardssonetal. applied

DFPMto eigenvalueproblems(S hrödingerequation)and ompared its omputationalspeedtostandardnumeri al

∗

Ele troni address: sverker.edvardssonmiun .se

librariessu hasARPACKandLAPACK.Forlargematrixsizesitwasseenthattheparti lemethodwasmu hmore

e ientin thedeterminationofafeweigenvaluesandeigenve tors. The onjugategradientmethodwastestedlater

dealingwiththesameS hrödingerproblemandfoundtobesome50%lesse ientthanDFPM[17℄. DFPMwasalso

demonstratedtoworkwellforthenon-linearS hrödingerequation[16℄. Convergen epropertieswereaddressed,and

in analogywithmany-parti le systemsin me hani s,theexisten eof apotentialminimumissu ientto guarantee

onvergen e.However,itwasre ognizedempiri allythatthereexistmany aseswhere onvergen estillo ursdespite

thela kofapotential. Thisfa t,andthedetailsrelatedto omputational omplexity,wasleftoutforfurtherstudies.

These detailswill bestudied in thepresentwork forthespe ial,but important, aseoflinearequationsystems. As

we shallsee theDFPM algorithm is attra tivedue to several reasons. The mostinterestingpoints arerelated to a

surprisingly generalappli abilityand robustness, omputational omplexity, easinessof implementation, interesting

Hamiltoniandynami sandstabilityofsymple ti integration.

Relatedmathemati al methods

Although the se ond order dynami al systemapproa hprovided by thepresentwork to solve linearproblems is

novel, there are related ideas given previouslyin mathemati s. Themost ommon method is the formulation of a

rstorderdynami alsystem,i.e., thesetup:

du/dt = F(u), u(0) = u 0

^{. As}

du/dt → 0

^{theoriginal problem}

F(u) = 0

is solved. Thus the idea is to solve a damped time dependent problem to identify the stationary solution. First

orderdynami alsystemshavefrequentlybeenappliedtosolvevariouskindsofequations

F(u) = 0

^{,bothasageneral}

approa handintendedforspe i mathemati alproblems. Anexampleofthisisthesolutionofanellipti PDEsu h

astheheatequation,seeSin ove andMadsen[18℄. Sin ethestationarystateissought,theevolutionofthesystem

is onsidered totakepla e inarti ial time. The on eptofarti ial timeis furtherdis ussedand analyzedin [19℄.

Otherworksareforexamplethedampedharmoni os illatorin lassi alme hani s,thedampedwaveequation,[20℄

andtheheavyballwithfri tion[21℄. Theseproblem settingsarespe i mathemati alexamplesofphysi al systems

andnotdevelopedtosolveequationsingeneral. In[22,23℄,iterativepro essestosolve,e.g.,eigenvalueproblemsare

onsideredas(gradientdriven)dynami alsystems. So alled titioustimeisusedin[24℄whereaDiri hletboundary

valueproblem ofquasilinearellipti equation is solvedbyusing the on eptof a titioustimeintegration method.

The inverse problem of re overing a distributed parameter model is onsidered in [19℄ using the rst order ODE

attainedfrom thene essaryoptimality onditionsoftheinverseproblem. Anotherapproa h isthat of ontinuation,

see [25℄for anintrodu tion, where

F(u) = 0

^{is embedded in afamilyofproblems dependingonaparameter}

s

^{, i.e.,}

F(u; s) = 0

^{. Thesolutionto}

F(u) = 0

^{isfoundbysolvingasequen eofproblemsforvaluesof}

s

^{de reasingfrom}

1

^to

0

^{. Further,seeNo edalandWright[26℄ foradis ussioninthe ontextof}optimization.

Outlineofthe paper

Thepaperstartswith anintrodu tionto thegeneralideasbehindDFPMand itsobvious onne tionsto lassi al

me hani s. The versatility of DFPM is emphasized. The paper then ontinues to treat spe i ally the

Ax = b

problem. The onvergen eand onvergen erateis investigated. This analysisis expe tedto beusefulalsolaterfor

otherrelatedproblems(e.g. eigenvalueproblems). The onne tionswiththeparti les hemeinphysi sisemphasized

throughout the arti le. It turns out that physi al arguments and many-parti leviews provide importantshort uts

in thederivations andanalyzes. A omparisonbetweenDFPMand therelatedrst orderdynami al systemisalso

made. Thearti le is on luded with numeri al examplesand omparisons with existent methods in the numeri al

linearalgebraliterature. Inorder toenhan e readabilityof thepaperwehaveprovidedanAppendix where several

relateddetailsaregiven.

II. THEDYNAMICALFUNCTIONAL PARTICLE METHOD

In the following we makea brief summary for the reader about DFPM as a quite general approa h for solving

equations. Let

F

^{beafun tional and}

u = u(x), u : X → R ^k , k ∈ N

^{and onsider theabstra t equation}

F(u) = 0

⁽¹⁾

that ouldbe,e.g.,afun tional,dierential,integralorintegro-dierentialequation. Further,atime parameter

t

^is

introdu edthat belongs tosome(unbounded)interval

T

^{andadynami alsystemin}

u = u(x, t), u : X × T → R ^k

^is

formedas

F(u) = η¨u + µ ˙u.

⁽²⁾

wherethedotsarethestandardnotationfortimederivativesin me hani s. Thesymbols

η = η(x, t)

^and

µ = µ(x, t)

are themass anddamping parameters,respe tively. The main ideais to solvethe originalequation (1) by instead

solving(2) in su h away that

˙u, ¨ u → 0

^when

t → T, T < ∞

^{, i.e.,}

u(x, t)

^{approa h}

u(x)

^as

t → T

^{. Further, two}

independentinitial onditionsareneededfor(2)tobewelldened. Inordertoobtainanumeri alsolutionof

u(x, t)

^,

equation(2) isdis retizedsu hthat

u i (t)

approximates

u(x i , t)

^and

µ i (t) = µ(x i , t), η i (t) = η(x i , t)

^for

i = 1, . . . , n

^.

Thedis retizationismadeherewithnitedieren es,butitispossibletousebasisfun tions,niteelements,orany

othermethodof dis retization.Afteranitedieren edis retizationwehavethefollowingequations

F i (u 1 . . . , u n ) = η i u ¨ i + µ i ˙u i , i = 1, . . . , n

⁽³⁾

orresponding to equation (2). This dis retized se ond order dynami al systemis a systemof ordinarydierential

equations. We all theapproa h for solving (1) using (3) the Dynami al Fun tional Parti le Method, DFPM. We

emphasize that the idea behind DFPM is quite versatile. There is for example no linearrestri tion built into the

methodsonon-linearequations(fun tionals)are ertainlypossibletoatta k. However,inthisworkweshallinvestigate

onlythespe ial,but important, aseofsystemsoflinearequations.

A. Conne tion betweenDFPMand lassi al me hani s

Itis learthatDFPMisaparti lemethodbelongingtothedomainof lassi alme hani s. Anaturalinterpretation

ofthefun tional

F(u)

^{isthat after}dis retization itmaybeviewedasave torfor e eld. Considerthe spe ial ase whenthedis retizedversionofalineardierentialequationredu esto

n

^{linearequations,i.e.,}

Ax = b

^{. One}possibility isthentowritethefun tionalas

F(x) = b − Ax

^{. Thisfor eeldisthus}

n

-dimensional. Adynami alparti lemethod anbe onstru tedbyviewingthisastheproblem to determinethepositionsof

n

^{parti leswhere thefor e iszero,}

i.e.,theequilibriumpointofa onservativemany-parti lesystem. Thatis,if

F(x)

^is onservativeweknowthatthere exists a orrespondings alar potential

Φ (x)

^{. The} many-parti le system will tend to equilibrate towardsone of its minimumpointsifdissipationispresent[27℄. Inthisparti ular asetheDFPMequationis simplygivenby

b − Ax = η¨x + µ ˙x

wherethedissipationtermis

µ ˙x

^{. Inthe aseofasymmetri matrix}

A

^,amany-parti lepotential

Φ (x)

^doesexistand

isexpli itlygivenby

Φ (x) = 1 2

X

i

A ii x ² _i + X

i<j

A ij x i x j − X

i

b i x i = 1

2 x ^T Ax − b ^T x

⁽⁴⁾

be auseea h omponent

k

^of

F(x)

^isgivenby

F k = − ∂Φ

∂x k

= b k − A kk x k − X

i6=k

A ki x i

whi his ompletely onsistentwithea h omponentof

F(x) = b − Ax

^{. At thesolution}

x

^{, allthe omponentfor es}

fulll

F ^k = 0

^{andthusalsothegradientof}

Φ (x)

^{iszero. This riti alpoint} orrespondstoaminimumif

A

^ispositive

deniteand amaximumif

A

^{isnegativedenite. Obviously,if}

A

^{isnegativedeniteone aninsteadusetheansatz}

F(x) = Ax − b

^{. However,if}

A

^isnonsymmetri ,apotentialfun tion doesnot exist[17℄, soapossible onvergen eof thedynami alsystemneedstobeanalyzedinfurther detail. Also,the onvergen erateisnotavailablebytheabove

potentialanalysis. Theseimportantissueswillbeaddressedin thefollowing.

III. DFPMFORTHE

Ax = b

^PROBLEM

Considerthefollowingdis retizedintera tionfun tional

F = b − Ax

^,where

A ∈ R ^n×n

^and

b ∈ R ⁿ

^{. Inthis ase}

F

orrespondtotheordinaryresidualinnumeri allinearalgebra. Tosimplifythepresentationwemaketheadditional

assumptionthattheeigenvalues

λ (A)

^{areallrealandpositive. Other aseswillbedealtwithaswegoon. Topro eed,}

onsiderthese ondorderdynami alsystem:

F(x ¹ , x 2 , . . . , x n ) = b − Ax = η¨x + µ ˙x

⁽⁵⁾

where

µ ˙x

^isthedissipationterm. Thereals alarparameters

η > 0

^(mass)and

µ > 0

^{(damping)arehereassumedto}

be onstants. Wethussuggestaparti lemethodwherethepositions

x 1 , x 2 , . . . , x n

^{areoptimizedinordertoapproa h}

theuniqueequilibriumpoint,i.e.,

F(x 1 , x 2 , . . . , x n ) → 0

^as

t → T

^{. Inpra ti e,onemaystartwitharandomansatz}

fortheve tor

x (0)

^andlet

˙x (0) = 0

^{. Theparti ularnumeri altime}integration,usingarealvaluedtimestep

∆t > 0

tomovetheve tors

x (t)

^and

˙x (t)

^{,isheremadebythe ostee tiveandstablesymple ti Euler[2830℄:}

(

˙x (t + ∆t) = ˙x (t) + 

1

η b − ¹ η Ax (t) − ^µ η ˙x (t) 

∆t x (t + ∆t) = x (t) + ˙x (t + ∆t) ∆t.

(6)

Firstly,itisimportanttoestablishthata onvergentresultsolvestheoriginalproblem

Ax = b

^{. If}

x (t + ∆t) ≈ x (t)

^,

the se ond equation of (6) gives that

˙x (t + ∆t) ≈ 0

^{. If}

˙x (t + ∆t) ≈ ˙x (t)

^{, the rst equation of (6) gives that}

µ ˙x (t) ≈ b − Ax (t)

^{. Sin e}

˙x (t) ≈ ˙x (t + ∆t)

^and

˙x (t + ∆t) ≈ 0

^{,weindeedhavethat}

Ax (t) ≈ b

^.

Admittedly, the algorithm (6) is not parti ularly a urate but its stability and a ura y is still superior to the

Eulermethod [29℄. Further,it shouldbenotedthathigh numeri ala ura yofthemany-parti lesystemduringits

evolutiontowardsthestationary stateis not ne essary. Theonlydesiredpropertyisthat theapproa htowardsthe

equilibriumpointismadeas fastaspossible.

Thefeatures ofinterestarethus: 1)goodnumeri alstabilityallowingalargetime step

∆t

^{,and 2)}appli ation of a heapalgorithminordertospendlittleCPUtimeperiteration. Thesefeatures areallprovidedbythesymple ti

integrationmethod[28,30,31℄. Othersymple ti algorithmsmaystillbeofinterest,amongwhi htheVerlet-Störmer

methodisoneofthemost ommon[32℄. However,althoughthis methodismorea urate,itwasfoundthatnogain

ouldbeobtained. Thenumberofiterationstowards onvergen eissimilar(ifnotthesame)butthe ostperiteration

ishigher.

Inorder tobeexibleintheanalysisbelow,bothparameters

η

^and

µ

^{arekeptthroughoutthe}derivations,i.e.,we shallnotlete.g.

η = 1

^{as suggestedbya ontinuumview. Thereasonforthisis that,possibly,adis retenumeri al}

algorithmthat isnonlinear(as forexamplehigherorder symple ti algorithmsorRunge-Kuttamethods) anresult

inperforman epropertiesdependingonboth

η

^and

µ

^{. Su hanexampleisprovidedlater(A.58). Belowwemakethe}

assumptionthatthematrix

A

^isdiagonalizablesin eitsimpliestheanalysistremendously.

A. Optimalparametersfor asingle os illator

Inorderto identifyoptimalparameters

µ, η

^and

∆t

^{in (6)letusrewrite(5)bymakinga hangeof basisinto}

A

^'s

eigenve tors. Thatis,

b → c, x → u, A → Λ ^A =

^diag

(λ 1 , λ 2 , ..., λ n )

^{. Wethenobtainthefollowingde oupledsystem}

ofequations:

c − Λ A u = η¨ u + µ ˙u.

⁽⁷⁾

Aftera hangeofvariables

w i = u i − c i /λ i

^wehave

−λ ¹ w 1 = η ¨ w 1 + µ ˙ w 1

−λ ² w 2 = η ¨ w 2 + µ ˙ w 2

.

.

. .

.

.

.

.

.

−λ ⁿ w n = η ¨ w n + µ ˙ w n .

(8)

Thereadershould notethat (8)is onlyusedfortheanalysis ofoptimalparameters,i.e.,(8) should notbe onfused

withthea tual omputationalalgorithm(6). Thegreatadvantageof(8)isthatthewholeproblemhasbeenredu ed

to theanalysis of

n

^damped os illatorsin lassi alme hani s. Now onsiderthe symple ti timeintegrationof one su hos illator,

−λw = η ¨ w + µ ˙ w

^:

(

˙

w (t + ∆t) = ˙ w (t) + 

− ^λ η w (t) − ^µ η w (t) ˙ 

∆t w (t + ∆t) = w (t) + ˙ w (t + ∆t) ∆t.

(9)

Theiterations an onvenientlybesummarizedontheform

z n+1 = Bz n

^:

z n+1 =  w n+1

˙ w n+1



=

"

1 − ^λ _η ∆t ² 

1 − ^µ _η ∆t 

∆t

− ^λ _η ∆t 1 − ^µ _η ∆t

#  w n

˙ w n



.

⁽¹⁰⁾

One analso write

z n+1 = B ⁿ z 0

^{. If weapply the similarity} transformationsu h that

B = C ⁻¹ Λ B C

^where

Λ B =

diag

(α 1 , α 2 )

^,wenotethat onvergen eisa hievedif

|α ^1,2 | < 1.

^Aminimization ofthemaximumeigenvaluegivesthe optimal onvergen erate. Thetwoeigenvaluesareexpli itlygivenby:

α 1,2 = 1 − λ

2η ∆t ² − µ

2η ∆t ± ∆t

2η pζ

⁽¹¹⁾

where

ζ = (µ + λ∆t) ² − 4ηλ

^{. There are three possible ases, but only twoare of interest. The overdamped ase}

(

ζ > 0

^{) is well known in me hani s [33℄ to onlyslowlyrestore thesystem tothe} equilibrium point and istherefore dis arded. Ideally,wehavethe riti allydamped ase(

ζ = 0

^{), but we annot expe t allthe}os illatorsin (8) to be riti allydampedsotheunderdamped ase(

ζ < 0

^{)willalsobeofinterest. Theproblemisthentominimize}

ζ

^forall

ofthem. Thebest solutionforasingleos illatorisgivenby

(µ + λ∆t) ² = 4ηλ

^{. Sin e}

λ > 0

^wendthat

µ + λ∆t = 2pηλ

⁽¹²⁾

In this ase, both eigenvalues are simply given by

α 1,2 = 1 − pλ/η∆t

^so

∆t opt = pη/λ

^and

∆t max = 2pη/λ

^.

The optimal hoi e of time step,

∆t opt

^{, makesthe algorithm to nish in just one iteration. This is great for this}

singleos illator. However,the problem (5) orrespondsto manyos illators asgiven by(8), somostos illators will

needmu hmorethanonesingleiterationto onverge. Inthe aseofunderdampingone ansimilarlyasaboveshow

that

|α 1,2 | ² = 1 − (µ/η) ∆t

^{, so an} underdamped os illator onverges if

0 < (µ/η) ∆t < 1

^{. We are now prepared}

to understand the generalbehaviorof

n

os illators. Eq. (12) annot befullled forall theos illators but weshall attemptto make(12)tobenearlyfullledforthem.

B. Optimalparametersfor

n

os illators

As mentionedabove,it is important foroptimal onvergen ethat none ofthe os illatorsis overdamped, be ause

su h asituation would giveaveryslow progresstowardsthe equilibrium point. The onvergen eproperties ofthe

whole system would then suer sin e all the other os illators would have to wait. We shall therefore pro eed by

studyingthe aseswhere only riti alandlightdampingareallowed.

First onsideronlytwoos illators. As seenin (6)theparameters

µ

^and

∆t

^{areglobalforall}os illators. Eq. (12) thusgivesthat

 1 λ 1

1 λ 2

  µ

∆t



= 2 √ η

 √ λ 1

√ λ 2



.

⁽¹³⁾

If

λ 1 6= λ 2

^{,thereisanuniquesolutionfor}

µ

^and

∆t

^sothatbothos illatorsevolvea ordingtothe riti allydamped ase. Thissolutionisgivenby

µ = 2√η √

λ 1 λ 2 / √ λ 1 + √

λ 2

and

∆t = 2√η/ √ λ 1 + √

λ 2

. Inthegeneral aseof

n

os illatorsthe ondition

µ + λ i ∆t ≤ 2 √

ηλ i

^{mustbefullledforall}

λ i

^{. This onditiononlyholdsif}

µ opt = 2 √ η p

λ min λ max / p

λ min + p λ max



(14)

∆t opt = 2 √ η/ p

λ min + p λ max



(15)

This anbeseenbyinserting

µ opt

^and

∆t opt

^{into the ondition}

µ + λ i ∆t ≤ 2 √

ηλ i

^{. Onethenndsthat}

 2 p λ i − p

λ min − p λ max

 2

≤ p

λ min − p λ max

 2

(16)