Simulatie + Contractie = Bewijs

(1)

komen, gewoonlijk teruggevallen op numerieke simulaties. Tegenwoordig kunnen zulke simulaties, dankzij hoogstaande al- goritmes, software en hardware vaak eenvoudig en snel ontwikkeld en ‘gedraaid’

worden. Echter, deze simulaties leveren in de meeste gevallen geen daadwerkelijk bewijs voor het bestaan van een oplossing of een garantie dat een numerieke methode een oplossing goed benadert.

Aan de basis van de methode die we in dit artikel beschrijven ligt het idee dat dit soort numerieke simulaties, ondanks deze beperkingen, toch gebruikt kunnen worden om een wiskundig exact bewijs te leveren. Om precies te zijn zullen wij laten zien dat in de buurt van een numeriek gesimuleerde ‘oplossing’ (dat wil zeggen, een uitkomst van een numerieke berekening waarvan we hopen/verwachten dat het een goede benadering van een oplossing is) een unieke wiskundig exacte oplossing bestaat. Sterker nog, onze methode vertelt ons expliciet hoe dicht de ware oplossing bij de gesimuleerde benadering ligt en hoe groot de omgeving is waarin deze uniek is.

De allereerste stap van de methode is om het probleem, of het nu afkomstig is uit bijvoorbeeld een stelsel gewone differentiaalvergelijkingen of bestaat uit een partiële differentiaalvergelijking met rand- voorwaarden, te herformuleren als

( ) ,

F x =0 (1)

waarbij F een afbeelding is tussen, laten we zeggen, twee vectorruimten. Die vector- ruimten zullen in het algemeen oneindig dimensionaal zijn, en dan kiezen we doorgaans Banachruimten, maar in eenvoudige gevallen kunnen ze ook eindig-dimensio- of enkele speciale oplossingen te vinden,

zoals periodieke oplossingen of verbindende banen. Door middel van zogenaamde forcing theorems kan het bestaan van dit soort speciale oplossingen gebruikt worden om stellingen te bewijzen die uitspraken doen over globale eigenschappen van het systeem, zoals de aanwezigheid van chaos of het bestaan van speciale families van oplossingen. Een beroemd voorbeeld is de stelling “periode drie impliceert chaos” voor iteraties in één dimensie. We komen daar later op terug.

De problemen in de dynamische systemen zijn eigenlijk zonder uitzondering niet-lineair. Kwalitatieve, grove beschrij- vingen van robuste eigenschappen van dit soort systemen kunnen vaak gegeven worden met behulp van topologische en variationele methoden. Bovendien kan in sommige gevallen een ‘kleine’ parameter uitgebuit worden om oplossingen in limiet- gevallen via asymptotische analyses te beschrijven. Maar door het niet-lineaire ka- rakter van de problemen kunnen bewijzen van stellingen die kwantitatieve uitspraken doen over dynamische systemen voor ‘gewone’ parameterwaardes in vrijwel geen enkel geval via analytische weg gevonden worden. Daarom wordt bij bestudering van dit soort problemen, die zowel in de wiskunde zelf als in toepassing erg vaak voor- Sinds het Manhattan project in de Tweede

Wereldoorlog, maar vooral vanaf de op- komst van de op transistors gebaseerde computers in de jaren zestig zijn computersimulaties niet meer weg te denken uit de wiskunde. In de toegepaste wiskunde (en ver daarbuiten) vormen simulaties vaak het eindstadium van een oplossings- strategie voor een gemodelleerd probleem.

Daarbij bestaat het model bijvoorbeeld uit een groot stelsel gekoppelde (differentiaal) vergelijkingen of een niet-lineaire partiële differentiaalvergelijking, zoals de Navier–

Stokes-vergelijking voor vloeistofstromin- gen. In de zuivere wiskunde dienen simulaties vooral om de wiskundige te helpen haar probleem te visualiseren en om inspi- ratie op te doen.

Een cruciale vraag die bij het gebruik van computersimulaties vrijwel altijd om de hoek komt kijken is hoe betrouwbaar deze resultaten zijn. Computer-assisted bewijzen zijn bij uitstek geschikt om die betrouwbaarheid onomstotelijk vast te stellen.

De problemen die wij met onze methode zullen behandelen zijn vrijwel allemaal afkomstig uit het vakgebied van de niet-lineaire dynamica, zoals het vinden van oplossingen van differentiaalvergelijkingen of banen van iteratieve afbeeldingen. In veel situaties kan het al lonen om slechts één

Simulatie + Contractie = Bewijs

In hoeverre is het verantwoord om de resultaten van een simulatie van een probleem als benadering van een oplossing van het probleem te interpreteren? Een van de manieren om deze vraag te benaderen is door middel van zogenaamde computer-assisted (door de computer ondersteunde) bewijsmethoden. In dit artikel presenteren Jan Bouwe van den Berg, Chris Groothedde en Ray Sheombarsing een methode die, met als uitgangspunt een zorgvuldig geproduceerde numerieke simulatie, de mogelijkheid biedt om het probleem op een wiskundig exacte manier op te lossen met behulp van een computer.

Jan Bouwe van den Berg

Afdeling Wiskunde Vrije Universiteit Amsterdam janbouwe.vanden.berg@vu.nl

Chris Groothedde

Afdeling Wiskunde Vrije Universiteit Amsterdam c.m.groothedde@vu.nl

Ray Sheombarsing

Afdeling Wiskunde Vrije Universiteit Amsterdam r.s.s.sheombarsing@vu.nl

(2)

optelling, vermenigvuldiging, et cetera). Dit betekent dat, zelfs in eindig veel dimensies, simulaties en numerieke berekenin- gen op basis van floating-point getallen afwijkingen met zich mee brengen. Om toch op een rigoureuze manier om te gaan met floats wordt gebruikt gemaakt van zogenaamde interval-aritmetiek. Interval- aritmetiek is een manier om, binnen een eindige verzameling van getallen F (zoals de floats), door stelselmatig naar boven én naar beneden af te ronden toch op ma- thematisch gegarandeerde resultaten uit te komen.

