overzicht   overview

Labrujère's probleem

In februari 1976 schreef T.E. Labrujère, toentertijd werkzaam bij het N.L.R. (Nationaal Lucht- en Ruimtevaartlaboratorium), een memorandum [NLR] met als titel: De "Eindige Elementen - Kleinste Kwadraten" Methode toegepast op de 2D Incompressibele Stroming om een Cirkel Cylinder. Om helemaal precies te zijn een onsamendrukbare en rotatievrije stroming. In de rapportage werd klaar en duidelijk vastgesteld dat een rechttoe rechtaan toepassing van de kleinste kwadraten eindige elementen methode, tegen alle verwachting in, niet goed blijkt te werken. Er spreekt een grote integriteit uit dit memorandum, hetgeen alle lof verdient. In de appendix "BASIC NLR" is alle benodigde programmatuur nagespeeld. Men dient deze programma's in de aangegeven volgorde te draaien, teneinde het probleem te reproduceren [FTP].

In december 1976 werd het probleem "opgelost" door G. de Vries en D.H. Norrie, beide verbonden aan de afdeling "Mechanical Engineering" van de universiteit van Calgary (in Canada). De samenvatting ("abstract") van het desbetreffende rapport [NdV] wordt hieronder weergegeven, door mij vertaald in het Nederlands:

De kleinste kwadraten eindige elementen methode wordt geformuleerd voor de twee-dimensionale rotatievrije stroming van een onsamendrukbare wrijvingsloze vloeistof. Vastgesteld wordt welke eisen van continuïteit aan de snelheids-componenten moeten worden opgelegd. Aangetoond wordt dat deze veel stringenter zijn dan tot nu toe werd aangenomen. Bij het ontwikkelen van een oplossingsprocedure is daarom uitgegaan van test-funkties van de vijfde orde voor beide componenten van de snelheid. Met gebruik van de kleinste kwadraten procedure wordt de stroming om een cirkel cylinder berekend. Getoond wordt dat de resultaten zeer nauwkeurig overeenstemmen met de theoretische oplossing. Einde vertaling.

Het is duidelijk dat bovenstaande "oplossing" van zuiver theoretisch belang is. De klaarblijkelijke noodzaak van vijfde orde interpolaties maakt de methode van Norrie en de Vries in de praktijk volledig onwerkbaar. De methode is voor dit eenvoudige geval al veel te ingewikkeld. Zodoende is een generalisatie naar minder geïdealiseerde situaties nauwelijks mogelijk. Uiteindelijk moeten 2-D en 3-D Navier-Stokes vergelijkingen worden opgelost, op een kromlijnig raster. Het vertrekpunt dient dus iets zijn wat vele malen eenvoudiger is. Met name het aantal onbekenden per knooppunt zou niet hoger moeten komen te liggen dan wat strikt nodig is: twee stuks.

Een onsamendrukbare rotatievrije stroming van een wrijvingsloze vloeistof wordt beschreven door het volgende stelsel PDV's: $$\begin{array}{ll}\Large \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} = 0 & \qquad \mbox{: onsamendrukbaar} \\ \Large \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} = 0 & \qquad \mbox{: rotatievrij} \end{array} $$ Hierin is $(x,y) =$ coördinaten, $(u,v) =$ snelheids-componenten.
Er bestaat voor dit stelsel vergelijkingen geen standaard Eindige Elementen methode, er is niet een of ander "natuurlijk" variatie-principe te bedenken, een procedure als in "Verenigde diffusie". De kleinste kwadraten methode vormt in zo'n geval een mogelijk alternatief. Kwadrateer de vergelijkingen zoals ze daar staan, tel de kwadraten bij elkaar op, integreer het geheel over het rekendomein, bepaal het minimum van de integraal als funktie van de onbekenden. Dit is de aanpak zoals beschreven in Zienkiewicz [OC] hoofdstuk 3.14.2. In ons geval: $$ \iint \left\{ \left[ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right]^2 + \left[ \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right]^2 \right\} \, dx.dy = \mbox{minimum} $$ Dit ziet er vrij eenvoudig uit, maar schijn bedriegt! Kleinste kwadraten is wel de meest verradelijke eindige elementen methode die ooit is uitgevonden. We hebben al gezien dat kleinste kwadraten niet werkt voor lineaire driehoeken. Mijn bedoeling is nu om te laten zien dat Labrujère's probleem ook "echt" kan worden opgelost. Daarmee bedoel ik: op een praktisch bruikbare, simpele manier. Daartoe moeten we echter naar het probleem kijken met een andere bril op: alsof het in wezen een Eindige Differentie, en geen eindige elementen probleem is. Het SaamhorigheidsBeginsel is in dit geval onontkoombaar.

