Twee technieken voor het oplossen van een stelsel niet-lineaire vergelijkingen

Citation for published version (APA):

Baaijens, F. P. T., & Brekelmans, W. A. M. (1985). Twee technieken voor het oplossen van een stelsel niet-lineaire vergelijkingen. (EUT report. WFW, vakgr. Fundamentele Werktuigbouwkunde; Vol. WFW-85.042). Technische Hogeschool Eindhoven.

Document status and date: Published: 01/01/1985

Document Version:

Publisher’s PDF, also known as Version of Record (includes final page, issue and volume numbers)

Please check the document version of this publication:

• A submitted manuscript is the version of the article upon submission and before peer-review. There can be important differences between the submitted version and the official published version of record. People interested in the research are advised to contact the author for the final version of the publication, or visit the DOI to the publisher's website.

• The final author version and the galley proof are versions of the publication after peer review.

• The final published version features the final layout of the paper including the volume, issue and page numbers.

Link to publication

General rights

Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain

• You may freely distribute the URL identifying the publication in the public portal.

If the publication is distributed under the terms of Article 25fa of the Dutch Copyright Act, indicated by the “Taverne” license above, please follow below link for the End User Agreement:

www.tue.nl/taverne Take down policy

If you believe that this document breaches copyright please contact us at: openaccess@tue.nl

providing details and we will investigate your claim.

Frank Baaijens. Marcel Brekelmans.

Aug. 1985.

ABSTRACT

non-linear equations by F.P.T. Baaijens

W.A.M. Brekelmans

Two iterative procedures for the solution of a set of non-linear normal equations are discussed. Both methods are based on the successive determination and evaluation of improved estimate solutions.

TIle Newton methods or modifications of them are first considered, including the merits and limitations. For a number of non-linear problems the quasi-Newton methods do clearly have advantages. The fundamentals of these methods are explained and from a general frame specific choices for updating formulas, the Huang-class and the Broyden-class, are made. A special case of the latter results in the well-known BFGS-method. Attention is paid to the numerical application of this method.

O. Inleiding.

Twee technieken voor het oplossen van een stelael niet-lineaire vergelijkingen.

In dit rapport bespreken we een tweetal iteratieve technieken waarmee een oplossing van een stelsel niet-lineaire vergelijkingen van de vorm

f:DcRn~Rn; f(x)=O (0,1)

~ -v... ...

*'

kan worden gevonden. Laat x een oplossing en ~k een schatting VOOI die

oplossing aan het begin van de k-de iteratie zijn Seide technieken doen een uitspraak over hoe bij een gegeven schatting ~k een richting bepaald kan

*'.

.

.

worden waarlangs mogelijk een betere schatting voor x te v1nden lS. Die

nieuwe schatting is te schrijven als

(0.2)

waarbij a

k een stapgrootteparameter en ~k de zoekrichting is. Naast de ma-nier waarop ~k te bepalen is willen we criteria formuleren op basis waarvan we een uitspraak kunnen doen over de kwaliteit van ~k en over de grootte van Ilk' Dit blijkt aIleen onder bepaalde omstandigheden mogelijk te zijn, Deze omstandigheden hangen nauw samen met het a1 dan niet mogelijk zijn om de kolom f(x) op te vatten als de gradi~nt kolom van een scalaire functie. De mogelijkheden daartoe worden voor een belangrijk deel bepaald door een

stelling uit de minimaliseringstheorie.

Naast de klassieke Newton methoden bespreken we de BFGS methode die de laatste jaren steeds veelvuldiger wordt toegepast bij het oplossen van stel-sels niet-l.ineaire vergelijkingen. De belangrijkste redenen hiervoor zijn dat deze techniek een groter convergentiedomein heeft dan de Newton methoden (d.w.z. zij levert 'betere' zoekrichtingen af) en dat zij vaak aanz,ienlijk

effici~nter is. V~~r een goed begrip van de BFGS methode is het noodzake-lijk een aantal teehnieken uit de minimaliseringstheorie te bespreken.

1. Stelling.

Onder bepaalde voorwaarden kan het oplossen van f(x)=O gezien worden a1s het minimaliseren van een object funetie F(x). Een belangrijke situatie waarin

'"