Als voorbeeld bekijken we de opera- tie vermenigvuldiging : R R) # "R. Door- gaans zal het niet zo zijn dat uit ,a b!F volgt dat ook a b) !F geldt. In plaats van de gewone vermenigvuldiging gebruiken processoren dus een andere, numerieke, vermenigvuldiging : F FU # "F. Deze numerieke vermenigvuldiging wordt door computerchips zo uitgevoerd dat deze de ware vermenigvuldiging benadert:

a b a) - Ub 1 e, waarbij e de ‘machi- ne precision’ wordt genoemd; gewoonlijk geldt dat e.10^-¹⁶. Met andere woorden, het getal a bU is een afronding van het ge- tal a b) . Het is echter a priori niet bekend of dit getal naar boven of naar beneden wordt afgerond, waardoor er nog meer informatie verloren gaat.

Door systematisch naar boven én naar beneden af te ronden, kunnen we dit soort verlies van informatie tegengaan. In plaats van te rekenen met losse floats, zullen we rekenen met intervallen met eindpunten in F:

{[ , ]a b R: ,a b F}.

I_F

r

^def= 1

r

!

r

Hierin is R

r

^def=R,{-3,+3} de uitgebreide reële lijn en zijn F

r

^def=F,{-3,+3}_de uitgebreide floating-point getallen. Voor elke x!R kunnen we dus een interval [ , ]a b !IF

r

vinden zodat x![ , ]a b.

Het idee is nu om een vermenigvuldiging : IU F

r

#IF

r

"IF

r

tussen intervallen te definiëren, zodanig dat voor elke I1 en I2 in IF

r

^{geldt dat}

en :

.

x x x I x I

I I R

I

1 2 1 1 2 2

1 2 ^F

) d d d

U d

1

r r

# -

In de praktijk kan de ondergrens van I₁UI₂ gevonden worden door al de uitein- den van I₁ en I₂ met elkaar te vermenig- vuldigen, deze naar beneden af te ronden (zodat we in F

r

belanden) en vervolgens het kleinste resultaat te kiezen. Analoog Net als bij de methode die we in dit ar-

tikel beschrijven berusten de bewijzen van bovengenoemde historische stellingen op het verifiëren van ongelijkheden. Om dit met een computer te kunnen doen wordt gebruik gemaakt van interval-aritmetiek.

Interval-aritmetiek is een methode waar- mee computersoftware op een geautoma- tiseerde manier afrondfouten kan bijhou- den, zodat ongelijkheden op een exacte manier gecontroleerd kunnen worden. Het belang van afrondfouten moet wel enigs- zins gerelativeerd worden. In de praktijk vallen ze in het niet in vergelijking met de analytische schattingen die vooral bij oneindig-dimensionale problemen de crux (en het struikelblok) voor het oplossen van het probleem zijn. Niettemin is het essen- tieel dat naast het deel van het bewijs dat door een brein van vlees en bloed is geconstrueerd, ook het deel van het bewijs dat aan een brein van silicium is uitbesteed foutloos is. Daarom besteden we in de volgende paragraaf kort aandacht aan de vraag hoe we door toepassing van interval-aritmetiek de computer kunnen inscha- kelen om ongelijkheden wiskundig exact te controleren.

De laatste jaren groeit de interesse in computer-assisted proofs in de dynami- sche systemen snel. Dit komt vooral doordat de ontwikkelingen in algoritmen, software en hardware ervoor gezorgd hebben dat we niet langer een supercomputer en veel geduld nodig hebben om de berekeningen en bewijzen uit te voeren, maar met een laptop al een heel eind komen.

Hierdoor kunnen nieuwe ideeën veel makkelijker uitgeprobeerd worden, en gaan de ontwikkelingen in het vakgebied overeen- komstig sneller. Aan het eind van dit artikel geven we enkele voorbeelden van problemen waar tegenwoordig de aandacht naar uitgaat (zie de laatste paragraaf).

Interval-aritmetiek

Een van de aspecten van een computer-as- sisted proof is de uitdaging dat computers niet zomaar kunnen werken met de reële getallen, maar zich moeten beperken tot een eindig aantal (binaire representeer- bare) getallen. Deze eindige verzameling getallen, die we noteren als F1R, wor- den de floating-point getallen of floats genoemd.