Laten we de details bekijken. Om te beginnen wordt de integraal opgesplitst in element-bijdragen: $$ \sum_E \iint \left\{ \left[ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right]^2 + \left[ \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right]^2 \right\} \, dx.dy = \mbox{minimum} $$ Het is altijd voordelig om Numerieke in plaats van "exakte" integratie uit te voeren: zie [OC] hoofdstuk 8.8. Dit betekent dat men funktiewaarden dient te bepalen in zogenaamde integratiepunten $p$. Met ieder integratiepunt is in het algemeen een zekere gewichtsfaktor $w_p$ verbonden. $$ \sum_E \sum_p w_p \left\{ \left[ \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} \right]_p^2 + \left[ \frac{\partial v}{\partial x} - \frac{\partial u}{\partial y} \right]_p^2 \right\} \, J_p = \mbox{minimum} $$ Hierin is $J_p$ de Jacobiaan (determinant) van een coördinaten transformatie: zie de appendices "Driehoeks/VierkantsAlgebra". Met behulp van Differentiatie matrices kunnen de uitdrukkingen in de integrand worden uitgewerkt.

Wat nu volgt is een kleine stap voor een mens (en niet echt een grote stap voor de mensheid). We voegen de sommaties over de elementen en de integratiepunten samen to één grote globale sommatie over integratiepunten $(I=E,p)$. Dit betekent dat optellen van alle element-bijdragen, samen met het optellen van alle bijdragen op de integratiepunten, gelijkwaardig is met optellen over alle integratiepunten in het gehele rekendomein, ale het ware in één ruk. Integratiepunten kunnen zodoende worden opgevat als elementairder dan de elementen zelf. En elementen met meer dan één integratiepunt kunnen worden beschouwd als een superpositie van elementair geïntegreerde elementen, met slechts één integratiepunt in elk ervan. Dergelijke beschouwingen weken de mens alvast een beetje los van allerhande conventionele opvattingen. Want als elementen superposities zijn, dan liggen ze ook over elkaar heen, en sluiten ze dus ook niet netjes op hun grenzen tegen elkaar aan.

Om de kleinste kwadraten methode te laten werken, moet het te behalen minimum een klein getal zijn. Een getal dat bovendien naar nul zal konvergeren naarmate de elementen kleiner worden gemaakt. Het is misschien niet zo'n dom idee om dan maar meteen te eisen dat het beoogde minimum van het begin af aan nul moet zijn. Maar in het laatste geval zou de kleinste kwadraten "variatie integraal" equivalent zijn met het volgende on-gekwadrateerde systeem van vergelijkingen. Immers een som van kwadraten van termen is nul, dan en slechts dan als elk van de afzonderlijke termen nul is: $$ \frac{\partial}{\partial x} u + \frac{\partial}{\partial y} v = 0 \qquad \frac{\partial}{\partial x} v - \frac{\partial}{\partial y} u = 0 \qquad \mbox{: voor elk integratiepunt} $$ In operatorvorm geschreven om te benadrukken dat $\partial/\partial (x,y)$ differentiaties voorstellen in gediskretiseerde vorm: differentiatie matrices. We hebben hier te maken met een algebraisch en geen analytisch stel vergelijkingen.

Word u ervan bewust dat ieder "integratiepunt" in feite niets anders doet dan een tweetal vergelijkingen aandragen. Alle integratiepunten samen dragen ertoe bij dat er op deze manier een stelsel vergelijkingen ontstaat. Er is niets op tegen om dit stelsel vergelijkingen een Eindige Differentie systeem te noemen. Laat ons nu "integratiepunt" vervangen door "eindige differentie vergelijking", en we zijn aangeland waar we wezen moeten:

Iedere goed werkende kleinste kwadraten eindige elementen methode is geheel en al gelijkwaardig met het nul stellen van de som van de kwadraten van alle vergelijkingen in een eindige differentie systeem, en heeft dus ook dezelfde oplossing als dit stelsel eindige differentie vergelijkingen.

En passant wordt hiermee de oplossing nabij gebracht voor een ander probleem: de vaak slechte konditie van kleinste kwadraten vergelijkingen. Merk op dat in het eindige differentie stelsel de gewichtsfaktoren $w_p.J_p$ niet zijn terug te vinden. Men kan ze met evenveel recht vervangen door anderssoortige schaalfaktoren. Dit wordt reeds gesuggereerd in [OC] hoofdstuk 3.14.2: "Wederom zou deze gewichtsfunktie zo kunnen worden gekozen dat men verzekerd is van een konstante verhouding van de verschillende element-bijdragen". Als dit de bedoeling is, dan doet bijvoorbeeld de oppervlakte van een element bij kleinste kwadraten methoden in het geheel niet ter zake. Zienkiewicz vervolgt met: "ofschoon dit nog nimmer in praktijk is gebracht". Iemand moet tenslotte de eerste zijn. Wij hebben de daad bij het woord gevoegd in BASIC programma-onderdelen getiteld "SpoorWegen". Het komt erop neer dan iedere vergelijking wordt gedeeld door zijn "lengte": de wortel uit de som van de kwadraten van de coëfficienten. Dit blijkt hetzelfde te zijn als de kleinste kwadraten eindige elementen matrix wegen met zijn spoor. Vandaar de naam.