dit mogelijk is wordt besehreven door de volgende stelling (Orthega en Rheinbolt (1910»

laat f:D C Rn"Rn met D een open convex gebied. Veronderstel dat f (x)

tenminste tweemaal continu diffentieerbaar is op D. Dan en sleehts dan als Vf

_....

T symmetrisch is voor aUe xeD, is f(x) gelijk aan de

gradi~nt

kolom g(x) . van een funetie F(x), waarbij f(x)=g(-x)=O een stationair

'"

..

punt van F(x) beschrijft .

_..

Veronderstel dat f voldoet aan de voorwaarden van deze stelling, dan kan op basis van die stelling de funetie Fex) in het algemeen niet bepaald worden. Dit is geen bezwaar omdat de in dit rapport te bespreken technieken geen gebruik makenvan de funetie F(x) .

_..

In het algemeen za1 g(x)=O geen eenduidige oplossing hebben, noeh zal de oplossing altijd een minimum van F beschdjven. Of een oplossing, zeg een loeaal minimum van F vastlegt kan worden geeontroleerd met behulp

*

x ,

...

van de Hessiaan van F: G., deze moet, als G. regulier is, in een minimum positief definiet zijn

T

*

y Q(x )y>O (1. 1 )

Vaak voldoet een stelsel vergelijkingen niet aan de eisen van de stelling. In dat geval kan het oplossen van f(x)=O toeh worden geformuleerd

_..

als een minimaliseringsprobleem door de funetie

(1.2)

te beschouwen. Eenvoudig kan worden nagegaan dat de globale minima van F(x) samenvallen met de oplossingen van f(x)=O. Het probleem is echter dat er locale minimaliseerders van F(x) bestaan die geen oplossing zijn van f(x)=O.

~ ~

Bovendien is het berekenen van de Hessiaan van F vaak uitermate gecompliceerd omdat daarin tweede orde afgeleiden van f een rol spelen. Deze problemen maken het onaantrekkelijk om het minimaliseren van F als al-ternatief te gebruiken voor het oplossen van f(x)=O.

~

2. Minimaliseringstechnieken.

Het doel van dit hoofdstuk is de introductie van een tweetal iteratieve technieken waarmee de positie van een Iocaal minimum (een minimaliseerder) kan worden bepaald: de Newton method en en de BFGS methode. Deze laatste methode val t binnen de klasse van de quasi -Newton' methoden, ook weI variabele met.rische methoden genoemd.

Dit hoofdstuk is als voigt opgebouwd. Ais eerste worden een viertal begrippen getntroduceerd die bij de toepassing van beide methoden een belangrijke rol spelen: de dal.i.ngseigenschap en daarmee gekoppeld de dalingsrichting, het lijnzoeken en de quadratische beeindiging. Daarna wor-den de werking en de voor en nadelen van de Newton methode besproken. De beperkingen van de Newton methode zijn aanleiding geweest andere mini-maliseringstechnieken toe te passen. Een belangrijke klasse van

alternatieven wordt gevormd door de quasi-Newton methoden waarvan we in dit rapport de BFGS methode bespreken.

Kenmerkend voor beide klassen van methoden is dat zij bij een gegeven schatting van de minimaliseerder ~k een zoekrichting ~k bepalen waarmee een nieuwe schatting van de minimaliseerder voIgt uit

daarbij treedt de constante ok op als stapgrootteparameter.

2.1 Palinqseiqenschap.

Bij voorkeur cre~ren we een rij van schattingen {~k} z.d.d.

(2.2)

pit wordt de dalingseigenschap genoemd. Een kolom ~k wordt een dalingsrich-ting genoemd indien voor voldoende kleine uk>O ~k+1=~k+ok~k voldoet aan de dalingseigenschap (2.2). Laat 0k~' dan voIgt uit de Taylorreeks

ontwikkel-ing van F(~k+1) rond ~k

(2.3)

dat bij 0k>O de kolom ~k een dalingsrichting is indien

(2.4)

Indien ~k=~ en ~k de functie F niet minimaliseert moeten hogere orde af-geleiden van F bestudeerd worden, hierop gaan we niet in.

2.2 Lijnzoeken.

Als ~k een dalingsrichting is zal niet voor iedere 0k>O voldaan worden aan de dalingseigenschap. Pit kan wel worden bewerkstelIigd door ok zodanig te kiezen d!t F minimaal is lan~s ~k'. Laat deze waarde van ok aangegeven worden door ok' dan moet voor Qk=a_k gelden

*

Het zoeken van a_k noemt men exact lijnzoeken. Uitwerken van (2.5) levert

*

dat voor ak=a_k moet gelden

T

~k+1~k=O (2.6)

In de praktijk wordt exact lijnzoeken zelden toegepast, meestal wordt a_k zodanig gekozen dat

T T

I ~k+1 ~k I (qll ~k~k I (2.7)

waarbij qlE;(O,1). Relatie (2.7) eist dat de helling van de kromme F(~k+ak~k)

in het punt ~k+1 voldoende moet zijn gedaald t.o.V. de helling in het punt

~k' zie figuur 2.1.

F

T T

I ~k+1 ~k I (!pI ~k~k I

Figuur 2.1 Niet exact lijnzoeken.

2.3 Quadratische be~indiginq.

Indien tijdens een iteratief pruces een quadratische functie in een eindig aantal stappen exact wOldt geminimaliseerd, dan noemen we dit quad-ratische beMndiging. Dat een algoritme deze eigenschap heeft is van belang omdat een te minimaliseren functie zich in de buurt van een locaal minimum

bij benadering quadratisch gedraagt. Dit is eenvoudig in te zien door de Taylorreeks ontwikkeling na de quadratische term af te breken.

2.4 Newton methoden.

De Newton methoden bereiken quadratische be~indiging in een iteratie. We maken onderscheid tussen de Newton methode (ook wel de full Newton-Raphson methode genoemd) en de gemodificeerde Newton methoden.

Als uitgangspunt voor het bepalen van de zoekrichting ~k bij de Newton methoden dient de Taylorreeks ontwikkeling van F(~k+~k) rond de laatst berekende schatting ~k

(2.8)

Indien ~k regulier is wordt de zoekrichting bepaald door te eisen dat deze quadratische benaded ng stationair is

(2.9)

De zoekrichting is een dalingsrichting indien

(2.10)

Hieraan wOldt zeker voldaan indien Q,k positief definiet is en 2kf~.

Er zijn een aantal nadelen aan de Newton methode verbonden. Ten eerste zal

Gx

niet overal positief definiet zijn, zodat ~k niet altijd een dalingsrichting is. Verder kan de quadratische benadering zeer slecht zijn omdat deze alleen het locale verloop van F(.) in rekening brengt. In beide gevallen kan een stapgrootte Q

k=1 slechte resultaten opleveren. Daarom wordt de Newton methode vaak in combinatie met een lijnzoekalgoritme toegepast. Daarbij moeten, indien ~k geen dalingsrichting is, n~gatieve

T

altijd toegepast kunnen worden: als ~k~k=O dan is F(.) stationair langs de lijn ~k+ak~k voor ak=O.

De Newton methode is verder een tamelijk dure methode. In iedere iteratie slag moet opnieuw de Hessiaan Yk worden bepaald, hetgeen vooral bij een groot aantal onbekenden veel tijd kost.

Om een aantal van de. nadelen te ondervangen zijn de gemodificeerde Newton method en qelntroduceerd. Daarbij wordt Yk vervanqen door een positief definiete matrix

a

k die op de een of andere manier een goede benadering is voor Yk; bijvoorbeeld ~k=Yj voor j<k waarbij ~k qedurende een aantal iteraties konstant wordt gekozen. Weliswaar levert deze werkwijze een traqer converqerend proces op maar vaak zijn de totale kosten lager dan bij toepassing van de Newton methode.

2.5 Quasi-Newton methoden.

Het belangrijkste probleem van aIle Newton methoden is dat zij een relatief klein converqentie domein hebben: slec~ts startwaarden ~O in een klein gebied rond een locale minimaliseerder x leveren een converqerend proces op. Daarom is er qezocht naar methoden die enerzijds robuuster zijn en anderzijds economischer zijn. Oit heeft geleid tot bijvoorbeeld de qeconjuqeerde qradi~nt methoden en de quasi-Newton methoden. De eerste wor-den hier niet besproken omdat Z.lJ bij het oplossen van stelsels

niet-lineaire vergelijkinqen een ondergeschikte rol spelen. Van de quasi-Newton methoden beperken we ons tot de Z.g. BFGS-methode, genoemd naar

Broyden-Fletcher-Goldfarb-Shanno. Deze werkwijze speelt de laatste jaren een steeds grotere rol bij het oplossen van stelsels niet-lineaire vergelijkingen.

Alvorens we deze methode bespreken moeten een aantal condities gefntroduceerd worden waarop de methode qebaseerd is, dat zijn achtereenvolqens: de quasi-Newton conditie, de geconjugeerdheidsconditie en daarmee samenhanqend de de erfelijkheidsconditie (eng.: heredity condition). De algoritmische structuur van de hier te bespreken quasi-Newton methoden is als voIqt.

k=1

*

Initialiseer bij een gegeven beginschatting :1 voor de minimaliseerder x de matrix

B

1. Voor B1 kan men bijvoorbeeld

I

of

~;1

kiezen.

while geen convergentie do

begin Stap 1: Bepaal een zoekrichting ~k volgens

(2.11)

b · . uT b d " -1 waar IJ ~k een ena erlng 15 van ~k .

Stap 2: De nieuwe schatting

(2.12)

wordt gevonden d.m.v. een lijnzoek algori.tme, waarbij a k zodanig wordt bepaald dat

T T

I

~k+1 ~k

I

<tp

I

~k~k

I

",e (O, 1 ) (2.13)

Stap 3: Daarna wordt Rk geherwaarderd (updated) door toepassing van

k=k+1

(2.14)

waarbij ~k wordt berekend met behulp van :k' :k+1' 2k' 2k+1 en

B

_{k ·}

2.6 Quasi-Newton conditie.

Veronderstel dat tot iteratie k aJle gegevens bekend Zl.JD. In

iteratiestap k wOldt ~k+1 berekend en vervolgens willen we dat

B

_{k+1 waarmee} de nieuwe zoekrichting Pk+1 wordt bepaald op de een of andere manier de krommingseigenschappen van de. object functie F(.) representeert. Een manier om dit te realiseren is door te voldoen aan de quasi-Newton conditie. Deze is als voIgt gedefinieerd. I.aat

dan voldoet B_{k+1 aan de Z.g. quasi-Newton conditie indien}

(2.17)

Hx+1 bevat dezelfde krommingsinformatie als de quadratische benadering van F ( .) 1angs Pk'

2.7 Gecon'juqeerdheidsconditielerfeHjkheidscondi tie.

Een set van kolommen ~i wordt onderling geconjugeerd genoemd met betrekking tot een positief definiete matrix y dan en slechts dan a1s

T

p.Q.p.=O

... 1 ... J

v

ifj

Hierbij wordt verondersteld dat Pif~ voor alle i.

(2.18)

Bewezen kan worden dat bij het sequentieel toepassen van geconjugeerde zoekrichtingen in combinatie met exact lijnzoeken, de minimaliseerder van een quadratische functie met positief definiete Hessiaan

y

in ten hQogste n (n is het aantal onbekenden) stappen wordt gevonden. In dit bewijs wordt

ondermeer gebruik gemaakt van het feit dat geconjugeerde kolommen onder ling onafhankelijk zijn. Het bewijs van deze eigenschappen wordt gegeven in

ap-pendix A.

Een belangrijk gevolg van het toepassen van geconjugeerde zoekrichtin-gen in combinatie met exact lijnzoeken bij het minimaliseren van een quadratische functie is dat .

v

j<k+1 (2.19)

Dit resultaat is te vinden in appendix A.

Met betrekking tot de matrices H. van de quasi-Newton method en eisen we

J

dat Zl) zoekrichtingen genereren die geconjugeerd zijn indien de methode wordt toegepast bij het minimaliseren van een quadratische functie met Hessiaan g waarbij exact lijnzoeken wordt gebruikt; m.a.w. in iteratiestap k eisen we dat Hx+1 een zoekrichting ~k+1 zal opleveren die-voldoet aan

v

j(k+1 (2.20)

T

Met ~k+1=-Rk+1~k+1 levert (2.20) _{dat Rk+1} moet voldoen aan

v

j<k+1 (2.21)

VerondersteI dat aIle voorgaande zoekrichtingen p., j<k+1, onderIing

gecon-. ~J

jugeerd zijn en dat exact lijnzoeken is toegepast, dan geldt

v

j<k+1 (2.22)

Hieruit blijkt dat aan (2.21) voldaan kan worden door

R

_k+₁ zodanig te kiezen dat

iL·+₁GP .=gp.

II. ~J ~J

v

j(k+1

(2.23)

Met a.p.=Ax. en a.Qp.=6g. kan (2.23) ook worden geschreven als

H."+ 1 Ag . =pAx .

A .. J ",J

v

j<k+1

(2.24)

Conclusie: door H.

k+1 aan (2.24) te laten voldoen worden geconjugeerde zoek-richtingen gegenereerd. Merk op dat voor j=k relatie (2.24) juist de quasi-Newton conditie representeert indien 0k=P wordt gekozen.

Naast de eigenschap da~ relatie (2.24) zorg draagt voor het genereren van geconjugeerde zoekrichtingen bij het minimaliseren van een quadratische objectfunctie kan zij in het algemeen op een andere manier worden geinterpreteerd. De eia (2.24) zorgt er namelijk voor dat de krommingsin-formatie uit vroegere iteraties wordt doorgegeven naar latere iteraties. Daarom wordt relatie (2.24) ook weI de erfelijkheidsconditie genoemd.

Een andere manier (naast de manier uit appendix A) om de quadratische

be~indigingseigenschap te laten zien iridien

H

k+1 voldoet aan de er-felijkheidsconditie is als voIgt. Beschouw het bijzondere geval

H +llJ.g. =pAx.

"""'ll .. J .. J

v

j<n+1

uit Q.A~k=A~k voor een quadratische functie voIgt

Laat dan geldt H+1Q.AX ·=oAx. =n .. J .. J V j(n+1 (2.25) (2.26) (2.27) (2.28)

Ais de kolommen IJ.x. onderling geconj'ugeerd zijn met betrekking tot G. dan .. J

zijn zij onderling onafhankelijk waardoor l:.X inverteerbaar is. Dan geldt

Met 0=1 blijkt dat deze method en , bij het bepalen van een minimaliseer-del van een quadratische functie,

quadratische object functie F(.)

juist de inverse van de Hessiaan van de hebben berekend. Dit is een bijzonder aantrekkelijke eigenschap van deze methoden omdat iedere functie F(.) zich in de buurt van een minimum quadratisch gedraagt.

2.8 De Huang klasse van 'updating' formules.

De essenti~le stap bij de quasi-Newton methoden is het defini~ren van de updating matrix Q_k, Deze moet zodanig zijn dat H

k+1=Hk+Qk tenminste aan de quasi-Newton conditie (Hx+1A2k=OkA2k) voldoet. Deze conditie legt echter geen eenduidige eis op aan Qk' Bij de in deze paragraaf te bespreken Huang-klasse van updating formules worden de matrices Q

k zodanig gekozen dat: a} Q_k ten hoogsterang 2 heeft.

b) Bx+1(=Rx+Qk) voldoet aan de erfelijkheidsconditie. Hierdoor wordt be-werkstelligd dat bij het minimaliseren van een quadratische functie

geconjugeerde zoekrichtingen worden gegenereerd indien exact lijnzoeken wordt toegepast.

c) voor het bepalen van Q_k aIleen schattingen van de minimaliseerder, de gradienten in die schattingen en Hk nodig zijn.

d)

Ox

uitsluitend wordt bepaald met behulp van grootheden ult de vorige en de huidige iteratie.

Eigenschap b) is van groot belang omdat vrijwel iedere objectfunctie zich in de buurt van een locaal minimum bij benadering quadratisch gedraagt.

De BFGS updating formules voor Q

k kunnen uit de klasse van Huang updat-ing formules worden afgeleid.

Toepassen van de quasi-Newton conditie op relatie (2.14) levert

Hierbij is verondersteld dat Q_k gedurende het gehele iteratieproces constant wordt gekozen, i.e. Qk=O. Aan (2.30) kan worden voldaan door Q_k te kiezen

volgens

(2.31)

T T

met !kA~k~O en :kA~k~O' In deze relatie spelen nog maar twee onbekende

kolommen een rol: !k en ~k' Door deze keuze wordt tevens voldaan aan er-felijkheidsconditie voor j=k.

Relatie (2.31) vormt het uitgangspunt voor de Huang klasse van 'updating' formules. Door het in rekening brengen van de

erfelijkheidscon-.

ditie en door gebruik te maken van het feit dat daardoor, onder bepaalde voorwaarden, geconjugeerde zoekrichtingen worden gegenereerd kunnen een aan-tal el.sen geformuleerd worden waaraan !k en ~k moe ten voldoen.

Opdat .Q_k volgens (2.31) een zodanige vorm heeft dat

H

_k+₁ voldoet aan de erfelijkheidsconditie (H:.'+1Ag'=QAX. voor alle j<k+1) volStaat de eis

II. .. J .. J

QkAg.=gAx·-RkAx . _{.. J .. J .. J}

v

j<k

Er van uitgaande dat Hk aan de erfelijkheidsconditie voldoet H:..Ag. =OAx .

II. .. J .. J

v

j<k

voIgt dus als extra voorwaarde aan .Q

k volgens (2.31)

v

j<k

Aan deze eis wordt ondermeer voldaan indien

T T ykAg'=ZkAg.=O _{.. ..J.. ..J}

v

j<k (2.32) (2.33) (2.34) (2.35)

Een geschikte keuze voor!k en ~k kan gemaakt worden op basis van de volgende beschouwing. Ooordat de voorgaande, bekende, zoekrichtingen

Pj

met j(k+1 onder ling geconjugeerd zijn indien een quadratische functie met Hessiaan ~ wordt geminimaliseerd waarbij exacte lijnminimalisering is

T T

toegepast, geldt a) ~kG~j=O voor aIle j<k en b) ~k+l~j=O voor aIle j<k+1. De eigenschap a} Ievert m.b.v. (2.15) en (2.16)

TI'!1"\ 0

>

Tr... T T 0 V J' <k

Pk~'= _{.. ..)} == akPk~'p, = _{.. .) .. )} AXkC~·P. _{.. ) ... J} = AXkAg. = _...) (2.36)

Nerk op dat als aan AXkTA9 .=0 voor alle j<k voldaan wordt, dat dit niet ... ..J

betekent dat de zoekrichting ~k geconjugeerd is ten opzichte van aIle voor-gaande zoekrichtingen. Oit geldt aIleen indien de te minimaliseren functie quadratisch is en indien daarbij exact lijnzoeken wordt toegepast.

Uit de eigenschap b) voIgt

v

j<k (2.31)

Omdat volgens (2.23) geldt PPj=[k+l~j voor aIle j<k+1 levert (2.37)

waarbij gebruik is gemaakt van (2.34).

Uit de re]aties (2.37) en (2.38) blijkt dat aan (2.35) voldaan kan wor-den door !k en:k als een Hneaire combinatie van

A~k

en HiA?k te kiezen. Met deze keuze is de meest algemene vorm van Huang klasse van 'updating' formules te schrijven als

(2.39)

waarbij ~k' ~k' 'k en w_k willekeurige constanten zijn zodanig dat de noemers ongelijk aan nul zijn.

2.9 De Sroyden klasse van 'updating' formules.

Een probleem van de Huang klasse is dat de matrices

B.

niet

nood-- ]

zakelijk symmetrisch zijn ook al is de initi~le matrix

BO

weI symmetrisch. Daarnaast zijn de matrices

B.

niet per definitie positief definiet. Deze

J

twee eigenschappen zijn bijzonder prettig bij de numerieke toepassing van de algoritmen. Vooral is te verwachten dat het positief definiet zijn van de benadering van de Hessiaan, of de inverse ervan, prettig is, immers ter plaatse van een minimum is Hessiaan positief definiet. Daarnaast levert een positief definiete matrix een dalingsrichting als zoekrichting at, hetgeen aangenaam is. Door een normalisering van de kolommen ~k en ~k kunnen twee van de onbekenden CPk' w_k' ~k en !lIk get!li~inieerd worden. Aan (2.31) is te zien dat deze normaliser.ing het resultaat voor Q.

k niet beinvloedt. Door een keuze te maken t.a.v. de twee resterende onbekenden kan de Broyden klasse van 'updating' formules worden afgeleid. Deze levert symmetrische benaderingen voor de Hessiaan (of de inverse daarvan) af. Een bijzonder geval van die klasse vormt de BFGS 'updating' tarmule. - Van deze kan onder bepaalde omstandigheden bewezen worden dat zij positief definiete benaderin-gen oplevert.

De procedure is als voIgt: normaliseer de kolommen !k en :k volgens (2.40) Dan kunnen CPk en w

k get!limineerd worden volgens

(2.41)

Substitutie van de bovenstaande vergelijkingen in (2.39) levert met de keuze J.ik=-"'k een symmetrische matrix op. De zo ontstane vergelijking vormt de Broyden klasse. Als onbekende parameter komt nog "'k v~~r.

De BFGS 'updating' formu]e kan gevonden worden door

"'k6~~62k=1

te kiezen. lndien lijnzoeken wordt toegepast volgens (2.17) en indien de

initi~le

waarde van!ik

posi~ief

definiet is zal

6~~62k)O

(zie appendix B).

Dan kan bewezen worden dat deze keuze leidt tot positief definiete matrices

~, zie Scales [1] pp. 92-94. Na enig schrijfwerk is !i

k+1 dan te schrijven als

De voorwaarden waaronder !i k+1

definiete

Hx

worden ge,tllustreerd aan

[2] beschreven equivalente formulering

met

v _{_k} =

(2.43)

positief definiet is bij positief de hand van de door 'Matthies en Strang

voor B._k+₁, deze luidt

(2.44)

(2.45)

(2.46)

T

Laat Ax.=(l.+~k~k)' dan is onmiddellijk in te zien dat B._{k+1 semi-p,os.itief} definiet is als !i

(2.47)

Indien Ak regulier is zal

H

k+1 positief definiet zijn. We kunnen een uitspraak doen over de regulariteit door de eigenwaarden van Ak te bestuderen. Het blijkt dat N-1 eigenwaarden 1 zijn en ~en eigenwaarde is gelijk aan

1+!~~k

wat

uitgewe~kt

(2.48)

opleverl. Hieruit volgt direct dat Ak zeker regulier is als

A~iA~k>O.

2.10 Numerieke qebruik van de BFGS methode.

De relaties (2.43) of (2.44) worden vrijwel nooit direct toegepast voor het berekenen van li

k+1. Met name bij toepassing in eindj.ge elementen methode berekeningen zou dit een bijzonder onaantrekkelijke werkwijze zijn omdat dan geen enkel gebruik van de structuur van de stijfheidsmatrix (die een overeenkomstige rol speelt a15 ~) zou worden gemaakt. Deze zou immers verloren gaan doordat b.v. de matrix A~kA~~ een volle matrix is. Daarom wordt over het algemeen voor een recursief schema gekozen waarmee een nieuwe zoekrichting wordt bepaald.

We geven hier een recursief schema waarmee een zoekrichting bepaald kan worden op basis van de updating formule van Matthies en Strang: relatie

(2.44).

Als eerste iteratie stap wordt vaak voar een standaard Newton stap -T gekozen. In essentie betekent dit dat als initi!le waarde van H1 de Hessiaan Y1=Y(~1) genomen. De zoekrichting ~1 wordt bepaald d.;.v. het oplossen van het stelsel Y1~1=-~1 via een directe methode (b.v. LDL decom-positie met terugsubstit.utie). De gedecomponeerde van ~1 wOldt bewaard. De

zoekrichtingen voor k>1 worden daarna bepaald door

~k=-li~2k

_{waarbij lik} 'ge-update' wordt m.b.v. (2.44). Zo wordt bijvoorbeeld ~J als voIgt berekend

In een recursief schema levert dit

Bereken d.m.v. terugsubstitutie : uit ~1:=~

T T T"

(hiermee is :=R1(I+~1~1) (1+~2~2)23 bepaald)

T fl=v1T. ... r=r+w1fl ~ '" . . trw ... T T T . T T (hiermee is :=(I+~2~2)(I+~1~1)R1(I+~1~1)(I+~2~2)23 bepaald)

Tenslotte volgt P3 uit P3=-!'

Het algemene recursieve schema ziet er ais voIgt uit. Voor k=1 wordt

~1 b;rekend door ~'~1=-21 op te lossen. Voor k>1 wordt" ~k bepaald m.b.v. ~k=-Hk2k' volgens

Initialiseer q=gk

_..

_..

for j=k-1 .§.lli -1 until 1 do

J,.."" .... .;... T

~ a=~j~

q=q+av.

.... ...J

end

bereken d. m . v. terugsubs ti tub e r uit 111 :=~

b . T ~ fI=v.r ~) ~ r=r+f:Jw . . . "'W . . ) De zoekrichting ~k is dan ~k=-:

V~~r de updating formule (2.43) is een soortgelijk recursief schema te formuleren.

3. ContinuOmsproblemen.

In dit hoofdstuk bespreken we de toepassingsmogelijkheden van de tech-nieken uit het voorgaande bij het oplossen van een stelsel niet-lineaire

vergelijkingen dat afkomstig is van een bepaalde klasse van

continuOmsproblemen.

We beperken ons tot quasi-statische problemen die beschreven kunnen worden door de volgende gewogen afwijkingen vergelijking

++ c ~ + + +

I

0: (vw) dQ = Jr.w dO + Jq.w dr (3.1)

Q Q

r

Na een geschikte discretisatie kan hieruit een stelsel vergelijkingen worden afgeleid van de vorm

waarbij

_Ei'

_Er

belasting volume T K =vf , a=i, f, q. ~ "" ... a (3.2)

en f de bijdragen van resp. de inwendige spanningen, de

.. q T

en de oppervlakte belasting representeren. Laat K=Vf en De matrix K wordt vaak de stijfheidsmatrix genoemd.

'*

De oplossing x=x van (3.2) willen we via een iteratief proces vinden dat een ana loge structuur heeft aan hetg~en in hoofdstuk 2 is gebruikt, dus:

'*

gegeven een zekere schatting ~k voor x willen we een zoekrichting ~k

'*

bepalen waarlangs we een betere schatting voor x verwachten volgens

De vragen die we ons stellen Zl)n: wat is een geschikte zoekrichting en hoe

.

*

bepalen we die, op grond waarvan vinden we de ene schatting voor x beter

_...

dan de andere en, direct daarmee samenhangend, wat is een geschikte keuze van uk' Ret zal blijken dat we in een aantal situaties nutt:ig gebruik kun-nen maken van de kennis rondom de minimaliseringstechnieken. Daarbij zal de stelling uit hoofdstuk 2 een.belangrijke rol spelen.

In het algemeen is ( niet symmetrisch en kan f dus niet gezien worden als de gradi~nt van een scalaire functie. In dit geval kunnen we geen gebruik maken van de stelling uit hoofdstuk 1: we kunnen het probleem f(x)=O ~ niet associ~ren met een minimaliseringsprobleem. Toch kunnen zowel de Newton als de quasi-Newton methoden worden toegepast. In beide gevallen wordt een nieuwe zoekrichting bepaald door f op de een of andere manier _... locaal door een lineair verband te benaderen: m.b.v. Newton wordt ~k bepaald volgens Kk~k=-!k en via de quasi-Newton methode wordt ~k bepaald d.m.v. via Pk=-I!Tkfk . Bij toepassing van de BFGS methode wordt als ihitHHe waarde voor

~ ... 1 T

Hx

vaak H(t2"(Ko+Ko) genomen. Er bestaan echter geen criteria op basis waar-van we ~k een goede of slechte richting kunnen noemen; het begrip dalingsrichting heeft hier b.v. geen enkele betekenis. Ook voor het kiezen van de stapgrootte uk ontbreekt een criterium. Vaak wordt het lijnzoeken toch toegepast met ais argument dat de component van !k+1 in de richting van ~k in absolute waarde voldoende gedaald moet zijn.

Echter, afhankelijk van de constitutieve vergelijking is Ie _-1 weI symmetrisch. Oaarnaast zijn Kf=Q en Kq=Q indien de belastingen geen functie van het verplaatsingsveld Zl)n. Als verwacht kan worden dat (=K_i positief

*

definiet is in de oplossing x=x dan kunnen de technieken uit het voorgaande zondermeer worden toegepast.

Indien het niet symmetrisch Zl)n van ( volledig wordt bepaald door

Kr

en Kq (dus

Ki

is. weI symmetrisch) en verwacht kan worden dat .K positief definiet is voor x=x dan kunnen we het begrip dalingsrichting en lijnzoeken ook met succes toepassen. Oit is als vo]gt in te zien. Vrijwel altijd wordt de belasting incrementeel aangebracht. Oaarbij is de verandering van de belasting als gevolg van de vormverandering van het lichaam gedur~nde een increment gering; m.a.w. zij beinvloedt het probleem in een geringe mate.

Dit betekent dat het oorsponkelijke probleem redelijk goed benaderd kan wor-den door het probleem waarbijde belasting· niet t.O.V. de huidige configuratie maar t.o.v. de laatst geschatte configuratie wordt genomen. Bij dit probleem behoort een symmetrische stijfheidsmatrix. ·Met behulp van dit probleem kunnen we een uitspraak do en over de kwaliteit van de zoek-richting en over de keuze van de stapgrootte. Omdat verwacht wordt dat

K

:t positief definiet is voor x=x als een minimaliseringsprobleem.

kunnen we dit aanverwante probleem opvatten Het oorspronkeli jke probleem gedraag.t zich naar verwachting op een analoge manier als dit minimaliseringsprobleem. Daarom is het gerechtvaardigd te eisen dat een zoekrichting een dalingsrichting is en tevens is het te rechtvaardigen dat lijnzoeken wordt

toegepast.

4. r.i teratuur.

[1] L. E. Scales. ' Introduction to nonlinear optimization I (1985) CYI,8SSCA

bsa.

[2] H. Matthies and G. Strang. 'The solution of nonlinear finite element equations't Int. Jrnl. Num. Meth. Engn. 14 (1919) 1613-1626.

[3] K.J. Bathe and A.P. Cimento. 'Some practical procedures for the solution of nonlinear finite element equations', Comp. Meth. Appl. Mech. Engn. 22

(1980) 59-85.

[4] H.Y. Huang. 'Unified aproach to quadratically convergent algori.thms for function minimization', Jrn. Opt. Theory. Appl. 5 (1970) 405-423.

Appendix A.

In deze appendix bewijzen we dat het sequentieel toepassen van gecon-jugeerde zoekrichtingen bij het minimaliseren van een quadratische functie quadratische be~indiging tot gevolg heeft.

Beschouw de functie F(x)-

_..

F(x) ₌ 1 T _{~ {i~+~ ~+c} T (A.1)

met positief definiete symmetrische {i. De gradient van F is gegeven door

g(x)

=

Qx+b (A.2)

IJaat ~k een zoekrichting. optreedt vaar

schatting voor de minimaliseerder van F en zij ~k een Eenvoud.ig kan worden nagegaan dat exacte lijnminimalisering

(A.3)

De gradienten in en een vorig punt x. zijn ais voIgt

.. J gerelateerd k I: a..{ip. . . 1 .. J. 1=J

v

j(k+1 (A.4)

Vermenigvuldig het getransponeerde van (A.4) met p. dan voIgt .. J lndien (A.5) p. met .. 1 i(k+1 k T

r

a.p.{ip. . . 1 .. 1 .. J l=J geconjugeerd , V j(k+1 (A.5)

is met p. met j<k+1,

_..

j:ji dan volgt uit

T T T

qk+1 P ·-q .P . _{~ .. J ~J .. J}

=

u.p.Qp. _{J .. ) .. )}

v

j(k+1

Met u_j voIqens (A.3) (exact lijnzoeken) resulteert dan

v

j<k+1

(A.6)

(A. 7)

Omdat de qeconjuqeerde zoekrichtinqen onderlinq onafhankelijk zijn (hetqeen eenvoudiq te bewijzen is) en omdat p.;O voor aIle

_..

j~n, vormen de kolommen

) ~

~1" "'~n een basis in Rn zodat de eniqe ko]om ~n+1 die voldoet aan (A.7)

Appendix B.

In deze appendix bewijzen we dat lijnzoeken volgens (2.17) bij een positief definiete Hessiaan !:i

k oplevert dat

6~~62k>O.

Irnmers "volgens (2.17) moet gelden dat

(B. 1 )

T T

Er zijn twee gevallen de onderscheiden: 1) 2k+1~k<O en 2) 2k+1~k>O.

Geval 1: (B.1) Ievert dan

(B. 2)

Hieruit voIgt direct

(B. 3)

of weI

(B.4) T

Omdat ~e(O,1) en 2k~k<O (volgens (2.4» levert (B.4) dat

(B.5)

Geval 2: uit

~~+1~k>O

en

~~~k<O

voIgt direct

(B.6)

T

gevallen hebben we gevonden dat 6~k~k>O. Omdat ak>O blijkt dat Q.E.D.