Bachelorthesis Wiskunde en Informatica
Scherpe Triangulaties in de Eindige Elementen
Methode
Raymond van Veneti¨e
15 juli 2014
Begeleiding: Dr. Jan Brandts & Dr. Leen Torenvliet
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
Instituut voor Informatica
Korteweg-de Vries Instituut voor Wiskunde
Samenvatting
We beginnen met een korte introductie in de theorie van de eindige elementen me-thode voor reactie-diffusieproblemen. Vervolgens laten we zien hoe we deze meme-thode (handig) kunnen implementeren. We tonen aan dat de exacte oplossing van een reactie-diffusieprobleem aan het maximum principe voldoet. Met een aantal (extra) eisen aan de triangulatie van het domein kunnen we laten zien dat de eindige elementen oplossing ook aan dit maximum principe voldoet. Een zo een eis is dat alle dihedrale hoeken van de simplices uit de triangulatie scherp zijn.
We bekijken vervolgens triangulaties van het vierkant waarbij de hoekpunten, van de driehoeken, op een regelmatig rooster in het vierkant liggen. We onderzoeken of we met behulp van de computer een scherpe roostertriangulatie van het vierkant kunnen vinden. We presenteren een aanpassing aan de methode die in het eerdere werk [13] gebruikt werd om dergelijke scherpe roostertriangulaties te vinden. We zijn hiermee in staat roostertriangulaties voor fijnere roosters van het vierkant te vinden.
We onderzoeken vervolgens of we een scherpe roostertriangulatie van de kubus kunnen vinden. Dit probleem blijkt een stuk complexer dan het vinden van dergelijke triangu-laties voor het vierkant. We presenteren een algoritme om de Conforme Verzameling van de kubus te construeren. Dit is een verzameling facetten die mogelijkerwijs in een scherpe roostertriangulatie van de kubus kunnen zitten. Deze verzameling blijkt leeg te zijn voor de eerste 29 regelmatige roosters van de kubus. Dit impliceert dat er voor deze roosters geen scherpe roostertriangulatie bestaat. Voor het 30ste regelmatige rooster van de kubus vinden we dat de Conforme Verzameling niet-leeg is. We hebben daarom het vermoeden dat er voor dit rooster een scherpe roostertriangulatie van de kubus bestaat.
Titel: Scherpe Triangulaties in de Eindige Elementen Methode Auteur: Raymond van Veneti¨e, r.v.venetie@gmail.com, 10004627 Begeleiding: Dr. Jan Brandts & Dr. Leen Torenvliet
Einddatum: 15 juli 2014
Korteweg-de Vries Instituut voor Wiskunde Universiteit van Amsterdam
Science Park 904, 1098 XH Amsterdam http://www.science.uva.nl/math
Inhoudsopgave
1. Inleiding 6
2. De eindige elementen methode 8
2.1. Definities en notatie . . . 8
2.2. Het eendimensionale geval . . . 9
2.2.1. Zwakke differentieerbaarheid . . . 10
2.2.2. Zwakke formulering . . . 10
2.2.3. De Galerkin methode . . . 11
2.2.4. De eindige elementen methode . . . 13
2.2.5. De volledig discrete eindige elementen methode . . . 14
2.2.6. Bovengrens fout eindige elementen oplossing . . . 15
2.2.7. Voorbeeld . . . 17
2.3. Het reactie-diffusieprobleem . . . 18
2.3.1. De zwakke formulering van het reactie-diffusieprobleem . . . 18
2.3.2. De Galerkin methode . . . 19
2.3.3. De eindige elementen methode . . . 20
2.3.4. De volledig discrete eindige elementen methode . . . 21
2.3.5. Bovengrens fout eindige elementen oplossing . . . 22
3. Implementatie van de volledig discrete eindige elementen methode 23 3.1. Het eendimensionale geval . . . 23
3.1.1. Discretisatie van de inproducten (f, ψj) . . . 24
3.2. Het tweedimensionale reactie-diffusieprobleem . . . 25
3.2.1. Nodale basisfuncties . . . 25
3.2.2. Opdelen matrices . . . 26
3.2.3. Transformatie nodale basisfuncties . . . 27
3.2.4. Samenstellen massamatrix . . . 28
3.2.5. Samenstellen stijfheidsmatrix . . . 29
3.2.6. Volledig discrete eindige elementen oplossing . . . 30
3.2.7. Constante reactieterm . . . 30
3.2.8. De norm van een functie uit S(T ) berekenen . . . 30
3.2.9. De uniforme verfijning van een partitie . . . 31
3.3. Het driedimensionale reactie-diffusieprobleem . . . 33
3.3.1. Samenstellen matrices . . . 34
3.3.2. De uniforme verfijning van een triangulatie . . . 35
3.3.3. Voorbeeld . . . 36
4. Het maximum principe 38 4.1. Het exacte maximum principe . . . 38
4.2. Het discrete maximum principe . . . 39
4.2.1. Herformulering van het discrete maximum principe . . . 40
4.2.2. Het discrete maximum principe en de Stieltjes matrix . . . 41
4.2.3. Het diffusieprobleem . . . 43
4.2.4. Het reactie-diffusieprobleem . . . 44
5. Scherpe roostertriangulaties van het vierkant 47 5.1. Voorbereiding . . . 47
5.2. De methode van Schukking en Weikamp . . . 49
5.3. De Conforme Verzameling . . . 50
5.3.1. Definitie . . . 50
5.3.2. De Conforme Verzameling construeren . . . 52
5.3.3. De Gezellige Verzameling . . . 53
5.4. De implementatie van het genereren van de Conforme Verzameling . . . . 54
5.4.1. Bepalen of een lijnstuk een conforme zijde is . . . 54
5.4.2. Het bepalen van de halfruimtes . . . 56
5.4.3. De Conforme Verzameling genereren . . . 56
5.4.4. De Gezellige Verzameling genereren . . . 58
5.5. Optimalisatie van de algoritme . . . 59
5.5.1. Het controleren op conformiteit . . . 59
5.5.2. Symmetrie¨en van het vierkant . . . 60
5.5.3. Het fundamenteel domein . . . 61
5.5.4. Implementatie van de symmetrie¨en . . . 63
5.5.5. Resultaat na de optimalisaties . . . 63
5.6. Scherpe roostertriangulaties vinden . . . 64
5.6.1. Preciezere beschrijving . . . 64
5.6.2. Implementatie en resultaat van de algoritme . . . 65
6. Scherpe roostertriangulatie van de kubus 67 6.1. Scherpe roostertetra¨eders . . . 67
6.1.1. Definities . . . 67
6.1.2. Ori¨entatie van de normaalvectoren . . . 68
6.1.3. Normaalvectoren op de facetten bepalen . . . 69
6.1.4. Scherpe roostertetra¨eders genereren . . . 71
6.1.5. Noodzaak andere algoritme . . . 72
6.2. De Conforme Verzameling . . . 73
6.2.2. De Conforme Verzameling construeren . . . 74
6.3. Implementatie van de algoritme . . . 75
6.3.1. Het controleren op conformiteit . . . 75
6.3.2. Datastructuur voor een verzameling roosterdriehoeken . . . 75
6.3.3. Resultaten . . . 77
6.4. Optimalisaties . . . 78
6.4.1. Scherpe roosterdriehoeken . . . 78
6.4.2. Projectie op de facetten . . . 79
6.4.3. Symmetrie¨en van de kubus . . . 81
6.4.4. Inline code . . . 82
6.4.5. De index matrix . . . 83
6.4.6. De fundamentele Conforme Verzameling . . . 84
6.4.7. De fundamentele Conforme Verzameling construeren . . . 86
6.4.8. Parallellisatie van de algoritme . . . 88
6.5. Resultaten op de DAS-4 . . . 88 6.6. Code in C . . . 90 6.7. Openstaande vragen . . . 91 7. Discussie 92 8. Populaire samenvatting 94 Bibliografie 97 A. Eindige elementen methode 99 A.1. Het eendimensionale geval . . . 99
A.2. Het tweedimensionale geval . . . 100
A.3. Het driedimensionale geval . . . 102
B. Triangulaties op een rooster 105 B.1. Scherpe roostertriangulaties van het vierkant genereren . . . 105
1.
Inleiding
In dit verslag gaan we op zoek naar een mooie triangulatie van de eenheidskubus [0, 1]3
bestaande uit scherpe tetra¨eders. In het artikel [1] uit 2009 is voor het eerst aangetoond dat zo een scherpe triangulatie van de eenheidskubus bestaat. We zoeken echter een scherpe triangulatie met een extra eigenschap. Laat p een natuurlijk getal zijn en bekijk het regelmatige p × p × p-rooster1 van [0, 1]3. Voor welke waarde van p bestaat er een triangulatie van [0, 1]3 die alleen uit scherpe tetra¨eders bestaat waarvan de hoekpunten op het p × p × p-rooster liggen? Wat is het grofste rooster (kleinste p) waarvoor zo een triangulatie bestaat? Kan deze gevonden worden met behulp van de computer?
Deze vragen bestuderen we in Hoofdstuk 6. Eerst motiveren we deze vragen door te kijken naar de eindige elementen methode. Deze methode wordt gebruikt om oplossin-gen van parti¨ele differentiaalverlijkingen te benaderen. In Hoofdstuk 2 bestuderen we de theorie van de eindige elementen methode toegepast op het reactie-diffusieprobleem. Een (handige) implementatie van de methode geven we vervolgens in Hoofdstuk 3.
Het blijkt dat de exacte oplossingen van het reactie-diffusieprobleem voldoen aan het maximum principe. Dit principe bespreken we in Hoofdstuk 4. Uit een voorbeeld zal blijken dat de eindige elementen oplossingen niet direct aan het maximum principe vol-doen. We zullen zien dat de eindige elementen methode, onder bepaalde voorwaarden aan de triangulatie van het domein, wel oplossingen geeft die aan het maximum principe voldoen. Een van de voorwaarden is dat de triangulatie van het domein scherp is. Een logische vervolgvraag is nu of zo een scherpe triangulatie voor het domein wel bestaat, en zo ja, hoe die gevonden kan worden. In Hoofdstuk 5 gaan we op zoek gaan naar scherpe triangulaties van het eenheidsvierkant [0, 1]2. We bekijken voor welke waarden van p we een triangulatie van [0, 1]2 kunnen vinden in scherpe driehoeken met hoekpunten uit het regelmatige p × p-rooster.
We bestuderen hiervoor het eerdere werk [13]. In dit werk wordt een algoritme voor-gesteld die voor een gegeven p alle scherpe triangulaties op het p × p-rooster kan bepalen. We presenteren in Hoofdstuk 5 een alternatieve methode om dergelijke triangulaties te vinden. Met deze methode hebben we een scherpe triangulatie op het 20 × 20-rooster gevonden2. Ter vergelijking, in het eerdere werk [13] was p = 9 de hoogste waarde van p waarvoor er scherpe triangulaties op het p × p-rooster berekend konden worden.
1Punten uit dit rooster hebben dus co¨ordinaten van de vorm k
p voor k ∈ {0, . . . , p}. 2Illustratie op de titelpagina.
In Hoofdstuk 6 gaan we op zoek naar een scherpe triangulatie van de eenheidskubus op het p × p × p-rooster, zoals beschreven in de eerste alinea. Het blijkt dat het vinden van een dergelijke triangulatie voor de kubus een stuk lastiger is dan voor het vierkant. Dit komt doordat het aantal scherpe tetra¨eders op het p × p × p-rooster sterk toeneemt als we p vergroten. De algoritme3 die we in Hoofdstuk 5 presenteren voor het vierkant kunnen we op handige wijze veralgemeniseren naar de kubus. We geven vervolgens een aantal essenti¨ele optimalisaties om deze algoritme te verbeteren.
Dankwoord
Graag wil Jan Brandts bedanken voor het begeleiden van mijn bachelorthesis. Hij heeft het onderwerp aangedragen en veel van de idee¨en uit Hoofdstuk 5 en 6 komen van zijn hand. De gesprekken met hem hebben niet alleen inzicht gegeven in de theorie, maar ook in het reilen en zeilen van de Universiteit van Amsterdam. Verder wil ik ook mijn andere begeleider, Leen Torenvliet, bedanken. Zijn bereidheid om te helpen met het informatica gedeelte heb ik zeer gewaardeerd.
3
Het woord ‘algoritme’ is afgeleid van de naam Mohammed ibn Moesa al-Chwarizmi. Daar dit een man is, is het juist om te spreken over ‘de algoritme’.
2.
De eindige elementen methode
De eindige elementen methode wordt gebruikt om oplossingen van bepaalde parti¨ele differentiaalvergelijkingen te benaderen. In Sectie 2.3 bestuderen we de methode toege-past op het zogenaamde reactie-diffusieprobleem. De methode wordt in Sectie 2.2 eerst ingeleid door naar een eendimensionale vorm van het reactie-diffusieprobleem te kijken.
2.1. Definities en notatie
We geven een overzicht van de notaties die we gebruiken. We schrijven Ω ⊆ Rdvoor het domein van de functies.
Definitie 2.1. De verzameling functies die k keer continu differentieerbaar zijn op Ω noteren we met Ck(Ω). De verzameling continue functies op Ω noteren we met C0(Ω) of C(Ω).
Definitie 2.2. Noteer L2(Ω) voor de verzameling Lebesgue meetbare functies die kwa-dratisch integreerbaar zijn op Ω. We identificeren twee functies u en v als u(x) = v(x) Lebesgue bijna overal.
Definitie 2.3. Zij f, g ∈ L2(Ω), dan noteren we het standaardinproduct op L2(Ω) met (f, g) =
Z
Ω
f (x)g(x) dx. De standaardnorm wordt ge¨ınduceerd door dit inproduct,
kf k2 =p(f, f ).
Verder geldt er dat alle continue functies Lebesgue meetbaar zijn, dus de bovenstaande definitie geeft ook een inproduct en norm op de ruimte C(Ω).
Definitie 2.4. We noteren Pk(Ω) voor de verzameling polynomen van graad k op Ω. Definitie 2.5. We noteren ∂Ω voor de rand van Ω. De afsluiting van Ω geven we aan met Ω.
Definitie 2.6. Bekijk een functie f : X → Y . De beperking van f tot een deelverzame-ling Ω ⊆ X noteren we met f |Ω, oftewel
Definitie 2.7. De Kronecker delta noteren we met δij, dus
δij =
(
0 als i 6= j 1 als i = j
Verder gebruiken we een aantal operatoren uit de vectorcalculus.
Definitie 2.8. De gradi¨ent van een een scalaire functie f (x1, x2, . . . , xn) geven we met
∇f = ∂f ∂x1 , . . . , ∂f ∂xn > .
Definitie 2.9. De divergentie van een continu differentieerbaar vectorveld F geven we met divF = ∇ · F = n X i=1 ∂Fi ∂xi .
Definitie 2.10. De Laplaciaan van een scalaire functie is gedefinieerd door ∆ : Ck(Rn) → Ck−2(Rn) : f 7→ ∆f = div(∇f ) = n X i=1 ∂2f ∂x2 i .
2.2. Het eendimensionale geval
Laat a, b ∈ R met a < b gegeven zijn en schrijf I = [a, b]. Laat verder ook een f ∈ C(I) gegeven zijn. Bekijk het volgende randwaardeprobleem: vind u ∈ C2(I) zodanig dat
−u00 = f in I,
u = 0 op ∂I. (2.1)
Dit is een eendimensionaal diffusieprobleem, in de literatuur ook wel bekend als de Poissonvergelijking [2, p. 152]. We zullen later zien dat dit een vorm van het reactie-diffusieprobleem is, waarbij de reactieterm gelijk aan nul is. Het is gemakkelijk in te zien dat dit probleem een unieke oplossing heeft.
Stelling 2.11. Het randwaardeprobleem (2.1) heeft een unieke oplossing.
Bewijs. Een oplossing van (2.1) kan gevonden worden door het tweemaal toepassen van de hoofdstelling van de integraalrekening,
u(x) = − Z x a Z y a f (t) dt dy.
Stel dat z ∈ C2(I) een andere oplossing is van (2.1). Dan geldt in I dat −(u − z)00= −u00+ z00= −f + f = 0.
Hieruit volgt dat u − z een polynoom van graad 1 is, dus u − z ∈ P1(I). Aangezien u en z beide oplossingen zijn van (2.1) weten we ook dat op de rand, ∂I = {a, b}, moet gelden dat u − z = 0 − 0 = 0. Een lineaire functie die de waarde 0 heeft in twee punten, is constant 0. Oftewel u − z = 0 op I en daarmee volgt dat u = z.
2.2.1. Zwakke differentieerbaarheid
Het diffusieprobleem (2.1) kunnen we in een equivalente formulering geven, de zwakke formulering. Deze zwakke formulering kunnen we vervolgens gebruiken om de eindige elementen oplossing te bepalen. Om deze formulering te kunnen geven hebben we een aantal definities nodig. Allereerst bekijken we een functieruimte die rijker is dan C1(I),
de Solobev ruimte.
Definitie 2.12. De Sobolev ruimte noteren we met H1(I). Deze ruimte bestaat uit de functies u ∈ L2(I), die Lebesgue bijna overal differentieerbaar zijn en waarvoor geldt u0 ∈ L2(I). Verder noteren we H1
0(I) = {v ∈ H1(I) : v = 0 op ∂I}.
Schrijf nu C01(I) = {v ∈ C1(I) : v = 0 op ∂I}, dan geldt er C01(I) ⊆ H01(I). We kunnen ook een inproduct geven op deze ruimte.
Definitie 2.13. Definieer [v, w] := (v0, w0) voor alle v, w ∈ H01(I). Lemma 2.14. [v, w] definieert een inproduct op H01(I).
Bewijs. Omdat v, w ∈ H01(I) geldt dat v0, w0 ∈ L2(I) en dus is de functie [·, ·] goed
gedefinieerd.
Uit het feit dat differenti¨eren een lineaire operator is volgt direct de bilineariteit en symmetrie van [v, w]. Verder geldt er
[v, v] = (v0, v0) ≥ 0, omdat v0 ∈ L2(I) en (·, ·) een inproduct is op L2(I).
Tot slot hebben we,
[v, v] = 0 ⇐⇒ (v0, v0) = 0 ⇐⇒ v0= 0 ⇐⇒ v = 0.
De laatste equivalentie komt uit het feit dat v0 = 0 impliceert dat v constant is. Combi-neren we dit met het feit dat v nul op de rand is, dan volgt dat v nul is.
2.2.2. Zwakke formulering
Met de definities uit de vorige sectie zijn we in staat de zwakke formulering van het randwaardeprobleem (2.1) te geven.
Definitie 2.15. De zwakke formulering van het randwaardeprobleem (2.1) is: Vind u ∈ H01(I) zodat [u, v] = (f, v) voor alle v ∈ H01(I).
We laten nu zien dat een oplossing van het randwaardeprobleem (2.1) ook een oplossing van de zwakke formulering is.
Stelling 2.16. Zij u de oplossing van (2.1) en laat v ∈ H1
0(I) zijn, dan geldt dat
Bewijs. Aangezien u de oplossing is van (2.1) geldt −u00(x) = f (x) voor a < x < b. We vinden dus (−u00, v) = Z b a −u00(x)v(x) dx = Z b a f (x)v(x) dx = (f, v). Anderzijds geldt vanwege parti¨ele integratie
(−u00, v) = Z b a −u00(x)v(x) dx =−u0(x)v(x)ba− Z b a −u0(x)v0(x) dx (v ∈ H01(I) , dus nul op de rand)
= Z b
a
u0(x)v0(x) dx = [u, v].
Zolang f ∈ C(I) geldt het omgekeerde van bovenstaande stelling ook. Een oplossing van de zwakke formulering geeft ons dus direct een oplossing van het randwaardepro-bleem (2.1).
Stelling 2.17. Zij u ∈ C2(I) de oplossing van het randwaardeprobleem (2.1). Laat w ∈ H01(I) een oplossing van de zwakke formulering 2.15 zijn, dan geldt er dat u = w. Bewijs. Uit de vorige stelling weten we dat u ook een oplossing is van de zwakke formu-lering. We hebben dus
[u − w, v] = [u, v] − [w, v] = (f, v) − (f, v) = 0 voor alle v ∈ H01(I).
Neem nu v = u − w, dan geldt er dat v ∈ H01(I). Verder volgt er nu dus [u − w, v] = [u − w, u − w] = 0. Omdat [·, ·] een inproduct is impliceert dit dat u = w.
Opmerking. In de zwakke formulering zoeken we naar een functie uit H01(I), terwijl we in het diffusieprobleem (2.1) een functie uit C2(I) zoeken. Het is a priori dus niet duidelijk dat een oplossing van de zwakke formulering ook in C2(I) zit. Dat dit wel het geval is blijkt uit de bovenstaande stelling.
2.2.3. De Galerkin methode
Het vinden van een oplossing van de zwakke formulering in Definitie 2.15 geeft ons direct een oplossing van het diffusieprobleem (2.1). De ruimte H01(I) is echter nog steeds oneindig dimensionaal is, dus dit is nog geen praktische manier om een oplossing te vinden.
Het idee is nu om een eindige lineaire deelruimte V ⊆ H01(I) te kiezen. Als we H01(I) in de zwakke formulering vervangen door V krijgen we namelijk een discreet probleem: vind uV ∈ V zodat
Het vinden van een oplossing voor dit probleem heet de Galerkin methode [3, p. 53]. In Sectie 2.2.6 zullen we zien dat we een bovengrens kunnen vinden voor de fout tussen de exacte oplossing u en de Galerkin oplossing uV.
Laat nu φ1, . . . , φn ∈ H01(I) lineair onafhankelijke functies zijn die een basis vormen
voor V , dus V = span{φ1, . . . , φn}.
Lemma 2.18. Laat u de (exacte) oplossing van het randwaardeprobleem (2.1) zijn. Schrijf uV voor de projectie van u op de deelruimte V ten opzichte van het inproduct [·, ·].
Er geldt dat uV de unieke oplossing van de Galerkin methode (2.2) is.
Bewijs. Er geldt dat u − uV ⊥ V , aangezien uV de projectie van u op de deelruimte V
is. Voor elke v ∈ V geldt er dus [u − uV, v] = 0, wat impliceert dat [u, v] = [uV, v]. Uit
Stelling 2.16 weten we dat u de zwakke formulering oplost, dus vinden we [uV, v] = [u, v] = (f, v) voor alle v ∈ V.
Hiermee is de existentie van het gediscretiseerde probleem (2.2) aangetoond.
Daar V een eindig dimensionale lineaire deelruimte is, kunnen we naar het ge¨ınduceerde inproduct [·, ·]V kijken. De uniciteit van uV volgt nu op eenzelfde manier als in het bewijs
van Stelling 2.17.
De oplossing van het gediscretiseerde probleem (2.2) is dus gelijk aan de projectie van de exacte oplossing op de deelruimte V . De crux van de Galerkin methode is nu dat we deze projectie uV kunnen berekenen zonder dat we de exacte oplossing u hoeven te
weten.
Om dit in te zien schrijven we uV als lineaire combinatie van de basisfuncties,
uV = n
X
i=1
αiφi.
Omdat uV de projectie is van u op V ten opzichte van [·, ·], geldt er in het bijzonder dat
u − uV ⊥ φj voor j ∈ {1, . . . , n}. Schrijven we dit uit dan krijgen we,
0 = [u − uV, φj] =⇒ n
X
i=1
αi[φi, φj] = [uV, φj] = [u, φj] = (f, φj) voor j ∈ {1, . . . , n}.
In matrix notatie is dit gelijk aan,
Aα := [φ1, φ1] · · · [φn, φ1] .. . . .. ... [φ1, φn] · · · [φn, φn] α1 .. . αn = (f, φ1) .. . (f, φn) . (2.3)
Als we de inproducten [φi, φk] en (f, φi) exact weten, dan kunnen we de de waarden
van αi bepalen . Daarmee kunnen we vervolgens uV, de projectie van u op V , bepalen
2.2.4. De eindige elementen methode
De eindige elementen methode is een speciaal geval van de Galerkin methode. In de eindige elementen methode wordt de deelruimte V ⊆ H01(I) gekozen als een ruimte van continue stuksgewijs lineaire functies ten opzichte van een partitie van het domein I. Definitie 2.19. Laat n+1 punten x0, . . . , xn gegeven zijn met a = x0 < x1 < · · · <
xn−1 < xn= b. Als J de verzameling van de intervallen Ij := [xj, xj+1] is, dan heet J
een partitie van I.
De grootte van de intervallen noteren we met hj = xj+1− xj en de fijnheid van de
partitie wordt gegeven door h = max{h0, . . . , hn−1}.
De deelruimte V van H01(I) die we nu gaan bekijken, is de ruimte van continue stuks-gewijs lineaire functies ten opzichte van een partitie J .
Definitie 2.20. De verzameling van continue stuksgewijs lineaire functies ten opzichte van J defini¨eren we met
S(J ) = {v ∈ C0(I) | ∀J ∈ J , v|J ∈ P1(J )}.
Elke functie v ∈ S(J ) ligt vast door zijn waarden in de punten xj. Hierdoor is de
ruimte dus een (n + 1)-dimensionale deelruimte van C0(I). Een (prettige) basis voor de ruimte S(J ) is de zogenaamde nodale basis.
Definitie 2.21. Laat ψj het unieke element uit S(J ) zijn met ψj(xi) = δij voor i ∈
{0, . . . , n}. De nodale basis van S(J ) wordt nu gegeven door B = {ψ0, . . . , ψn}.
Voor v ∈ S(J ) hebben de eigenschap dat
v = a0ψ0+ · · · + anψn ⇐⇒ aj = v(xj) voor alle j ∈ 0, . . . , n.
Uit deze eigenschap volgt direct dat B daadwerkelijk een basis van S(J ) is. Zie Figuur 2.1 voor een voorbeeld van nodale basisfuncties.
Met S(J ) vinden we een deelverzameling van H1
0(I) door te kijken naar
S0(J ) = S(J ) ∩ C01(I).
Dit is de ruimte van continue stuksgewijs lineaire functies die tevens nul op de rand zijn. Als een functie uit S(J ) nul op de rand is, dan moet deze functie de waarde nul aannemen in de punten x0 en xn. De nodale basis voor deze ruimte wordt dus
gegeven door B0 = {ψ1, . . . , ψn−1}. De eindige elementen oplossing wordt nu gegeven
als de projectie van u op S0(J ). Met de methode die in de vorige sectie beschreven werd
kunnen we deze projectie bepalen. Dit doen we door het volgende matrixstelsel op te lossen, Aα := [ψ1, ψ1] · · · [ψn−1, ψ1] .. . . .. ... [ψ1, ψn−1] · · · [ψn−1, ψn−1] α1 .. . αn−1 = (f, ψ1) .. . (f, ψn−1) (2.4)
0.0 0.2 0.3 0.7 0.9 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Nodale basis functie ψ0
Nodale basis functie ψ2
Nodale basis functie ψ3
Figuur 2.1.: Een voorbeeld van 3 nodale basisfuncties op een partitie van [0, 1].
2.2.5. De volledig discrete eindige elementen methode
Om de eindige elementen oplossing te bepalen moeten de inproducten (f, ψj) exact
be-rekend worden, dat kan lastig zijn. Bij de volledig discrete eindige elementen methode wordt de berekening van deze inproducten ook gediscretiseerd. Dit kan door f te ver-vangen door zijn nodale interpolant Πf , gedefinieerd door
Π : C(I) → S(J ) : w 7→ Πw =
n
X
i=0
w(xj)ψj.
De volledig discrete eindige elementen methode is nu gedefinieerd als het vinden van uΠV ∈ S0(J ) zodat
[uΠV, v] = (Πf, v) voor alle v ∈ S0(J ). (2.5)
In de volgende sectie geven we een bovengrens voor de fout tussen de exacte oplossing u en de volledig discrete eindige elementen oplossing uΠV. De volledig discrete eindige elementen methode is eigenlijk een herformulering van de eindige elementen methode, met Πf in plaats van f . Dit reduceert daarom weer tot het oplossen van een matrixver-gelijking zoals in (2.4). In dit geval vinden we voor de rechterkant van deze vermatrixver-gelijking
(Πf, ψj) = n
X
i=0
f (xi)(ψi, ψj) voor j = 1, . . . , n − 1.
In matrix notatie wordt de inproductvector dus gegeven door (Πf, ψ1) .. . (Πf, ψm−1) = (ψ0, ψ1) · · · (ψn, ψ1) .. . . .. ... (ψ0, ψn−1) · · · (ψn, ψn−1) f (x0) .. . f (xn)
Voor het bepalen van de volledig discrete eindige elementen oplossing moeten we dus het volgend systeem oplossen,
[ψ1, ψ1] · · · [ψn−1, ψ1] .. . . .. ... [ψ1, ψn−1] · · · [ψn−1, ψn−1] α1 .. . αn−1 = (ψ0, ψ1) · · · (ψn, ψ1) .. . . .. ... (ψ0, ψn−1) · · · (ψn, ψn−1) f (x0) .. . f (xn) . In het kort noteren we dit als,
Aα = M ˆf . (2.6)
Hier wordt A de stijfheidsmatrix genoemd en M de massamatrix. Merk op dat de dimensies van deze matrices niet hetzelfde zijn. Dit komt doordat de nodale interpolant Πf niet gelijk aan nul hoeft te zijn op de rand van het interval.
2.2.6. Bovengrens fout eindige elementen oplossing
We kunnen een aantal bovengrenzen geven voor de fout van eindige elementen oplossing. Definieer hiertoe de energienorm op H01(I) met,
kvkA=p[v, v] voor alle v ∈ H01(I).
Volgens Lemma 2.18 is de oplossing uV van de Galerkin methode gelijk aan de projectie
van u op de deelruimte V ten opzichte van [·, ·]. In de energienorm betekent dit dat uV
de beste benadering van u in V is, oftewel
ku − uVkA≤ ku − vkAvoor alle v ∈ V. (2.7) In de eindige elementen methode namen we de ruimte V gelijk aan S0(J ). Als we
opmerken dat Πu ∈ S0(J ), dan vinden we de volgende bovengrens
ku − uVkA≤ ku − ΠukA= O(h).
Hierin is h de fijnheid van de partitie J zoals gedefinieerd in Definitie 2.19. De eerste afschatting volgt direct uit (2.7). De tweede afschatting volgt uit de approximatietheorie, zie bijvoorbeeld [17, Theorem 3.3].
We kunnen de energienorm relateren aan de standaardnorm op C01(I) middels het volgende lemma dat ook wel de ongelijkheid van Poincar´e-Friedrich wordt genoemd [3, Thm 1.5].
Lemma 2.22. Voor alle v ∈ C1
0(I) hebben we dat
kvk2≤ 1 2
√
Bewijs. Aangezien v(a) = 0 volgt uit de hoofdstelling van de integraalrekening dat v(x)2= Z x a v0(t) dt 2 ≤ Z x a 12dt Z x a (v0(t))2dt Cauchy-Schwarz op C[a, x] ≤ (x − a) Z b a (v0(t))2dt = (x − a)kvk2A.
Als we dit invullen in de standaard norm van v vinden we ||v||2 = s Z b a v(x)2dx ≤ s Z b a (x − a)kvk2Adx = kvkA 1 2 √ 2(b − a)
Voor de volledig discrete eindige elementen oplossing vinden we nu de volgende bo-vengrens.
Stelling 2.23. Laat uV de eindige elementen oplossing zijn en uΠV de volledig discrete
eindige elementen oplossing, dan hebben we kuV − uΠVkA≤
1 2
√
2(b − a)||f − Πf ||2 = O(h).
Bewijs. We hebben voor alle v ∈ V dat
[uV, v] = (f, v) en [uΠV, v] = (Πf, v).
Dit impliceert dat [uV − uΠV, v] = (f − Πf, v) voor alle v ∈ V . In het bijzonder geldt nu
voor v = uV − uΠV ∈ V dat kuV − uΠ Vk2A= [uV − uΠV, uV − uΠV] = (f − Πf, uV − uΠV) ≤ ||f − Πf ||2||uV − uΠ V||2 Cauchy-Schwarz ≤ ||f − Πf ||21 2 √ 2(b − a)kuV − uΠVkA Lemma 2.22
Delen door kuV − uΠVkA levert nu het resultaat (als kuV − uΠVkA = 0 valt er niets te
bewijzen).
De afschatting kf − Πf k2 = O(h) volgt uit de approximatietheorie, zie bijvoorbeeld
[17, Theorem 3.3].
Combineren we deze stelling met de eerdere afschatting ku−uVkA= O(h), dan vinden
we
2.2.7. Voorbeeld
We zijn nu in staat een implementatie in Python te geven van de volledig discrete eindige elementen methode. Neem I = [0, 1] en f (x) = 2, dan zien we eenvoudig in dat de exacte oplossing van (2.1) wordt gegeven door u(x) = −(x − 1)x. Neem partitie J = (0, 0.25, 0.5, 0.75, 1). Met de code uit Appendix A.1 kunnen we de volledig discrete eindige elementen oplossing uΠV berekenen voor de deelruimte S0(J ). Het resultaat kun
je zien in onderstaande Figuur 2.2. De implementatiedetails bekijken we in Sectie 3.1.
0.00 0.25 0.50 0.75 1.00 0.00 0.05 0.10 0.15 0.20 0.25
Eindige elementen oplossing Exacte oplossing
Figuur 2.2.: Volledig discrete eindige elementen benadering van u(x) = −(x − 1)x, de oplossing van het eendimensionale randwaardeprobleem (2.1).
2.3. Het reactie-diffusieprobleem
In de voorgaande secties hebben we de theorie achter de eindige elementen methode gezien voor het (speciale) eendimensionale diffusieprobleem (2.1). Het idee achter de bovenstaande theorie kun je echter in een algemenere vorm toepassen op elliptische differentiaalvergelijkingen [3, Chapter 1]. De elliptische differentiaalvergelijkingen die wij zullen bekijken zijn de reactie-diffusieproblemen.
Gegeven d ≥ 1, laat Ω ⊆ Rd een begrensd polytoop1 zijn met Lipschitz rand2 ∂Ω. Zij verder f, g ∈ C(Ω) gegeven met g ≥ 0. Het reactie-diffusieprobleem is het vinden van een ug ∈ C2(Ω) zodat
−∆ug+ gug = f in Ω
ug = 0 op ∂Ω (2.8)
Hierin is ∆ de Laplaciaan uit Definitie 2.10. De term gug wordt de reactieterm van het probleem genoemd. Als we d = 1 en g = 0 kiezen dan zien we dat het eendimensionaal randwaardeprobleem (2.1) een instantie van het reactie-diffusieprobleem is.
We zullen nu de stellingen en definities uit Sectie 2.2 generaliseren, zodat we de eindige elementen methode kunnen toepassen op dit algemenere probleem.
2.3.1. De zwakke formulering van het reactie-diffusieprobleem
Net als het eendimensionale geval kunnen we het reactie-diffusieprobleem (2.8) in een zwakke formulering geven. Schrijf daartoe H1(Ω) voor de Sobolev ruimte met (meer)di-mensionale functies die nu op Ω leven. Verder noteren we weer H01(Ω) voor de functies die tevens nul zijn op de rand, H1
0(Ω) = {v ∈ H1(Ω) : v|∂Ω = 0}. In het
eendimensi-onale geval hebben we een inproduct [·, ·] op de Sobolev ruimte gedefinieerd. Voor het algemenere geval kunnen we ook een dergelijk inproduct defini¨eren.
Definitie 2.24. Het inproduct dat geassocieerd is met het reactie-diffusieprobleem (2.8) wordt gegeven door
a(g; ·, ·) : H01(Ω) × H01(Ω) → R : a(g; w, v) = (∇w, ∇v) + (gw, v).
Hier is (·, ·) het standaard inproduct op L2(Ω) en ∇ de gradi¨ent operator, zie Definitie 2.8. Dat a(g; ·, ·) daadwerkelijk een symmetrische bilineaire vorm is, is makkelijk na te gaan. Voor het bewijs dat het ook een inproduct definieert verwijzen we naar de literatuur, zie bijvoorbeeld [4, Lemma 4.4].
De zwakke formulering van het reactie-diffusieprobleem (2.8) is nu het vinden van een ug ∈ H1
0(Ω) zodat
a(g; ug, v) = (f, v) voor alle v ∈ H01(Ω). (2.9) In dit algemenere geval vinden we de volgende verbanden tussen een ug en de oplossing van zwakke formulering.
1
Veralgemenisering van het begrip polygoon.
Stelling 2.25. Zij ug een oplossing van het reactie-diffusieprobleem (2.8), dan is ug ook een oplossing van de zwakke formulering (2.9).
Bewijs. Volgens de divergentiestelling3 geldt voor alle v ∈ H1
0(Ω) dat Z Ω div(v∇ug) dS = Z ∂Ω v∇u · da.
Het rechterlid van deze vergelijking is gelijk aan nul, want v is nul op de rand ∂Ω. Als we het linkerlid uitschrijven met behulp van de productregel voor de divergentie, dan vinden we Z Ω div(v∇ug) dS = Z Ω ∇v · ∇ug+ v div(∇ug) dS. Aangezien het linkerlid nul is en div(∇ug) = ∆ug geldt er nu dat
−(∆ug, v) = (∇ug, ∇v). (2.10)
Anderzijds is ug een oplossing van het reactie-diffusieprobleem (2.8), dus er geldt dat −∆ug+ gug = f in Ω. Maar dit impliceert dat −(∆ug, v) + (gug, v) = (f, v) voor alle
v ∈ H01(Ω). Gebruiken we nu de gelijkheid (2.10), dan vinden we dat −(∇ug, ∇v) + (gug, v) = (f, v) voor alle v ∈ H01(Ω). Oftewel ug is een oplossing van de zwakke formulering (2.9).
Stelling 2.26. De zwakke formulering (2.9) heeft een unieke oplossing.
Bewijs. Dit volgt uit het feit dat a(g; ·, ·) een inproduct is. Het bewijs gaat op eenzelfde manier als het bewijs van Stelling 2.17.
Opmerking. Uit de twee bovenstaande stellingen volgt nu direct dat het reactie-diffusie-probleem (2.8) een unieke oplossing heeft.
2.3.2. De Galerkin methode
De zwakke formulering van het reactie-diffusieprobleem (2.9) kunnen we ook voor een deelruimte V ⊆ H01(Ω) geven. De Galerkin methode is dan het vinden van een ugV ∈ V zodat
a(g; ugV, v) = (f, v) voor alle v ∈ V. (2.11) Net als in Sectie 2.2 vinden we het volgende resultaat. Voor een bewijs verwijzen we naar de literatuur, bijvoorbeeld [4, Proposition 4.1].
Stelling 2.27. Zij ug de unieke oplossing van het reactie-diffusieprobleem (2.8). De
unieke oplossing van de Galerkin methode (2.11) is gelijk aan de projectie van ug op de deelverzameling V in het inproduct a(g; ·, ·).
3
Zie een dictaat over verctorcalculus, bijvoorbeeld http://mathworld.wolfram.com/ DivergenceTheorem.html.
Stel dat een basis voor V wordt gegeven door de n onafhankelijke functies φ1, . . . , φn.
Analoog aan het eendimensionale geval kunnen we de oplossing van de Galerkin methode schrijven als ugV = Pn
i=1αiφi, waarbij de co¨efficienten αi gevonden worden door het
oplossen van het matrixstelsel:
Aα := a(g; φ1, φ1) · · · a(g; φn, φ1) .. . . .. ... a(g; φ1, φn) · · · a(g; φn, φn) α1 .. . αn = (f, φ1) .. . (f, φn) . (2.12)
2.3.3. De eindige elementen methode
De eindige elementen oplossing van het reactie-diffusieprobleem verkrijgen we nu door V gelijk te nemen aan de verzameling van continue stuksgewijs lineaire functies op een partitie van Ω. De partitie deelt Ω op in een eindig aantal subdomeinen, elementen. De opdeling van een interval [a, b] altijd bestaat uit (sub)intervallen. In hogere dimensies hebben we meer keuze voor het type elementen, zo kunnen we een tweedimensionaal domein bijvoorbeeld opdelen in vierkanten of driehoeken. We zullen kijken naar een opdeling in simplices, de generalisatie van een driehoek/tetra¨eder.
Definitie 2.28. Een k-simplex is een k-dimensionaal polytoop, verkregen uit het con-vexe omhulsel van k+1 affien onafhankelijke4punten. De (k−1)−dimensionale zijvlakken worden de facetten genoemd.
Om continue stuksgewijs lineaire functies te defini¨eren op een opdeling in simplices moet de partitie van Ω aan bepaalde eisen voldoen.
Definitie 2.29. Laat T = {T1, . . . , TK} een partitie van Ω ⊂ Rd zijn in simplices Ti.
De partitie T heet face-to-face of een conforme triangulatie als aan de volgende eisen is voldaan:
• Ω = ∪K i=1Ti.
• Het inwendige van de simplices Ti is paarsgewijs disjunct.
• (face-to-face) Elke facet van een simplex uit de partitie ligt ´of op de rand ∂Ω ´of is wederom de facet van een andere simplex uit T .
Met een triangulatie bedoelen we altijd een conforme triangulatie, tenzij anders vermeld. Een triangulatie heeft ook een waarde h. Deze waarde is gerelateerd aan de geometri-sche grootte van de simplices en geeft een maat voor de ‘fijnheid’ van de partitie. Definitie 2.30. Laat T een triangulatie zijn. De grootte hi van een simplex Ti ∈ T
wordt gegeven door zijn (euclidische) diameter, hi = max
x,y∈Ti
kx − yk2.
De waarde h die de fijnheid van de triangulatie T aangeeft, wordt gegeven door h = max{h1, . . . , hK}.
Laat T nu een triangulatie zijn van het domein Ω en noteer de vertices van T met x1, . . . , xm. Order de vertices zodat x1, . . . , xn in Ω liggen en xn+1, . . . , xm op de rand
∂Ω liggen. We bekijken nu de ruimte S(T ) ⊂ H1(Ω) van continue stuksgewijs line-aire functies ten opzichte van T . Voor deze ruimte kunnen we weer een nodale basis defini¨eren.
Definitie 2.31 (Nodale basis). Laat ψjhet unieke element uit S(T ) zijn met ψj(xi) = δij
voor i ∈ {1, . . . , m}. De nodale basis van S(T ) wordt nu gegeven door B = {ψ1, . . . , ψm}.
De nodale basis van S0(T ) := S(T ) ∩ H01(Ω) wordt gegeven door B0= {ψ1, . . . , ψn}.
Als we in de Galerkin methode de deelruimte S0(T ) gebruiken, dan noemen we dit
de eindige elementen methode. De co¨efficienten van de eindige elementen oplossing, ugV =Pm
i=1αiφi, vinden we door het matrixstelsel (2.12) op te lossen met de basisfuncties
ψ1, . . . , ψn.
2.3.4. De volledig discrete eindige elementen methode
Compleet analoog aan het eendimensionale geval kunnen we nu ook de volledig discrete eindige methode bekijken. Definieer daartoe de nodale interpolatie operator met
Π : C(Ω) → S(T ) : w 7→ Πw =
m
X
i=1
w(xi)ψi.
De volledig discrete methode bestaat dan uit het vinden van een uΠgV zodat a(g; uΠgV , v) = (Πf, v) voor alle v ∈ S0(T ).
De inproducten (Πf, ψi) kunnen we nu gemakkelijk uitrekenen,
(Πf, ψj) = m
X
i=1
f (xi)(ψi, ψj) voor j = 1, . . . , n.
Het vinden van de volledig discrete eindige elementen oplossing reduceert zich dus tot het oplossen van het volgende matrixstelsel,
a(g; φ1, φ1) · · · a(g; φ1, φn) .. . . .. ... a(g; φn, φ1) · · · a(g; φn, φn) α1 .. . αn = (φ1, φ1) · · · (φ1, φm) .. . . .. ... (φn, φ1) · · · (φn, φm) f (x1) .. . f (xm) , (2.13) of kortweg, Aα = M ˆf .
De n × n-matrix A heet de stijfheidsmatrix en de n × m-matrix M noemen we de massamatrix.
2.3.5. Bovengrens fout eindige elementen oplossing
De energienorm voor het reactie-diffusieprobleem wordt gegeven door, kvkA=pa(g; v, v) voor alle v ∈ H01(Ω).
De waarde h van een triangulatie T geeft de fijnheid aan. We kunnen hiermee de volgende bovengrens vinden,
kug− ugVkA= O(h).
Voor een bewijs verwijzen we naar de literatuur, bijvoorbeeld [3, Chapter 7]. We kunnen ook een bovengrens voor de fout tussen de eindige elementen oplossing en de volledig discrete eindige elementen oplossing geven [5, Proposition 2.1],
kugV − uΠgV kA= O(h).
Combineren we deze twee bovengrenzen dan vinden we de volgende afschatting, kug− uΠgV kA= O(h).
3.
Implementatie van de volledig
discrete eindige elementen methode
In het vorige hoofdstuk zagen we de theorie achter de eindig elementen methode. Het vinden van de eindig elementen oplossing voor het reactie-diffusieprobleem (2.8) redu-ceert zich tot het oplossen van de matrix vergelijking
Aα := a(g; φ1, φ1) · · · a(g; φn, φ1) .. . . .. ... a(g; φ1, φn) · · · a(g; φn, φn) α1 .. . αn = (f, φ1) .. . (f, φn) .
In dit hoofdstuk zullen we bekijken hoe we dit stelsel (handig) numeriek kunnen oplos-sen. Voor de implementatie zullen we gebruik maken van Python1 met de bibliotheek NumPy2.
3.1. Het eendimensionale geval
We bekijken eerst weer het eendimensionale diffusieprobleem (2.1). Laat J een partitie van het interval I zijn met S(J ) = span{ψ0, . . . , ψn}. Voor het bepalen van de eindige
elementen oplossing moeten we de stijfheidsmatrix A uit (2.4) berekenen. Schrijf daartoe Ij = [xj−1, xj] en hj = xj − xj−1, zoal in Definitie 2.19.
Bedenk dat ψj de unieke functie uit S(J ) is zodat ψj(xi) = δij. Een voorbeeld van
zo een functie is gegeven in Figuur 2.1. We kunnen deze functies ook middels expliciete formules geven voor j ∈ {1, . . . , n − 1},
ψj(x) = x − xj−1 xj − xj−1 = x − xj−1 hj x ∈ Ij xj+1− x xj+1− xj = xj+1− x hj+1 x ∈ Ij+1 0 anders . (3.1)
De stijfheidsmatrix A bestaat uit de inproducten [ψi, ψj] = (ψi0, ψj0) met i, j ∈ {1, . . . , n−
1}. Over deze inproducten kunnen we het volgende opmerken:
1Programmeertaal, zie https://www.python.org/ voor meer informatie. 2
Een bibliotheek die ondersteuning biedt aan matrices met lineaire algebra met lineaire algebra func-tionaliteiten, zie ook www.numpy.org/.
• Voor j = 1, . . . , n − 1 is het inproduct [ψj, ψj] gelijk aan de integraal van (ψj0)2
over Ij∪ Ij+1. Op de rest van het domein is ψj immers constant nul.
• Voor j = 1, . . . , n − 1 is het inproduct [ψj−1, ψj] gelijk aan de integraal van ψj−10 ψj0
over Ij. Dit volgt direct uit het functievoorschrift (3.1).
• Voor |i − j| ≥ 2, geldt dat [ψi, ψj] = 0, ook dit volgt direct uit (3.1).
• Het inproduct is symmetrisch.
Uit deze opmerking volgt dat alleen de elementen [ψj, ψj] en [ψj−1, ψj] niet triviaal zijn.
Voor j = 1, . . . , n − 1 geldt er [ψj, ψj] = Z Ij∪Ij+1 (ψj0)2 = Z Ij 1 h2j + Z Ij+1 1 h2j+1 = 1 hj + 1 hj+1 en voor j = 2, . . . , n − 1 hebben we [ψj−1, ψj] = Z Ij ψj−10 ψj0 = Z Ij −1 hj 1 hj = − 1 hj . We kunnen nu de gehele stijfheidsmatrix A geven:
A = 1 h1 + 1 h2 − 1 h2 0 · · · 0 −h1 2 1 h2 + 1 h3 1 h2 . .. . .. .. . 0 . .. . .. . .. . .. 0 .. . . .. . .. −h1 n−2 1 hn−2 + 1 hn−1 − 1 hn−1 0 · · · 0 − 1 hn−1 1 hn−1 + 1 hn .
3.1.1. Discretisatie van de inproducten (f, ψj)
Bij de volledig discrete eindige elementen methode vereenvoudigen we de inproducten (f, ψj) door f te vervangen voor zijn nodale interpolant Πf . De inproducten (Πf, ψj)
kunnen we bepalen met behulp van de massamatrix M uit (2.6). De elementen van M worden gegeven door de inproducten (ψi, ψj), die we met behulp van soortgelijke
opmerkingen als hierboven expliciet kunnen berekenen. Als we dit doen, dan vinden we dat de massamatrix gelijk is aan
M = 1 6 2(h1+ h2) h2 0 · · · 0 h2 2(h2+ h3) h2 . .. ... 0 . .. . .. . .. 0 .. . . .. hn−2 2(hn−2+ hn−1) hn−1 0 · · · 0 hn−1 2(hn−1+ hn) .
Voor een partitie J kunnen we nu de stijfheidsmatrix A en de massamatrix M bepa-len. Als we vervolgens het stelsel (2.6) oplossen, dan wordt de volledig discrete eindige elementen oplossing gegeven door α1ψ1+ · · · + αn−1ψn−1. In Appendix A.1 is een
im-plementatie in Python gegeven. De matrices A en M worden in deze imim-plementatie iteratief opgebouwd, de logica hierachter zal in de volgende sectie duidelijk worden.
3.2. Het tweedimensionale reactie-diffusieprobleem
We bekijken het reactie-diffusieprobleem (2.8) met als domein een polygoon Ω ⊆ R2. Neem eerst g = 0 zodat we een diffusieprobleem in twee veranderlijken krijgen. Laat T = {T1, . . . , Tn} een triangulatie van Ω zijn en schrijf verder x0, . . . , xN voor de vertices
in T . Deze triangulatie T wordt in de implementatie gegeven door:
• Een matrix N waarin kolom j gelijk is aan de co¨ordinaten van vertex xj ∈ R2,
• Een matrix T waarin kolom j gelijk is aan de getallen p, q, r als de driehoek Tj
hoekpunten xp, xq, xr heeft,
• Een vector G waarin positie j aangeeft of hoekpunt xj op de rand ligt of niet.
Voor het bepalen van de volledig discrete eindige elementen oplossing moeten we nu de stijfheidsmatrix A en de massamatrix M uit (2.13) bepalen. We hebben geen reactieterm, dus de inproducten uit A laten zich vereenvoudigen: a(g; ψi, ψj) = (∇ψi, ∇ψj).
De methode die in de vorige sectie werd gebruikt om de inproducten (ψi, ψj) te bepalen
veralgemeniseert niet direct naar het tweedimensionale geval. We kunnen immers niet op voorhand stellen dat (ψi, ψj) = 0 voor zekere vaste |i − j|, daar er in twee dimensies
geen natuurlijke ‘volgorde’ voor de basisfuncties ψi is.
3.2.1. Nodale basisfuncties
De deelruimtes die we bekijken in de eindige elementen methode zijn S(T ) en S0(T ).
Noteer m = dim(S0(T )) en m+l = dim(S(T ). Voor deze ruimtes gebruiken we de nodale
basis, die wordt gegeven door de unieke ψ1, . . . , ψm+l met ψi(xj) = δij. Zie Figuur 3.1
voor een voorbeeld van een nodale basisfunctie in twee dimensies.
We gaan nu de berekening van de inproducten (ψi, ψj) en (∇ψi, ∇ψj) vereenvoudigen.
Daartoe defini¨eren we de referentiedriehoek ˆT met hoekpunten (0, 0), (1, 0) en (0, 1) en de functies,
v1(x, y) = 1 − x − y, v2(x, y) = x, v3(x, y) = y,
met gradi¨enten
∇v1 = −1 −1 , ∇v2 = 1 0 , en ∇v3= 0 1 .
Op de driehoek ˆT geldt nu dat v1, v2 en v3 lineaire functies zijn met de waarde 1 in ´e´en
hoekpunt en waarde 0 in de andere twee hoekpunten. Merk op dat de nodale basisfuncties ψizijn samengesteld uit kopie¨en van getransformeerde functies v1, v2 en v3. Dit feit gaan
−1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 0.0 0.2 0.4 0.6 0.8 1.0
Figuur 3.1.: Voorbeeld van een nodale basisfunctie
3.2.2. Opdelen matrices
Een manier om de massamatrix M en stijfheidsmatrix A te bepalen, is door eerst de (m + l) × (m + l) matrices ˜A en ˜M voor de grotere ruimte S(T ) uit te rekenen, preciezer
˜
Aij = (∇ψi, ∇ψj) en ˜Mij = (ψi, ψj) voor i, j ∈ {1, . . . , l + m}.
Uit ˜A kunnen we vervolgens de stijfheidsmatrix A verkrijgen door de l rijen en kolommen te verwijderen die bij de vertices horen die op de rand ∂Ω liggen. Op eenzelfde wijze kunnen we de massamatrix M construeren uit ˜M .
Een handige manier om ˜M samen te stellen komt voort uit de volgende observatie, (ψi, ψj) = Z Ω ψiψj = n X k=1 Z Tk ψiψj = n X k=1 (ψi, ψj)Tk, waarbij (ψi, ψj)Tk := R
Tkψiψj. We kunnen de matrix ˜M dus schrijven als een sommatie,
˜ M = n X k=1 ˜ Mk= n X k=1 (ψ0, ψ0)Tk · · · (ψ0, ψ0)Tk .. . . .. ... (ψ0, ψn)Tk · · · (ψn, ψn)Tk .
Op dezelfde wijze kunnen we ˜A schrijven als sommatie, ˜ A = n X k=1 ˜ Ak= n X k=1 (∇ψ0, ∇ψ0)Tk · · · (∇ψ0, ∇ψ0)Tk .. . . .. ... (∇ψ0, ∇ψn)Tk · · · (∇ψn, ∇ψn)Tk .
Met de defini¨erende eigenschap van de nodale basis kunnen we het volgende opmerken over deze basisfuncties beperkt tot ´e´en driehoek uit de triangulatie.
Opmerking. Zij Tk ∈ T een driehoek met hoekpunten xp, xq, xr ∈ R2. Op de driehoek
Tk zijn alleen de basisfuncties ψp, ψq en ψr niet identiek nul.
Uit deze opmerking volgt dat de matrices ˜Mk en ˜Ak slechts 3 × 3 elementen hebben
die niet nul zijn. Voor de matrix ˜Mk zijn dit de elementen
Ek= (ψp, ψp)Tk (ψp, ψq)Tk (ψp, ψr)Tk (ψq, ψp)Tk (ψq, ψq)Tk (ψq, ψr)Tk (ψr, ψp)Tk (ψr, ψq)Tk (ψr, ψr)Tk . (3.2)
Deze elementen staan niet in deze volgorde in de matrix M˜k, maar op de posities
{p, q, r} × {p, q, r}. In de implementatie worden deze posities gegeven door de matrix T.
3.2.3. Transformatie nodale basisfuncties
De matrices ˜Ak en ˜Mk bestaan uit de inproducten beperkt tot Tk ∈ T . Met behulp
van de substitutiestelling voor integralen kunnen we de berekening van deze beperkte inproducten vereenvoudigen.
Stelling 3.1 (Substitutiestelling voor integralen). Laat U en V open deelverzamelingen van Rn zijn en laat F : V → U een C1-diffeomorfisme3 zijn. Dan geldt
Z U f (x) dx = Z F (V ) f (x) dx = Z V (f ◦ F )(y)|detDF (y)| dy Bewijs. Zie bijvoorbeeld [7, Theorem 6.6.1].
Om in te zien hoe we deze stelling kunnen gebruiken bekijken we een driehoek Tk∈ T .
Schrijf p, q en r voor de hoekpunten van deze driehoek en laat ψp, ψq, ψrde bijbehorende
nodale basisfuncties zijn. Er bestaat een C1-diffeomorfisme tussen de referentiedriehoek ˆ T en Tk, namelijk Fk : ˆT → Tk: x y 7→p1 p2 +q1− p1 r1− p1 q2− p2 r2− p2 | {z } Bk x y , 3Bijectie met f ∈ C1 en f−1∈ C1.
met een constante totale afgeleide, DFk(x) = q1− p1 r1− p1 q2− p2 r2− p2 = Bk. (3.3)
We zien dat Fkhoekpunten van ˆT naar hoekpunten van Tkstuurt. Als we de basisfuncties
ψp, ψq en ψr beperken tot Tk dan geldt er,
v1= ψp◦ Fk, v2 = ψq◦ Fk en v3= ψr◦ Fk. (3.4)
3.2.4. Samenstellen massamatrix
Als we de substitutiestelling gebruiken met Fk dan vinden we dat
(ψp, ψq)Tk = Z Tk ψp(x)ψq(x) dx = Z Fk( ˆT ) ψp(x)ψq(x) dx = Z ˆ T ψp(Fk(x))ψq(Fk(x))|detDFk(x)| dx Substitutiestelling = |detBk| Z ˆ T v1(x)v2(x) dx identiteit (3.4) = |detBk|(v1, v2)Tˆ.
Een dergelijke relatie kunnen we afleiden voor (ψi, ψj)Tk met i, j ∈ {p, q, r}. Met deze
relaties zijn we in staat een expliciete uitdrukking te geven voor de matrix Ek, de 3 ×
3-niet triviale elementen van ˜Mk. Er geldt namelijk,
Ek= |det(Bk)|EM, (3.5) waarbij, EM := (v1, v1)Tˆ (v1, v2)Tˆ (v1, v3)Tˆ (v2, v1)Tˆ (v2, v2)Tˆ (v2, v3)Tˆ (v3, v1)Tˆ (v3, v2)Tˆ (v3, v3)Tˆ = 1 24 2 1 1 1 2 1 1 1 2 .
De laatste gelijkheid kan gevonden worden door de inproducten (vi, vj)Tˆ expliciet uit te
rekenen. Een algoritme om de massamatrix ˜M samen te stellen vinden we nu als volgt. Data: De triangulatie: N, T, G.
Result: (m + l) × (m + l) massamatrix. Initialiseer de massamatrix ˜M op nul; foreach Tk∈ T do
Bepaal de indices {p, q, r} van de basisfuncties die niet triviaal zijn op Tk;
Vind de co¨ordinaten van deze indices en bepaal hiermee Bk;
Bereken de matrix Ek aan de hand van (3.5);
Tel de waarden uit Ek op bij ˜M op posities {p, q, r} × {p, q, r};
end
3.2.5. Samenstellen stijfheidsmatrix
De substitutiestelling kunnen we ook gebruiken om de inproducten uit ˜Akte
vereenvou-digen. Hiervoor merken we eerst op dat de relatie tussen de totale afgeleide D en de gradi¨ent ∇ gegeven wordt door ∇v = (Dv)>. In (3.4) vonden we,
v1= ψp◦ Fk, v2 = ψq◦ Fk en v3= ψr◦ Fk.
Gebruiken we de kettingregel en de relatie tussen D en ∇, dan vinden we ∇v1= ∇(ψp◦ Fk) = [D(ψp◦ Fk)]> = [(Dψp◦ Fk)(DFk)]> = (DFk)>(Dψp◦ Fk)> = (DFk)>(∇ψp◦ Fk) en evenzo, ∇v2= (DFk)>(∇ψq◦ Fk) en ∇v3 = (DFk)>(∇ψr◦ Fk).
Bekijk nu het inproduct (∇ψp, ∇ψq)Tk, dan vinden we met de substitutiestelling dat
(∇ψp, ∇ψq)Tk = Z Fk( ˆT ) (∇ψp(x))>(∇ψq(x)) dx = Z ˆ T (∇ψp◦ Fk(x))>(∇ψq◦ Fk(x))|detDFk(x)| dx = |detBk| Z ˆ T ((DFk)−>∇v1(x))>((DFk)−>(∇v2(x))) dx = |detBk| Z ˆ T (∇v1)>(DFk)−1(DFk)−>(∇v2) dx = |detBk|(∇v1)>(DFk)−1(DFk)−>(∇v2) Z ˆ T 1 dx = |detBk|(∇v1)>(Bk)−1(Bk)−>(∇v2) 1 2.
De vijfde gelijkheid volgt uit het feit dat ∇vi(x) en DFk(x) constant zijn. Voor (∇ψi, ∇ψj)
met i, j ∈ {p, q, r}, vinden we een soortgelijke relatie, (∇ψp, ∇ψp)Tk · · · (∇ψp, ∇ψr)Tk .. . . .. ... (∇ψr, ∇ψp)Tk · · · (∇ψr, ∇ψr)Tk = 1 2| det Bk| ∇v1> ∇v> 2 ∇v3> (Bk)−1(Bk)−>[∇v1 ∇v2 ∇v3]
Deze vergelijking geeft ons de 3×3 niet-triviale elementen uit ˜Ak. Met een aantal (kleine)
aanpassingen aan Algoritme 1 uit de vorige sectie kunnen we nu de stijfheidsmatrix ˜A bepalen.
3.2.6. Volledig discrete eindige elementen oplossing
Met de methodes uit de vorige secties kunnen we de matrices ˜A en ˜M bepalen. De stijfheidsmatrix A vinden we nu door de l kolommen en rijen uit ˜A te verwijderen die horen bij de vertices die op de rand ∂Ω liggen. De massamatrix M vinden we door de l rijen uit ˜M te verwijderen die horen bij de vertices die op de rand ∂Ω liggen. De co¨efficienten αi van de volledig discrete eindige elementen oplossing,Pmi=1αiψi, worden
vervolgens gegeven door het matrixstelsel Aα = M ˆf uit (2.13) op te lossen.
3.2.7. Constante reactieterm
In de vorige secties zagen we een manier om de volledig discrete eindige elementen oplossing van een tweedimensionaal diffusieprobleem te berekenen. We kunnen deze me-thode eenvoudig uitbreiden naar een reactie-diffusieprobleem met constante reactieterm, g(x) = c ∈ R. Het inproduct a(g; ·, ·) laat zich dan immers schrijven als
a(g; v, w) = (∇v, ∇w) + (gv, w) = (∇v, ∇w) + c(v, w).
De bijbehorende stijfheidsmatrix Ag van de eindige elementen oplossing kunnen we dus schrijven als
Ag = A + c · Mg.
Hier is A de stijfheidsmatrix uit de vorige sectie en Mg is een deelmatrix van de massa-matrix ˜M . Preciezer, Mg is de matrix verkregen uit ˜M door de l kolommen en rijen te verwijderen die bij de vertices horen die op de rand ∂Ω liggen.
Met behulp van de vorige secties kunnen we nu de volledig discrete eindige elementen oplossing berekenen van het tweedimensionale reactie-diffusieprobleem met een constante reactieterm. Een implementatie in Python is te vinden in Appendix A.2.
3.2.8. De norm van een functie uit S(T ) berekenen
Neem de reactieterm nog steeds gelijk aan een constante c. Voor alle v ∈ S(T ) kunnen we nu eenvoudig de energienorm kvk2A= a(g; v, v) berekenen. Schrijf hiertoe v uit in de nodale basis, v = m+l X i=1 γiψi, met γ = γ1 .. . γm+l = v(x1) .. . v(xm+l) .
Als we deze sommatie in het inproduct a(g; v, v) invullen en de uitdrukking uitschrijven, dan vinden we
kvk2
A= a(g; v, v) = γ
>( ˜A + c ˜M )γ,
waarbij ˜A en ˜M de volledige (m + l) × (m + l) matrices zijn voor de ruimte S(T ). Op analoge wijze vinden we voor de standaardnorm dat
Als ug de exacte oplossing van het reactie-diffusieprobleem is en uΠgV de volledig dis-crete eindige elementen oplossing is, dan vinden we middels de driehoeksongelijkheid,
kug− uΠgV kA≤ kug− ΠugkA+ kΠug− uΠgV kA.
Uit de approximatietheorie volgt dat kug − Πugk
A = O(h). Aangezien Πug en uΠgV
beide elementen uit S(T ) zijn, kunnen we de norm kΠug − uΠgV kA op bovenstaande
manier uitrekenen. Als we hier vinden dat deze norm ook O(h) is, dan volgt dus dat kug− uΠg
V kA= O(h).
3.2.9. De uniforme verfijning van een partitie
Het idee achter de eindige elementen methode is dat de benadering ‘beter’ wordt als de partitie wordt verfijnd4. Een manier om een tweedimensionale partitie T te verfijnen is
middels uniforme verfijning.
Definitie 3.2 (Uniforme verfijning). Laat T een triangulatie van het domein Ω zijn. De uniforme verfijning T0 verkrijgen we nu door iedere driehoek Tk ∈ T op te delen in vier
congruente subdriehoeken. Preciezer, door een hoekpunt op het midden van iedere zijde van Tk te plaatsen en de drie lijnstukken die deze nieuwe hoekpunten verbinden toe te
voegen als zijdes.
Een voorbeeld van een uniforme verfijning wordt gegeven in onderstaand figuur.
−1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 Triangulatie (a) Triangulatie T −1.0 −0.5 0.0 0.5 1.0 −1.0 −0.5 0.0 0.5 1.0 Triangulatie
(b) Uniforme verfijning van T
Voor de implementatie is de volgende opmerking handig.
Opmerking. Een nieuw hoekpunt op het midden van een zijde −−→xixj ligt in het inwendige
van Ω dan en slechts dan als de zijde −−→xixj onderdeel is van twee verschillende driehoeken
in de triangulatie T .
De implementatie uit Appendix A.2 is in staat een triangulatie gegeven door N, T, G uniform te verfijnen.
3.2.10. Voorbeeld
We geven nu een voorbeeld van een tweedimensionaal reactie-diffusieprobleem. Laat Ω = [0, 1] × [0, 1] en f (x, y) = sin(πx) sin(πy)(2π2+ 1) zijn, bekijk het probleem van het vinden van een u ∈ C2(Ω) zodat,
−∆u + u = f in Ω u = 0 op ∂Ω. Merk op dat in dit geval de exacte oplossing u gelijk is aan
u(x, y) = sin(πx) sin(πy).
We bekijken nu de volledig discrete eindige elementen oplossing van dit probleem. Laat de beginpartitie T van [0, 1] × [0, 1] bestaan uit twee rechthoekige driehoeken die samen het vierkant vormen. Met de implementatie uit Appendix A.2 kunnen we deze partitie een aantal keer verfijnen en hierop vervolgens de volledige eindige elementen methode loslaten. De verfijnde partities met de eindige elementen oplossingen zijn weergegeven in Figuur 3.4. In Figuur 3.3 zien we dat de volledig discrete eindige elementen oplossing lijkt te convergeren naar de nodale interpolant van de exacte oplossing.
1 2 3 4 5 10-3 10-2 10-1 100 ||Πu−uΠ V ||A ||Πu−uΠ V ||2
Figuur 3.3.: Grafiek die het verschil tussen Πugen uΠV weergeeft. Op de x-as is het aantal verfijning weergegeven. De triangulaties uit Figuur 3.4 komen overeen met x = 2, 3, 4.
0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.60.8 1.0 0.00.1 0.20.3 0.40.5 0.6 0.7 0.8 0.9 uh 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Triangulatie 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.81.0 0.0 0.2 0.4 0.6 0.8 1.0 uh 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Triangulatie 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.60.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 uh 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Triangulatie
Figuur 3.4.: Links de eindige elementen oplossing, rechts de partitie van [0, 1] × [0, 1].
3.3. Het driedimensionale reactie-diffusieprobleem
In de vorige sectie zagen we hoe we de volledig discrete eindige elementen methode voor het tweedimensionale geval kunnen implementeren. Het driedimensionale probleem kan op analoge wijze worden ge¨ımplementeerd. We hebben in dit geval een domein Ω ⊂ R3 met een opdeling T = {T1, . . . , Tn} bestaande uit tetra¨eders. Noteer weer
m + l = dim S(T ) en m = dim S0(T ). De nodale basis van S(T ) wordt gegeven door
3.3.1. Samenstellen matrices
Evenals in het tweedimensionale geval berekenen we eerst de stijfheidsmatrix ˜A en de massamatrix ˜M voor S(T ). We vinden dat deze matrices te schrijven zijn als
˜ A = n X i=1 ˜ Ak en ˜M = n X i=1 ˜ Mk,
waarbij ˜Ak en ˜Mk de inproducten beperkt tot Tk bevatten. Een tetra¨eder heeft vier
hoekpunten, dus er zijn vier nodale basisfuncties niet triviaal op Tk. De matrices ˜Ak en
˜
Mk bestaan dus uit 4 × 4 elementen die niet nul zijn.
Laat ˆT de referentietetra¨eder zijn met hoekpunten (0, 0, 0), (1, 0, 0), (0, 1, 0) en (0, 0, 1). Op ˆT defini¨eren we de basisfuncties,
v1= 1 − x − y − z, v2= x, v3= y en v4 = z.
Schrijf {p, q, r, s} voor de indices van de hoekpunten van Tk, dan wordt een C1
-diffeomor-fisme tussen ˆT en Tk gegeven door
Fk: ˆT → Tk : x y z 7→ xp yp zp + xq− xp xr− xp xs− xp yq− yp yr− yp ys− yp zq− zp zr− zp zs− zp | {z } Bk x y z .
Met dit diffeomorfisme vinden we de volgende verbanden,
v1 = ψp◦ Fk, v2= ψq◦ Fk, v3 = ψr◦ Fk en v4= ψs◦ Fk.
Gebruiken we deze samenstellingen in combinatie met de substitutiestelling voor Fk dan
vinden we dat de 4 × 4 niet triviale elementen van ˜Mk worden gegeven door
1 120| det Bk| 2 1 1 1 1 2 1 1 1 1 2 1 1 1 1 2 .
Op eenzelfde wijze als voor het tweedimensionale geval kunnen we dergelijke samenstel-lingen ook vormen voor ∇ψi. We vinden dan dat de 4 × 4 niet triviale elementen van
˜
Ak worden gegeven door
1 6| det Bk| −1 −1 −1 1 0 0 0 1 0 0 0 1 (Bk)−1(Bk)−> −1 1 0 0 −1 0 1 0 −1 0 0 1 .
Met deze verbanden voor ˜Ak en ˜Mk kunnen we een implementatie geven voor de
volledig discrete eindige elementen methode in het driedimensionale geval. In Appendix A.3 is een dergelijke implementatie in Python gegeven.
3.3.2. De uniforme verfijning van een triangulatie
Een driedimensionale partitie kunnen we uniform verfijnen, in de literatuur heet dit ook wel red refinement [14]. Bij deze verfijning wordt elke tetra¨eder T ∈ T in 8 subtetra¨eders verdeeld. Dit gebeurt door een hoekpunt op het midden van iedere zijde van T toe te voegen. Met deze nieuwe vertices kunnen we nu direct vier subtetra¨eders vormen in de ‘punten’ van T , zie Figuur 3.5. In het midden van de tetra¨eder vormen de nieuwe vertices een octra¨eder. Deze octa¨eder kunnen we vervolgens opsplitsen in vier subtetra¨eders. Voor deze opsplitsing moet een diagonaal in de octa¨eder worden gekozen, het is dus geen unieke opdeling.
Figuur 3.5.: De tetra¨eder wordt gesplitst in vier subtetra¨eders en ´e´en octa¨eder. De octa¨eder wordt vervolgens langs een diagonaal gesplitst in vier subte-tra¨eders. Bron: [15].
In Appendix A.3 is deze verfijning ge¨ımplementeerd, zie Figuur 3.6. We kiezen een willekeurige diagonaal om de octa¨eder langs te splitsen. Verder gebruiken we bij de verfijning een functie die expliciet retourneert of een punt op de rand ∂Ω ligt. Dit kan vermeden worden door gebruik te maken van een datastructuur die bijhoudt welke zijdes er op de rand liggen.
Figuur 3.6.: Implementatie van de uniforme verfijning van een tetra¨eder. Links het or-gineel, rechts de verfijning.
3.3.3. Voorbeeld
We geven nu een voorbeeld van een driedimensionaal reactie-diffusieprobleem. Als do-mein nemen we Ω = [0, 1]3 en als functie f (x, y, z) = sin(πz) sin(πx) sin(πy)(3π2+ 1), bekijk nu het probleem van het vinden van een u ∈ C2(Ω) zodat,
−∆u + u = f in Ω u = 0 op ∂Ω. Merk op dat in dit geval de exacte oplossing u gelijk is aan
u(x, y, z) = sin(πx) sin(πy) sin(πz).
Als beginpartitie T nemen we een opdeling van de kubus [0, 1]3 in vijf tetra¨eders, zie Figuur 3.7. Met de implementatie uit Appendix A.3 kunnen we nu de partitie verfijnen en de volledig discrete eindige elementen oplossing berekenen. In Figuur 3.8 lijkt deze eindige elementen benadering te convergeren naar de nodale interpolant van de exacte oplossing.
2 3 4 10-2 10-1 100 ||Πu−uΠ V ||A ||Πu−uΠ V ||2
Figuur 3.8.: Fout tussen de eindige elementen benadering en de nodale interpolant van de exacte oplossing. Op de x-as staat het aantal uniforme verfijningen van de partitie. In de eerste twee verfijningen liggen er geen vertices in het inwendige van de kubus, de eindige elementen oplossing is dan dus triviaal.
4.
Het maximum principe
In Hoofdstuk 2 hebben we de eindige elementen methode bestudeerd voor het reactie-diffusieprobleem. In dit hoofdstuk zullen we eerst zien dat exacte oplossingen van dit probleem voldoen aan het zogenaamde maximum principe. Vervolgens bekijken we wan-neer de eindige elementen oplossing ook aan dit principe voldoet. Er zal blijken dat dit samenhangt met de dihedrale hoeken van de simplices in de triangulatie van het domein.
4.1. Het exacte maximum principe
In Hoofdstuk 2 keken we naar het reactie-diffusieprobleem (2.8). We kunnen ook kijken naar een algemenere versie van dit probleem. Hiervoor hebben we eerst de volgende definitie nodig.
Definitie 4.1. Een re¨ele n × n matrix M is positief-definiet als x>M x > 0 voor elke niet nulvector x ∈ Rn.
Laat een co¨efficient matrix A(x) = {aij(x)}ni,j=1 gegeven zijn met aij ∈ C1(Ω) en
A(x) symmetrisch1 positief-definiet. Laat verder f, g ∈ C(Ω) gegeven zijn met g ≥ 0. Het algemene reactie-diffusieprobleem is het vinden van een u ∈ C2(Ω) zodat
−div(A∇u) + gu = f in Ω
u = 0 op ∂Ω (4.1)
Als we voor A de identiteitsmatrix nemen, dan vinden we het reactie-diffusieprobleem (2.8) uit Hoofdstuk 2.
Definitie 4.2 (Maximum principe). Het algemene reactie-diffusieprobleem voldoet aan het maximum principe als,
f ≤ 0 =⇒ u ≤ 0 op Ω
Met andere woorden zegt dit dat u zijn maximum aanneemt op de rand ∂Ω.
Bestudeer nu eerst het eendimensionale geval van (4.1). Het domein bestaat dan uit een interval I = [a, b] ⊂ R en we zoeken een u ∈ C2(I) zodat
−(A(x)u0(x))0+ g(x)u(x) = f (x) voor a < x < b en u(x) = 0 op ∂I. (4.2) In de literatuur wordt dit ook wel een Sturm-Liouville randwaardeprobleem genoemd [2, p. 363].
1a
Lemma 4.3. Het Sturm-Liouville randwaardeprobleem (4.2) voldoet aan het maximum principe.
Bewijs. De functie u is continu, dus hij neemt zijn maximum aan op het interval [a, b]. We leiden nu een tegenspraak af. Stel dat f ≤ 0 en dat het maximum wordt aangenomen voor x0 ∈ (a, b), waarbij u(x0) > 0.
Omdat u continu is kunnen we ξ1, ξ2∈ [a, b] vinden met ξ1 < x0< ξ2 zodat er geldt
u(x0) ≥ u(x) > 0 voor x ∈ (ξ1, ξ2) en u(ξ1) = u(ξ2) = 0. (4.3)
Pas nu de middelwaardestelling [6, Theorem 29.3] toe op de functie u en de intervallen [ξ1, x0] en [x0, ξ0]. We vinden dan een x1∈ (ξ1, x0) en een x2 ∈ (x0, ξ2) zodat geldt
u0(x1) = u(x0) − u(ξ1) x0− ξ1 = u(x0) x0− ξ1 > 0 u0(x2) = u(ξ2) − u(x0) ξ2− x0 = −u(x0) ξ2− x0 < 0.
Aangezien A(x) positief definiet is geldt er dat A(x) > 0 en dus volgt er dat (Au0)(x1) > 0
en (Au0)(x2) < 0.
Passen we nu nogmaals de middelwaardestelling toe op de functie Au0 en het interval [x1, x2] (dit kan omdat A en u0differentieerbaar zijn). Dan is er een ξ ∈ (x2, x3) waarvoor
geldt dat (Au0)0(ξ) = Au 0(x 2) − Au0(x1) x2− x1 < 0. Merk op dat ξ ∈ (ξ1, ξ2), dus volgt uit (4.3) dat u(ξ) > 0.
Anderzijds is u de exacte oplossing van de differentiaalvergelijking (4.2), dus (Au0)0(ξ) = g(ξ)u(ξ) − f (ξ) ≥ 0.
De ongelijkheid volgt uit het feit dat f ≤ 0 en g(ξ)u(ξ) ≥ 0.
We hebben nu een tegenspraak gevonden voor (Au0)0(ξ). Er moet dus gelden dat u ≤ 0, oftewel het eendimensionale algemene reactie-diffusieprobleem (4.2) voldoet aan het maximum principe.
Stelling 4.4. Het algemene reactie-diffusieprobleem (2.8) voldoet aan het maximum principe.
Zie bijvoorbeeld [3, Theorem 2.1] voor een bewijs.
4.2. Het discrete maximum principe
We bekijken (weer) het reactie-diffusieprobleem (2.8). Schrijf ugvoor de exacte oplossing van dit probleem. In Hoofdstuk 2 hebben we gezien hoe we de bijbehorende eindige elementen oplossing ugV kunnen bepalen. We bestuderen nu of ugV aan het discrete maximum principe voldoet.
Definitie 4.5 (Discrete maximum principe). De eindige elementen benadering ugV vol-doet aan het discrete maximum principe als,
f ≤ 0 =⇒ ugV ≤ 0 op Ω.
Uit Stelling 4.4 volgt dat ug aan het maximum principe voldoet. Het blijkt dat dit
niet direct garandeert dat zijn eindige elementen benadering ugV hier ook aan voldoet. Hieronder staat een voorbeeld van een eindige elementen oplossing waar het maximum principe geschonden wordt. Het voorbeeld is overgenomen uit [5], voor een verklaring waarom het maximum principe geschonden wordt verwijzen we naar dit artikel.
Voorbeeld 4.6. Bekijk het volgende reactie-diffusieprobleem: we zoeken ugj zodat,
−∆ugj + g
jugj = fj in (0, 1)
ugj = 0 op {0, 1},
waarbij fj(x) = −2j(2x − 1)2 en gj(x) = 2j. Neem de partitie J = [0, 0.25, 0.5, 0.75, 1].
Passen we nu de eindige elementen methode toe voor j ∈ {0, . . . , 15}, dan verkrijgen we het Figuur 4.1. Hierin is duidelijk te zien dat het discrete maximum principe geschonden wordt in het punt x = 0.5. We hebben immers fj ≤ 0, maar we vinden ook eindige
elementen oplossingen met ugj(0.5) > 0.
0.00 0.25 0.50 0.75 1.00 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0.0 0.1 0.2
Figuur 4.1.: De lijnen geven de eindige elementen oplossing weer, hogere intensiteit blauw is voor hogere waarde van j.
4.2.1. Herformulering van het discrete maximum principe
Uit het bovenstaande voorbeeld blijkt dat het maximum principe niet direct behouden blijft bij de eindige elementen methode. Om te onderzoeken wanneer dit wel gebeurt
bekijken we nogmaals de methode. Laat hiertoe T de triangulatie van het domein Ω zijn met vertices x1, . . . , xn+m. De nodale basis voor S(T ) noteren we met ψ1, . . . , ψn+m.
Hernoem de vertices zodat ψ1, . . . , ψnde nodale basis van S0(T ) is. De eindige elementen
oplossing, ugV = Pn
i=1αiψi, kan nu gevonden worden door het matrixstelsel uit (2.12)
op te lossen:
Agα = F, waarbij Agij = a(g; ψi, ψj) en Fj = (f, ψj).
Om te te bepalen of de eindige elementen oplossing ugV aan het discrete maximum prin-cipe voldoet merken we het volgende op.
Opmerking. Voor de nodale basisfuncties geldt ψj ≥ 0, verder hebben we per aanname
dat f ≤ 0. Hieruit volgt dat Fj = (f, ψj) =
R
Ωf (x)ψj(x) dx ≤ 0. Oftewel we vinden dat
f ≤ 0 =⇒ F ≤ 0.2
Anderzijds vinden we, omdat ugV een lineaire combinatie van ψi is, dat ugV zijn
ex-trema aanneemt op de vertices van T . Het discrete maximum principe kunnen we nu herformuleren in lineair algebra¨ısche termen.
Lemma 4.7. De eindige elementen benadering ugV voldoet aan het discrete maximum principe als,
F ≤ 0 =⇒ α ≤ 0 waarbij α = (Ag)−1F. (4.4) Bewijs. De functie ugV neemt zijn extrema aan op de vertices van T . De waarden in deze punten worden gegeven door ugV(xi) = αi voor i ∈ {0, . . . , n}. De ongelijkheid α ≤ 0
impliceert dus dat ugV ≤ 0.
Met de aanname en de eerder gemaakte opmerking vinden we nu de volgende reeks implicaties:
f ≤ 0 =⇒ F ≤ 0 =⇒ α ≤ 0 =⇒ ugV ≤ 0.
4.2.2. Het discrete maximum principe en de Stieltjes matrix
Met bovenstaand lemma hebben we hebben we het discrete maximum principe geredu-ceerd tot een lineair algebra¨ısch probleem. Uit Lemma 4.7 en het feit dat α = (Ag)−1F vinden we dat
(Ag)−1 ≥ 0 (4.5)
een voldoende voorwaarde is voor het discrete maximum principe. Als we dus vinden dat (Ag)−1 ≥ 0, dan weten we dat ugV aan het discrete maximum principe voldoet.
We introduceren het begrip Stieltjes matrix om te bekijken wanneer de stijfheidsma-trix Ag de eigenschap (Ag)−1≥ 0 heeft.
Definitie 4.8 (Stieltjes matrix). Een re¨ele matrix (Mij)ni,j=1met Mij ≤ 0 voor alle i 6= j
is een Stieltjes matrix als M symmetrisch en positief-definiet is.