Laten we nu Labrujère's probleem nogmaals tegen het licht houden. Opdat het eindige elementen systeem een oplossing bezit, moet in bijbehorend systeem van eindige differentie vergelijkingen het aantal onbekenden $N$ gelijk zijn aan het aantal vergelijkingen $M$. Is dit niet het geval, dan noemen we het stelsel overbepaald, en dan is het nog maar de vraag of het kleinste kwadraten minimum (snel genoeg) naar nul gaat. Eenvoudig tellen van de driehoeken laat al vlug zien dat we hier inderdaad te maken hebben met een behoorlijk overbepaald stelsel. Het aantal elementen overtreft het aantal knooppunten met bijna een faktor twee. Dit betekent dat er ongeveer twee keer zoveel "ongekwadrateerde" vergelijkingen zijn dan onbekenden. Onafhankelijk van ingewikkelde argumenten zoals hogere orde continuïteit en dergelijke, hebben we hier een heel ander soort kwestie.

Norrie en de Vries hebben deze kwestie echter wel degelijk opgelost. Ze hielden om te beginnen vast aan de mesh van driehoeken. Als kompensatie voor het teveel aan elementen moesten zij heel wat extra variabelen toevoegen. Het is nu duidelijk welke andere benadering tot de mogelijkheden behoort. Er hoeft immers alleen maar een balans gevonden te worden tussen een aantal onbekenden en een aantal vergelijkingen: $M$ moet gelijk worden aan $N$. In plaats van het aantal onbekenden uit te breiden, kan men ook het aantal elementen verminderen. Een zuiver algebraisch argument, veel minder "geavanceerd" dan polynomen van de vijfde orde.

De kunst is alleen om het juiste element te vinden. De voorkeur gaat uit naar iets wat net zo eenvoudig is als driehoeken. Bestaat dat? Jazeker! Gerefereerd wordt naar de appendix "VierkantsAlgebra". Naast de lineaire driehoek bestaat er ook een lineaire vierhoek. Dit element heeft knooppunten op de middens van de zijden, en gaat vergezeld van extra vergelijkingen voor de onbekenden, die overigens geen beperking van de algemeenheid inhouden. Om te kunnen vaststellen of dit element eventueel geschikt is voor ons doel, hoeven we alleen maar een optelsom te maken van vergelijkingen en onbekenden. Ga uit van een raster van $I$ bij $J$ van zulke vierhoeken. We maken een staatje op: $$ \begin{array}{ccc} \mbox{aantal vergelijkingen bulk} &=& 2 I.J \\ \mbox{aantal randvoorwaarden} &=& 2I + 2J \\ \mbox{de extra vergelijkingen} &=& 2 I.J \\ \mbox{Bij elkaar opgeteld} &=& 4 I.J + 2I + 2J \\ \mbox{Totale aantal onbekenden} &=& 2 \{ (I+1)J + I(J+1) \} \end{array} $$ Het aantal vergelijkingen is dus exakt gelijk aan het aantal onbekenden. Onze kleinste kwadraten eindige elementen methode is equivalent met een stelsel van eindige differentie vergelijkingen. Het kleinste kwadraten minimum is om te beginnen al gelijk aan nul.

Doet deze kleinste kwadraten methode het nu ook werkelijk? Ik laat het over aan de lezer om dit eigenhandig te kontroleren. Als appendix opgenomen is een serie BASIC programma's "BASIC HdB". Men dient de programmatuur [FTP] in de aangegeven volgorde te laten werken.

Heel wat details zijn weliswaar geprogrammeerd, maar zullen hier niet worden behandeld. (Dit komt helaas wel vaker voor bij het ontwikkelen van software.) Ik noem het gebruik van "repeteerbare elementen": het valt te bewijzen dat de gebezigde element-matrices onafhankelijk zijn van rotaties en schaalfaktoren, waardoor alle matrices identiek zijn. Het plotten van stroomlijnen is in 2-D hetzelfde als het plotten van contouren van de stroomfunktie. De laatste wordt eveneens met een kleinste kwadraten methode uitgerekend. Het eerste, grafische representatie van contouren, is een verhaal apart, over hoe men in de praktijk bijvoorbeeld te maken krijgt met Klassifikatie (Equivalentie Klassen: diskrete ondergrond).