Het grootste nadeel van deze getallen is dat er afrondfouten optreden wanneer er bewerkingen op worden uitgevoerd (zoals naal zijn. Het is bijzonder lastig om met

een computer gelijkheden zoals in (1) te controleren. In plaats daarvan zullen we daarom gebruik maken van een dekpuntstelling. Dit doen we door het probleem

( )

F x = om te schrijven naar een dek-0 puntprobleem door middel van een variant van de methode van Newton. In de methode die we hier beschouwen, vervult de Ba- nach contractiestelling de rol van dekpuntstelling, maar er bestaan ook aanpakken die werken met de dekpuntstellingen van Brouwer en Schauder, zie bijvoorbeeld [5].

We zullen vervolgens laten zien dat op een goed gekozen omgeving van de numeriek gesimuleerde oplossing (benadering) aan de voorwaarden van de dekpuntstelling van Banach is voldaan. Dan garandeert deze stelling namelijk het bestaan en de uniciteit van een (echte) oplossing.

Via pen-en-papier-wiskunde zijn de voorwaarden van de Banach-contractiestelling te formuleren als een eindig aantal on- gelijkheden, waarin de numerieke simu- latie en analytische eigenschappen van het probleem verenigd worden. Om te laten zien dat er aan deze ongelijkheden voldaan wordt, zal er gewerkt moeten worden met de data die eerder door de computer geproduceerd zijn. Omdat gesimuleerde resultaten honderden of zelfs duizenden datapunten bevatten, is het praktisch om deze ongelijkheden met een computer te verifiëren. Het is cruciaal dat ongelijkheden (die contractie garande- ren) met een computer veel makkelijker te controleren zijn dan gelijkheden zoals (1). Niettemin betekent het dat de com- puter op een rigoureuze manier (zonder afrondfouten, of althans met rigoureuze controle van de afrondfouten) met de data moet omgaan.

Al sinds de jaren tachtig worden dit soort computer-assisted proofs toegepast binnen het gebied van dynamische systemen, zowel voor discrete problemen (zoals iteraties van afbeeldingen) als voor continue problemen (zoals gewone en partiële differentiaalvergelijkingen). Een bekend voorbeeld is het eerste bewijs van het Feigenbaum-vermoeden over de univer- saliteit van de Feigenbaum-constante [8].

Ook binnen de continue dynamische systemen werd in de jaren negentig al gebruik gemaakt van computer-assisted proofs, hetgeen uitmondde in het bewijs van het bestaan van de chaotische attractor in het beroemde Lorenz-systeem [12, 13].

(3)

, , . F x x x

f x x f x x

f x x

1 2 3

1 2

2 3

3 1

def=

- - -

m

m m m

J

L KKKK

^ KKK ^

^

^ N

P OOOO

h h OOO

hh

(3)

In dit geval is het niet meer mogelijk om de tussenwaardestelling te gebruiken. In plaats daarvan combineren we het gebruik van interval-aritmetiek met de methode van Newton–Kantorovich (zie Stelling 1 ver- derop). Het grote voordeel van deze aanpak is dat deze makkelijk te generaliseren is naar willekeurige, en in het bijzonder oneindige, dimensies.

De methode van Newton–Kantorovich We concentreren ons eerst op het eindig- dimensionale geval: we willen het bestaan van een nulpunt van een afbeelding

:

F Rⁿ"Rⁿ aantonen met behulp van de computer. De eerste stap bestaat uit het vinden van een goede numerieke benade- ring voor een nulpunt van F (bijvoorbeeld met de klassieke methode van Newton).

De tweede stap bestaat uit een combinatie van analyse op papier en interval-aritmetiek. Het idee is om te bewijzen dat er een exact nulpunt van F bestaat in de buurt van de numerieke benadering door de con- dities van de stelling van Newton–Kantoro- vich te controleren.

In het Newton-algoritme wordt een nulpunt benaderd met het iteratieve proces

( )

xk=N xk 1- , k!N. Hierbij is :N Rⁿ"Rⁿ de klassieke Newton-afbeelding gegeven door

( ) ( ),

N x ^def= -x 6DF x^ h@^-¹F x (4) en x₀!Rⁿ is een benadering (een gok) voor het nulpunt van F.

specifieke >m mc aantonen dat de dynamica van de logistische afbeelding chaotisch is door te bewijzen dat er zo’n periodieke baan van periode drie bestaat. Het is in het algemeen lastig (of zelfs onmogelijk) om een periodieke baan van een niet-lineair probleem met een pen-en-papier- analyse te vinden, en dit is typisch een si- tuatie waarin numerieke simulaties handig zijn en computer-assisted bewijzen goed tot hun recht komen.

Een voor de hand liggende strategie om een periodieke baan van periode drie te vinden is het oplossen van de vergelijking

( )

g xm =0, waarbij

( ) ,

g xm ^def=f f f xm^ ^ ^ hhhm m -x

zie Figuur 1(b). Om een eenvoudig voor- beeld van een computer-assisted bewijs te geven merken we op dat het probleem

( )

g xm =0 kan worden opgelost met behulp van interval-aritmetiek én de tussenwaardestelling: door aan te tonen dat g xm^ h₁ !I₁ en g xm^ h₂ !I₂ voor zekere x₁<x₂, waarbij I₁1[ , ]03 en I₂1[-3, ]0 (of andersom), garandeert de tussenwaardestelling het bestaan van een nulpunt van gm in ,6x x₁ ₂@. In dit speciale geval zijn die berekeningen zelfs met enig geduld nog wel met de hand te doen. Een nadeel van deze specifieke aanpak is dat de methode gebruik maakt van puur een- dimensionale argumenten die lastig te generaliseren zijn naar hoger dimensionale problemen.

Een alternatief voor de bovenstaande methode is om een nulpunt te vinden van de afbeelding F 0 1_m:[ , ]³"[ , ]0 1³ gegeven door

kan de bovengrens gevonden worden door naar boven af te ronden en het grootste resultaat te kiezen.

Op soortgelijke wijze kunnen we de andere elementaire operaties , , :5 6 { I_F

r

#I_F

r

"I_F

r

definiëren. Het op deze wijze rekenen met intervallen noemen we interval aritmetiek. Minder elementaire func- ties zoals exp, log, sin, et cetera, kunnen met behulp van bijvoorbeeld de stelling van Taylor zo geïmplementeerd worden dat ze aan de cruciale inclusie-eigenschap

voor alle

( ): ( )

f x x I f I I

I I

F F

d 1 d

!

r r

" ,

(2) voldoen. Het is goed om te beseffen dat operaties in interval aritmetiek niet uniek gedefinieerd zijn, maar dat dit geen probleem is. Voor verschillende implementa- ties van interval aritmetiek kan ( )f I een verschillend resultaat opleveren, maar de inclusie (2) moet altijd gegarandeerd zijn.

Op die manier kunnen we interval aritmetiek dus gebruiken om ongelijkheden te bewijzen.

Voor het gebruik van interval aritmetiek zijn verscheidene software-pakketten in omloop. Veel gebruikt zijn bijvoorbeeld de Intlab library voor Matlab en de Profil/

BIAS library voor C++. Daarnaast bieden ook symbolische software pakketten als Maple en Mathematica ondersteuning voor interval aritmetiek.

Voorbeeld: de logistische afbeelding Als voorbeeld bekijken we de in Figuur 1(a) gerepresenteerde logistische afbeelding f_m:[ , ]0 1 "[ , ]0 1, gegeven door

( ) ( )

f xm ^def=mx 1-x met m![ , ]0 4 een parameter. De dynamica x_k+₁=f xm( )_k gegenereerd door fm is relatief eenvoudig te beschrijven voor m!60,m_c@, waarbij m_c.3 57, . Namelijk, alle banen convergeren naar een evenwichtstoestand of een periodieke baan. Voor >m m_c is de dynamica niet meer eenvoudig te beschrijven: het systeem vertoont chaotisch gedrag. Maar hoe bewijs je dit?

Zoals besproken in de introductie be- staan er stellingen, forcing theorems, die het bestaan van speciale banen gebruiken om uitspraken te doen over het globaal gedrag van een dynamisch systeem. Een beroemd voorbeeld van een forcing theorem voor discrete dynamische systemen op de reële lijn is dat het bestaan van een periodieke baan van periode drie chaos impliceert [11].

Met andere woorden, we kunnen voor een

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

-1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1

x gλ(x)

(a) (b)

Figuur 1 (a) Het beroemde bifurcatiediagram van Feigenbaum. Het asymptotisch gedrag van de logistische afbeelding is in kaart gebracht (te zien op de verticale as) voor m![ , ]0 4 met behulp van de computer. De numerieke simulatie doet vermoeden dat de dynamica chaotisch is voor m>m_c.3 57, . (b) De grafiek en nulpunten van gm voor m=3 95, . De De open cirkels zijn de twee evenwichtspunten van fm. De kruisjes en zwarte punten corresponderen met twee verschillende (niet triviale) periodieke banen van periode drie.

(4)

den bepaald middels een afschatting van het residu ( )eT x

t

-x

t

e. De contractie sterkte kan worden bepaald door ;DT x( ); te begrenzen op B x_r( )

t

. Omdat we a priori niet weten voor welke straal r de afbeelding T een contractie is, houden we deze va- riabel in de analyse op papier. We maken vervolgens gebruik van de computer om r te bepalen.

De straal r moet namelijk enerzijds groot genoeg zijn om ervoor te zorgen dat het exacte nulpunt bevat is in B x_r^ h

t

. An- derzijds moet de straal klein genoeg zijn om de contractiviteit van T te garanderen (zie Figuur 2). Door op zoek te gaan naar een zo klein mogelijke straal waarvoor T een contractie is kunnen we op een wiskundig exacte manier uitspraken doen over de kwaliteit van de numerieke simulatie. Aan de andere kant, door op zoek te gaan naar een zo groot mogelijke straal kunnen we uitspraken doen over de grootte van het gebied waarin de exacte oplossing uniek is.

De precieze formulering van de hierbo- ven beschreven methode is samengevat in de volgende stelling en staat bekend als de geparametriseerde Newton–Kantorovich- methode (het bewijs is niet ingewikkeld, zie bijvoorbeeld [7]).

Stelling 1. Laat T een differentieerbare afbeelding van een Banachruimte X naar zichzelf zijn, en laat x

t

!X. Neem aan dat er een getal Y> en een functie 0

:

Z R⁺"R⁺ bestaan zodat

( ) ,

( ) ( ).

sup

T x x Y

DT x rv Z r

( ) v B 01

e e

; ;

#

# -

+

!

t t

t

Als er een straal r

t

>0 bestaat zodat

( ) ,

Y Z r r+

t t t

<r (5) dan is :T B x_r

t

( )

t

"B x_r

t

( )

t

een contractie.

Aangezien dekpunten van T correspon- deren met nulpunten van F kunnen we dus aantonen dat F een uniek nulpunt in de buurt van de numerieke benadering x

t

heeft door de ongelijkheid te controleren.

Het werk is nu ‘verplaatst’ naar het vin- den van goede schattingen Y en ( )Z r ; in het bijzonder moet minstens gelden dat Z strikt kleiner dan 1 is.

We illustreren dit aan de hand van het voorbeeld uit de vorige paragraaf.

We kiezen m=3 95, en proberen een nulpunt te vinden van Fm, gegeven in (3). We gebruiken een benadering voor niet-lineaire afbeeldingen F door-

gaans moeilijk, omdat de analyse van de inverse [DF x( )]^-¹ in (4) lastig is vanwege de afhankelijkheid in x.

Omdat de contractiesterkte van de Newton-afbeelding blijft toenemen naarmate we dichter in de buurt van het nulpunt komen, hebben we de vrijheid om kleine wijzigingen in N aan te brengen zonder de contractie-eigenschap te verlie- zen. In het bijzonder kunnen we de inverse [DF x( )]^-¹ vervangen door een numerieke benadering A.[DF x( )]

t

^-¹ en daarmee de afhankelijkheid in x vermijden (hier evalu- eren we DF dus in een vast gekozen punt x

t

). Dit levert een nieuwe Newton-achtige afbeelding :T Rⁿ"Rⁿ op, gegeven door

( ) ( ).

T x ^def= -x AF x

Onder de aanname dat de matrix A injectief is zullen de vaste punten van T precies de nulpunten van F zijn. De contractiesterkte van T is weliswaar zwakker dan die van N, maar het gebruik van deze afbeelding maakt het mogelijk om de analyse daadwerkelijk uit te voeren.

Zoals al eerder aangegeven is het doel nu om een omgeving U=B x_r^ h

t

te vinden (een gesloten bol om x

t

met straal r> ) 0 waarop T een contractie is. Hier hebben we twee ingrediënten voor nodig: de kwa- liteit van de numerieke benadering x

t

^en

een Lipschitz-constante voor T. De kwali- teit van de numerieke benadering kan wor- Een belangrijke aanname in de methode

van Newton is dat de afgeleide DF van F in het exacte nulpunt bestaat, en continu en inverteerbaar is. Dit betekent in het bijzonder dat DF x inverteerbaar is in een ( ) kleine omgeving van het exacte nulpunt.

In deze omgeving geldt dat ( )F x = dan 0 en slechts dan als ( )N x = . Met andere x woorden, nulpunten van F zijn vaste pun- ten van N en andersom.

De methode van Newton is in praktijk bijzonder effectief. Dit komt doordat de Newton-afbeelding een contractie is in een omgeving van het exacte nulpunt. Boven- dien is de contractiesterkte van N ‘onein- dig’ goed in de buurt van het exacte nulpunt. Om preciezer te zijn, een afbeelding

:

T X"X op een genormeerde vectorruim- te X is een contractie op U1X als er een constante L!( , )0 1 bestaat zodat

( ) ( ) , ,

T x -T y #L x y- voor allex y!U waarbij $e e de norm op X is. De constante L wordt doorgaans een Lipschitz-constante voor T op U genoemd. In het geval dat T differentieerbaar is (en U convex) kan de middelwaardestelling worden gebruikt om een Lipschitz-constante te bepalen: kies L zodat ;DT x( ); #L voor alle x!U, waarbij

( ) DT x

; ; de operator-norm van de lineaire operator DT x is. Hoe kleiner L, hoe ster-( ) ker de contractie.

De dekpuntstelling van Banach stelt dat een contractie op een volledige ruimte een uniek dekpunt x⁾ heeft. Met andere woor- den, als U volledig is en T een contractie zodat ( )T U 1U, dan bestaat er een uniek element x⁾!U zodat T x^ h⁾ =x⁾. In het bijzonder geldt lim_k" T x^k =x)

3 ^ h voor elke startwaarde x!U, zie Figuur 2.

Om nu terug te komen op de methode van Newton, die werkt zo goed omdat DN x( )⁾ =0. Dit heeft als gevolg dat de Lipschitz-constante voor N willekeurig klein wordt als we maar dicht genoeg in de buurt van x⁾ zijn.

Stel nu dat we een goede numerie- ke benadering x

t

!Rⁿ voor een nulpunt van F hebben gevonden. Het idee is om deze numerieke benadering te gebruiken om te bewijzen dat er een uniek nul- punt van F bestaat in een omgeving van x

t

middels de Banach-contractiestelling.

In het ideale geval zouden we een expli- ciete Lipschitz-constante voor N uitreke- nen om aan te tonen dat de afbeelding een contractie is in een omgeving van x

t

^.

Echter, een dergelijke analyse van N is

Figuur 2 Een schematische weergave van de dek- punt-stelling van Banach. We starten met een omgeving U=B x_r^{^ h}

t

van een numerieke benadering x

t

. Omdat T een contractie is die U in zichzelf afbeeld is U₁^def=T U^ h1U een verzameling met een kleinere diameter. De verzameling U₁ wordt vervolgens door T afgebeeld op een nog ‘kleinere’

verzameling U₂^def=T U_{^ h}₁ 1U₁, enzovoort. Uit het uniform kleiner worden van het rijtje verzamelingen U_k (en de vol- ledigheid van U) volgt het bestaan van een uniek punt x⁾!U (het dekpunt van T) zodat lim_k" T x^k =x)

3 ^ h

voor elke x!U.

(5)

:

( ) ( , ) ( , , )

( , , )

. F

a a a a

a

a x a b a a b a a a b a a a

ka b a a

2 3

k k k k

0 1 2 3

0 0

1 0 0

2 1 0 1

3 2 0 1 2

1 0 1

7 h h

h f h - - - -

- - -

J

L KKKK KKKK KKKK KKKK KKK

J

L KKKK KKKK KKKK KKKK KKK N

P OOOO OOOO OOOO OOOO OOO

N

P OOOO OOOO OOOO OOOO OOO

(10)

Er zijn in de praktijk (door het kiezen van andere bases dan de machten t^k) nog vele andere keuzes voor dit soort rijtjesruimtes mogelijk. Zo wordt er ook vaak gebruikt gemaakt van de Fourier-, Legendre- of Che- byshev-basis. De precieze keuzes die worden gemaakt voor de afbeelding F X: "X’ en de ruimten X en ’X verschillen sterk per probleem en hangen vaak nauw samen met hoe de bijna-inverse A geconstrueerd wordt.

Al deze methoden hebben één ding gemeen: F is nu een afbeelding tussen oneindig-dimensionale ruimten. Op het eerste gezicht zou dat weinig moeten ver- anderen, want Stelling 1 is ook waar voor oneindig-dimensionale ruimten. Echter, omdat we nu te maken hebben met oneindig-dimensionale afbeeldingen is het door- gaans niet meer mogelijk om de functie F expliciet uit te rekenen. In het geval dat F een afbeelding is tussen rijtjesruimten is dit meteen duidelijk: F heeft immers on- eindig veel termen. Ook het kiezen van de benadering A voor de inverse DF x( )

t

^-¹^is

minder makkelijk.

Niettemin, in veel gevallen kunnen we F schrijven als F=Fⁿ+F³, waarbij Fⁿ al- leen afhangt van een n-dimensionale deel- ruimte van X en afbeeldt op een n-dimen- sionale deelruimte van ’X . Hierbij willen we Fⁿ zodanig kiezen dat het al het ‘interes- sante’ gedrag van F bevat. Als verder de staart F³ voldoende ‘klein’ is, dan kunnen we dus een bijna-nulpunt x

t

van F vinden simpelweg door een bijna-nulpunt van de eindig-dimensionale afbeelding Fⁿ te berekenen. Voor de verdere analyse van het probleem willen we de oneindig-dimensi- onale component F³ zo kiezen dat deze analytisch ‘makkelijk’ te bestuderen is. Het maken van een juiste keuze is doorgaans moeilijk wanneer F niet-lineaire termen be- vat. Er moet in alle niet-lineaire gevallen (dus eigenlijk in alle interessante gevallen) dus goed nagedacht worden over de ge- volgen van de keuzes die gemaakt worden voor de functie F, de ruimtes X en ’X , en de projecties op de eindige dimensionale deelruimtes (waarin de numerieke simulatie wordt uitgevoerd).

, r , r , .

14 05 ²-0 94 +5 7 10$ ^-³<0 (8) Met interval-aritmetiek is dan eenvoudig aan te tonen dat voor elke r in het in- terval [ ,0 007 0 06 aan (8) wordt voldaan, , , ] en daarmee is het bestaan van een unieke periodieke baan van periode drie in B xr( )

t

bewezen. Dus de logistische afbeelding is chaotisch voor m=3 95, op grond van de forcing theorem.

Oneindige dimensies

Zoals we gezien hebben in bovenstaand voorbeeld is het mogelijk om met behulp van Stelling 1 en interval-aritmetiek het bestaan van nulpunten van :F Rⁿ"Rⁿ aan te tonen. Om deze methode ook op differentiaalvergelijkingen toe te kunnen passen moeten we deze uitbreiden naar een oneindig-dimensionale context.

Als voorbeeld bekijken we het beginwaardeprobleem

’( ) ( ( )) met [ , ],

( ) .

x t f x t t

x x

0 1

0 ₀

= ! ) =

Allereerst zullen we dit herschrijven als een probleem van de vorm ( )F x = op een 0 oneindig-dimensionale ruimte. Een natuurlijke keuze is om de differentiaalvergelijking te integreren en dan de afbeelding

: ([ , ]) ([ , ])

F C 0 1 "C 0 1 te definiëren door

( )( ) ( ) ( ( )) .

F x t x t x f x s ds

t 0

0

def= - -

#

⁽⁹⁾

Immers, ( )( )F x 0 =x( )0 - en x0 F x( ) ( )’ t =

’( ) ( ( ))

x t -f x t , en dus geldt dat ( )F x = 0 dan en slechts dan als ’( )x t =f t( ) én

( ) x 0 =x0.

Een alternatieve strategie is om x als een machtreeks ( )x t =

/

_k³₌₀a t_k ^k te schrijven. Dit is een tweede natuurlijke keuze in het geval f reëel analytisch is. In dat geval kunnen we ( ( ))f x t represente- ren met Taylor-coëfficiënten ( , , )b ak 0fak, zodat

( ( )) ( , , ) .

f x t f a tk k b a a t

k k k k

0 k 0

0

= = f

3 3

= =

f

/

p

/

Omdat ’( )x t =

/

³_k₌₀(k+1)a_k₊₁t^k^{, is het} beginwaardeprobleem equivalent met het oneindige stelsel van algebraïsche vergelijkingen

(k+1)ak+1=b ak( ,0f, ),ak k=0 1 2, , ,f, aangevuld met a₀=x₀. In dit geval kunnen we F dus definiëren door

, , , x

0 19 0 60 0 95

t

= J

L KKKK KKK

N

P OOOO OOO

van een periodieke baan van periode drie en een bijna-inverse

, , ,

. A

0 47 0 17 0 13

0 62 1 54 0 17

0 17 0 43 0 33

def= -

- -

- J

L KKKK KKK

N

P OOOO OOO

De volgende stap is om de begrenzingen Y en ( )Z r uit te rekenen.

Om het residu te begrenzen is het voldoende om in te zien dat

( ) ( ).

T x

t

- = -x

t

AF xm

t

De bovenstaande term kan worden be- grensd met behulp van interval aritmetiek.

Dit levert Y 5 7 10^def= , $ ^-³ als maat voor het residu (hoe kleiner hoe beter) van de nu- merieke benadering x

t

. Om de afschatting

( )

Z r te bepalen, laat v!R³ met v #1, waarbij $e e de Euclidische norm is. Dan geeft de driehoeksongelijkheid

( ) ( )

( ( ))

[ ( ) ( )] .

DT x rv I ADF x rv I ADF x

A DF x rv DF x

#

+ = - +

-

+ + -

m m

t t

t

t t (6)

De eerste term in (6) is opnieuw eenvoudig te begrenzen met behulp van interval-aritmetiek, en we vinden ,0 06 als maat voor de kwaliteit (weer: hoe kleiner hoe beter) van de bijna-inverse A. De tweede term in (6) kunnen we niet direct afschatten met behulp van de computer, want die heeft een functionele afhankelijkheid van r. Er moet dus eerst enige analyse op papier worden uitgevoerd:

( ) ( ) .

DF x rv DF x r

v v

v

2 0

0 0 0

0 0

1 2

3

+ - = - m

m

t

m

t

J

L KKKK KKK

N

P OOOO OOO (7) Omdat we hebben aangenomen dat ve e #1, geldt

. v

v v 0 0

0 0

0

0 1

1 2

3

# J

L KKKK KKK

N

P OOOO OOO

We kunnen dus de tweede term in (6) afschatten met

( ( ) ( )) ,

A DF x rvm t+ -DF xm t # m2 r A en met behulp van interval aritmetiek een begrenzing voor A (en dus de bo- venstaande term) uitrekenen. We vinden dat we 2m A kunnen afschatten met

,

14 05. Door de bovenstaande berekeningen te combineren concluderen we dat

( ) , ,

Z r^def=14 05r+0 06, en de ongelijkheid wordt

(6)

voegen daarom een extra conditie aan (12) toe, een zogenaamde fase-conditie, om de tijd-parametrisatie van de periodieke baan uniek vast te leggen. We eisen daartoe dat de baan op t= in een vooraf gespecifi-0 ceerd vlak in de fase-ruimte R³ ligt.

Aangezien het gaat om een periodieke oplossing ligt het voor de hand om een startformulering ( )F x = , vergelijk (1), te 0 kiezen in termen van Fourier-coëfficiënten:

( ) ,

x t a e

y t b e

z t c e

k ik t k

Z

=

!

~

/ / /

waarbij ~=2r/L. Het invullen van de bovenstaande expansies in (12) en de fase-conditie levert een oneindig stelsel van algebraïsche vergelijkingen op voor de on- bekende rijtjes coëfficiënten , ,a b c en de frequentie ~. We schrijven dat stelsel hier niet expliciet op, maar het is in hoge mate analoog aan (10). We passen nu de in de vorige paragraaf ontwikkelde methode toe op dit probleem en kunnen daarmee bijvoorbeeld bewijzen dat de kromme (het resultaat van een simulatie) in Figuur 3 gegarandeerd een periodieke baan van (12) representeert.

Een radiaal symmetrische PDV

We kunnen de methode ook toepassen om oplossingen van partiële differentiaalvergelijkingen (PDVs) te vinden. Als voorbeeld kijken we naar gelokaliseerde oplossingen van de stationaire niet-lineaire Schrödinger vergelijking

( )x ( )x ( )x ³ 0,

} } }

D - + = (13)

waarbij }: R³"R. Wanneer we voor deze PDV op zoek gaan naar een radiaal symmetrische oplossing schrijven we

( )x f x(| |)

} = voor een : [ , )f 0 3 "R met (waarvan A dan weer een bijna inverse

is). Het vinden van een scherpe afschatting van het verschil tussen DF x( )

t

en zijn benadering is het lastigste gedeelte van de hele analyse. Bovendien is het helaas niet voor elk probleem zo eenvoudig om het probleem op een diagonaal dominante manier te formuleren (zoals in ons voorbeeld). En zelfs als dat wel kan, is het nog maar de vraag of er een eindig-dimensionale projectie van het probleem gevonden kan worden die groot genoeg is om het essentiële gedrag van het oneindig-dimensionale niet-lineaire probleem te ‘vangen’

(zodat we schattingen kunnen vinden die aan de voorwaarden van Stelling 1 voldoen), maar die anderzijds klein genoeg is om mee te kunnen rekenen op een computer. Tot slot van dit artikel beschrijven we nog enkele toepassingen waarin het gelukt is om deze uitdagingen het hoofd te bieden.

Toepassingen

In deze laatste sectie geven we twee (om- wille van de leesbaarheid relatief eenvoudige) voorbeelden van toepassingen. Er is een heel scala aan problemen uit de dynamische systemen waar momenteel aan gewerkt wordt: evolutionaire partiële dif- ferentiaalvergelijkingen [1], delay vergelij- kingen [9], sterk indefiniete problemen [6], verbindende banen [2, 10], domeindecom- positie [4], driedimensionale periodieke patronen, en natuurlijk het ontwikkelen van algemeen bruikbare software.

Periodieke banen in het Lorenz-systeem In dit voorbeeld zoeken we periodieke banen in het klassieke Lorenz systeem:

( )

( ) , [ , ],

( ) ( ), ( ) ( ), ( ) ( ).

x y z

y x

x z y

xy z t L

x x L y y L z z L

10

28 0

0 0 0

38 !

=

-

- -

-

= = =

o o o

J

L KKKK KKK

J

L KKKK KKKK N

P OOOO OOO

N

P OOOO OOOO Z

[

\ ]]]]]]

]]]]

(12)

Zowel de periodieke baan t7( ( ), ( ), ( ))x t y t z t als de bijbehorende periode L> zijn on-0 bekend (dus het bepalen van de periode is onderdeel van het probleem).

Als t7( ( ), ( ), ( ))x t y t z t een oplossing is van (12), dan beschrijft de ‘verschoven’ afbeelding t7( (x t+x), (y t+x), (z t+x)) de- zelfde periodieke baan (voor elke x!R), en is dus ook een oplossing. Door deze translatie-symmetrie is de oplossing niet lokaal uniek, dus een contractiestelling is in deze formulering niet toepasbaar. We De lineaire operator A X: ’"X, die

dienst doet als bijna-inverse van DF x( )

t

^,

heeft in het oneindig-dimensionale geval ook een aparte behandeling nodig. Ook die moet worden gesplitst in een eindig-di- mensionaal deel Aⁿ (een matrix) en een on- eindig-dimensionaal ‘staartstuk’ A³. Aan- gezien Fⁿ geacht wordt al het interessante gedrag van F te bevatten, kunnen we voor de matrix Aⁿ een numerieke benadering voor de inverse van de Jacobiaan DF xⁿ( )

t

kiezen. Het oneindig-dimensionale deel A³ is (hopelijk) klein, maar mag niet een soort nulafbeelding zijn, want A moet wel injec- tief gekozen worden om correspondentie tussen dekpunten van T en nulpunten van F te garanderen.

In het voorbeeld van de afbeelding ge- definieerd in (10) is de term ka_k dominant voor grote k. Daarom kunnen we proberen om DF x( )

t

te benaderen door

( )

. DF x

DF x n

n 0

0 1

n

j

. +

t J t

L KKKK KKKK KK

N

P OOOO OOOO OO (11) In dit geval kunnen we A als volgt con- strueren. Eerst berekenen we (numeriek) een bijna-inverse voor DF xⁿ( )

t

, namelijk

( )

Aⁿ.DF xⁿ

t

^-¹. Vervolgens kiezen we

. A

A n 0 n

1 0

11

n

def

j

=

+ J

L KKKK KKKK KKKKK

N

P OOOO OOOO OOOOO

Er geldt dan dat A DF x( )

t

een goede benadering van de identiteit is (hoe goed, dat moet expliciet worden afgeschat).

Grote delen van de eindig-dimensionale analyse uit de vorige paragraaf kunnen nu worden gegeneraliseerd naar het oneindig-dimensionale geval. Elk van de af te schatten termen omvat nu wel twee aparte analyses, namelijk die van het n-di- mensionale gedeelte (dat net als in het eindig-dimensionale geval interval aritme- tisch berekend kan worden) en die van het staartstuk. Dat laatste gedeelte is oneindig-dimensionaal (dus ongeschikt voor een directe computerberekening) en moet met behulp van een combinatie van functi- onaal-analytische ingrediënten en interval- aritmetiek afgeschat worden.

De kern van het probleem zit hem echter toch vooral in de benadering van de (oneindig-dimensionale) afgeleide DF x( )

t

door een eenvoudigere lineaire operator

15 20 5 10 0

x -5 -10 -20 -15 -40 -20 0 20 5 10 15 20 25 30 35 40 45

40

y

z

Figuur 3 Een gevalideerde periodieke baan van periode ,

L.25 03 in het Lorenz systeem.

(7)

gevonden simulatie w

t

, waarbij voor een uniform rooster is gekozen met m=1000 punten, is afgebeeld in Figuur 5, evenals de bijbehorende functie f

t

. Het bewijs voor het bestaan van een (echte) oplossing van ( )F w = vlakbij w0

t

, en daarmee van een radiaal-symmetrische gelokaliseerde oplossing van de PDV (13), volgt nu, na enige zorgvuldige analyse en met behulp van een computationeel gedeelte [3], uit de methode zoals beschreven in de vorige

paragraaf. s

sing f van de GDV wordt gegeven door

( ) ( ) ( )

f r =_r¹e^-^r 1-e^-^r w e^-^r, zie Figuur 5.

We merken op dat (15) op hoofdlijnen de- zelfde vorm heeft als de formulering (9).

Voor dit probleem kiezen we als ein- dig-dimensionale benadering w

t

van de op- lossing een stuksgewijs lineaire spline. Elke functie op [ , ]C 0 1 kan worden opgedeeld in zo’n spline-deel met m roosterpunten, en een restant, waarbij het restant kleiner wordt naarmate het aantal roosterpunten m toeneemt, zie Figuur 4. De numeriek de eigenschap dat lim_r" 3f r( )=0. Dit re-

duceert de PDV (13) tot

’’( ) ’( ) ( ) ( ) , f r +2r f r -f r +f r ³=0 (14) een singuliere niet-autonome gewone differentiaalvergelijking (GDV) op ( , )0 3 . Nu- meriek rekenen op een oneindig domein is lastig, maar door deze GDV te transforme- ren naar het interval ( , )0 1 en de methode van Green toe te passen, kan aangetoond worden dat (14) een begrensde oplossing heeft wanneer er een functie w!C 0 1[ , ]) bestaat die voor alle t![ , ]0 1 aan de volgende integraalvergelijking voldoet:

( ) ( )

( )

( ) ( )

( )

( )( )

( ) .

log

log F w w t

t t

s

s s w s ds

t s

s s s

w s ds

21 1 1

21

11 1 1

0

t

2 2

3 3

0

3

2

4 3

1

def= - + -

- - + -

=

#

(15) De relatie tussen de oplossing w van de integraalvergelijking (15) en de oplos-

0.0 0.2 0.4 0.6 0.8 1.0t

1 2 3 4 5 w

0.0 0.5 1.0 1.5 2.0 2.5 3.0r 1

2 3 4 f

0.0 0.2 0.4 0.6 0.8 1.0

2 4 6 8

=

0.0 0.2 0.4 0.6 0.8 1.0

2 4 6 8

+

_0.2 _0.4 _0.6 _0.8 _1.0

-1.0-0.5 0.51.0

Figuur 4 Een afbeelding op [ , ]0 1 kan worden opgedeeld in een spline deel (met 6 roosterpunten in dit voorbeeld) en een restant.

Figuur 5 Links de functie ( )w t

t

^{op [ , ]}^{0 1}. Rechts de corresponderende functie ( )f r

t

^{op ( , )}^{0 3}^.

1 G. Arioli en H. Koch, Integration of dissipa- tive partial differential equations: a case study, SIAM J. Appl. Dyn. Syst. 9(3) (2010), 1119–1133.

2 J. B. van den Berg, A. Deschenes, J. P. Les- sard en J. D. Mireles-James, Stationary Coex- istence of Hexagons and Rolls via Rigorous Computations, SIAM J. Appl. Dyn. Syst. 14(2) (2015), 942–979.

3 J. B. van den Berg, C. M. Groothedde en J. F. Williams, Rigorous computation of a radially symmetric localized solution in a Ginzburg–Landau problem, SIAM J. Appl.

Dyn. Syst. 14 (2015), no. 1, 423–447.

4 J. B. van den Berg en R. S. S. Sheombarsing, Rigorous numerics for ODEs using Cheby- shev series and domain decomposition, 2015. Preprint.

5 CAPD: Computer assisted proofs in dynamics, a package for rigorous numerics, http://

capd.ii.uj.edu.pl.

6 R. Castelli, M. Gameiro en J.-P. Lessard, Rig- orous numerics for ill-posed PDEs: periodic orbits in the Boussinesq equation, 2015.

Preprint.

7 S. Day, J.-P. Lessard en K. Mischaikow, Val- idated continuation for equilibria of PDEs, SIAM J. Numer. Anal. 45(4) (2007), 1398–

1424.

8 O. E. Lanford, A computer-assisted proof of the Feigenbaum conjectures, Bull. Amer.

Math. Soc. 6 (1982), 427–434.

9 J. P. Lessard, Recent advances about the uniqueness of the slowly oscillating periodic solutions of Wright’s equation, J. Differential Equations 248(5) (2010), 992–1016.

10 J. P. Lessard, J. D. Mireles-James en C. Rein- hardt, Computer assisted proof of trans- verse saddle-to-saddle connecting orbits for first order vector fields, J. Dynamics and Dif- ferential Equations 26(2) (2014), 267–313.

11 T. Li en J. A. Yorke, Period three implies chaos, The American Mathematical Monthly 82 (1975), no. 10, 985–992.

12 K. Mischaikow en M. Mrozek, Chaos in the Lorenz equations: a computer assisted proof, Bull. Amer. Math. Soc. 32(1) (1995), 66–72.

13 W. Tucker, The Lorenz attractor exists, C. R.

Acad. Sci. Paris 328(12) (1999), 1197–1202.

Referenties