De Verende Slinger

Visualisatie van mate van voorspelbaarheid van het onderzochte systeem in het bereik x = [0, 2.5] y = [-0.5, 1]

Versie 2

Mathijs de Bruin 0420891

(bij eerder versie: Els Hoogeveen 0373583)

Datum: 27-9-2005
Docent: Wolter Kaper 
Vak: Mathematica 
Bijlage: batch-bestanden en pdf-files met script voor uitvoering onder unix

Samenvatting

Bij dit project heb ik een numeriek fysisch model onderzocht bestaande uit een veer die vrij kan draaien om één as met een gewicht eraan onder invloed van een zwaartekracht die naar beneden trekt. Hierbij werden als randvoorwaarden de startpositie, het gewicht (voorgesteld als puntmassa) aan de veer genomen en de initieële snelheid van 0.

Vervolgens keken we naar het verschil in uitkomst als gevolg van veranderingen in de randvoorwaarden, zowel met behulp van animaties van specifieke randvorwaarden als met afbeeldingen die de mate van voorspelbaarheid visualiseren. Hierbij bleken er in deze afbeeldingen interssante patronen te zitten welke zich naarmate de tijd vordert ontwikkelen.

Inleiding

Een springveer is opgehangen aan een horizontaal rond balkje, zodat vrije rotatie rond het balkje mogelijk is. Aan het andere uiteinde van de springveer hangt een massa (zie figuur). We veronderstellen dat de veer zowel bij uitrekken als bij indrukken recht blijft.
We bekijken beginsituaties waar de springveer zich in een vlak loodrecht op het balkje bevindt. Het systeem heeft dan 2 vrijheidsgraden: de hoek θ ten opzichte van de richting van de zwaartekracht, en de lengte L', die kan verschillen van de evenwichtslengte L van de ontspannen veer.

Net als we deden voor de starre slinger, kunnen we Newtons wet  F = m a toepassen op deze situatie, want we kennen de zwaartekracht ((Subscript[F, z]) = m g) en we kennen de veerkracht (Wet van Hooke: (Subscript[F, v]) = k u. Hier is u de uitrekking van de veer). Dat leidt nu tot twee differentiaalvergelijkingen, met twee onbekende functies van t erin (naar keuze θ(t) en L'(t) of de gewone coordinaten x(t) en y(t)). 
Bij sommige begincondities doet zich bij het oplossen een probleem voor: de gevonden oplossing blijkt erg gevoelig voor kleine verschillen in begincondities, zodat een willekeurig kleine onnauwkeurigheid in de beginvoorwaarde uiteindelijk tot een heel andere baan leidt. Omdat begincondities (in de fysica) altijd maar in een beperkt aantal decimalen bekend zijn, betekent dit, dat deze slinger onder bepaalde condities principieel onvoorspelbaar is, op den duur. (Het weer is om dezelfde redenen onvoorspelbaar).
We onderzoeken wanneer deze onvoorspelbaarheid zich voordoet.

[Bron: Opgave 'Project: de verende slinger']

Werkwijze

Allereerst zijn er een aantal schetsjes getekend van de werkende krachten. In het geval van dit model spelen alleen de veerkracht (Subscript[F, v]) en de zwaartekracht (Subscript[F, z]) een rol. Verder heb ik ter naslag (bij de vorige versie) 'Differentiaalvergelijkingen: Slinger over de kop', behorende tot het Mathematica-practicum, gebruikt om een beeld te scheppen van een mogelijke aanpak.

Ik kwam na het nodige overleg tot de conclusie dat de dubbele afgeleide van respectievelijk de x en de y positie van het gewicht (de versnelling dus) te berekenen is als alle krachten op het gewicht bekend zijn. Dit gegeven als handrijking genomen te hebben hoefde ik slechts de formule voor de netto kracht te schrijven in zowel de x als y richting en dit om te schrijven naar de versnelling.

Dit alles is mogelijk dankzij het verband van kracht en versnelling voorvloeiend uit \(TraditionalForm\`F\&_ = m  a\&_\). Verder is er een veerkracht gegeven door (Subscript[F, v]) = uk waarvoor geldt u = (Sqrt[(x^2) + (y^2)]). Tezamen met (Subscript[F, z]) = mgkunnen we nu alle krachten op de massa bepalen.

Dit is simpelweg in te vullen in een differentiaalvergelijking die wij Mathematica numeriek kunnen laten oplossen. Vervolgens heb ik een functie geschreven welke deze differentiaalvergelijking eenvoudig oplost als functie van x en y, danwel evalueert als functie van t, x en y. Met andere woorden: als alleen de x en y als randvoorwaarden gegeven zijn komen hier twee functies uit terug, als ook de t gegeven is komt hier de positie van de slinger op de gewenste coördinaten uit.

Om het geheel wat te versimpelen hebben we de gravitatieconstante g op 1 gesteld, net als de veerconstante k, de massa m en de (oorspronkelijke) lengte van de veer L.

Met de nu verkregen functie kunnen we nu met behulp van ParametricPlot eenvoudig animaties maken van de hele weg die de puntmassa doorloopt. Hierna ben ik een behoorlijke stap verder gegaan door het berekenen van de mate van invloed van de startcoördinaten op de eindcoördinaat na een bepaalde tijd. Dit hebben we gedaan door een schatting te maken van de afgeleide naar x respectievelijk y en deze bij elkaar op te tellen en hiervan de lengte te nemen.

Om deze schatting te kunnne maken heb ik de volgende formule voor de afgeleide gebruikt, waarbij voor h een waarde van 0.0001 is gebruikt, welke proefondervindelijk afdoende klein lijkt te zijn. Van deze formule hebben we de lengte genomen en deze met behulp van DensityPlot tot een grafiek verwerkt.

δf[t, x + y]/δx + δf[t, x, y]/δy ≈ ((f[t, x + h, y] - f[t, x - h, y])/(2  h)) + ((f[t, x, y + h] - f[t, x, y + h])/(2  h))

Op deze manier kunnen we een eenvoudig interpreteerbare afbeelding laten zien die de mate van voorspelbaarheid ('chaos') in dit systeem aanduidt. Een dergelijk plaatje is al te zien aan het begin van dit project. Hierbij stellen de witte gebieden punten voor waar de mate van chaos groot is en hier zal dus een kleine wijziging in de randvoorwaarden een grote uitwerking op de uiteindelijke positie hebben.

Vervolgens heb ik de plaatjes eens iets beter bekeken.

Onder andere heb ik gekeken of de mate van detail naarmate er een kleiner gebied wordt bekeken (bij een gelijkblijvende resolutie) afneemt, zoals je dit bij een normale afbeelding zou verwachten. Daarnaast heb ik gekeken naar de manier waarop het systeem zich ontwikkeld als functie van de tijd door een plaatje op een beperkte resolutie na verschillende tijden te berekenen.

Omdat het berekenen van de plaatjes nogal veel tijd kost heb ik ervoor gekozen om deze uit te voeren op andere computers met behulp van batch-jobs. Zodoende heb ik deze plaatjes simpelweg geplakt binnen dit notebook. Als het goed is zouden echter de gebruikte batch-jobs met het resultaat hiervan in PDF-formaat met dit project meegestuurd moeten worden.

Berekeningen en Resultaten

De Schrijver (functies)

Off[General::spell]

<< Graphics`Animation`

m = 1;
L = 1;
k = 1;
g = 1;

(Subscript[F, x]) = ((k X[t] (L - (Sqrt[X[t]^2 + Y[t]^2])))/(Sqrt[X[t]^2 + Y[t]^2]));
(Subscript[F, y]) = ((k Y[t] (L - (Sqrt[X[t]^2 + Y[t]^2])))/(Sqrt[X[t]^2 + Y[t]^2])) - g m;

dv1 = X''[t] == ((Subscript[F, x])/m);
dv2 = Y''[t] == ((Subscript[F, y])/m);

Commentaar: omdat t in de animaties hier maar tot 50 loopt hoeft de differentiaalvergelijking ook maar tot hier berekend te worden; dit scheelt tijd. Het spreekt dus vanzelf dat bij een animatie die verder loopt dan tot 50 ook deze tijd voor de differentialvergelijking verandert moet worden, en dat deze functie opnieuw geëvalueerd moet worden.

\(func[time_, x_, 
        y_] := \((
data\  = \ {\ 
            dv1, 
\ \ \ \ \ \ \ \ \ \ \ \ \ dv2, 
\ \ \ \ \ \ \ \ \ \ \ \ \ X[0] == 
              x, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \(X'\)[
                0] == 
              0, 
\ \ \ \ \ \ \ \ \ \ \ \ \ Y[0] == 
              y, 
\ \ \ \ \ \ \ \ \ \ \ \ \ \(Y'\)[
                0] == 0}; 

        result = {X[t], 
              Y[t]}\  /. \ \(NDSolve[data, {X[t], Y[t]}, {t, 0, 50}, \ 
                MaxSteps -> 
                  Infinity]\)[\([1]\)]; 
Return[
          result\  /. \ t -> time]
)\)\)

h = 0.0001;

diffuncx[t_, x_, y_] := 
      Norm[func[t, x + h, y] - func[t, x - h, y]]/(2 h)
diffuncy[t_, x_, y_] := 
      Norm[func[t, x, y + h] - func[t, x, y - h]]/(2 h)
diffunc[t_, x_, y_] := 
      diffuncx[t, x, y] + diffuncy[t, x, y]

\(animate[path_, start_, \ stop_, \ 
        step_] := \((
images\  = \ {}; 
Do[
images\  = \ 
            Append[images, 
              ParametricPlot[
path, \ {t, start, 
                  current + step}, \ 
                PlotRange -> {{\(-10\), 10}, {\(-10\), 10}}, 
                DisplayFunction -> 
                  Identity]], 
{current, start, 
            stop - step, \ step}]; 
ShowAnimation[
          images]; 
SelectionMove[EvaluationNotebook[], 
            All, GeneratedCell]
\(FrontEndTokenExecute["OpenCloseGroup"];\)
Pause[
            0.01]
FrontEndTokenExecute["SelectionAnimate"];
)\)\)

On[General::spell]

Het Boek (animatie)

Noot: zie het commentaar bij de definitie van de functie func.

anim  =  func[t, 3, 3];
animate[anim, 0, 50, 0.5];

De Bibliotheek (chaos in kaart)

Hier heb ik geprobeert om een duidelijke weergave te maken van de mate van chaos in het systeem. De eerste berekening heb ik, meer bij wijze van voorbeeld, op deze computer berekend. De andere afbeeldingen zijn in batches berekend vanwege de lange rekentijden die hiermee gemoeid zijn. Hierbij heb ik zo goed mogelijk geprobeert aan te geven wat er precies wordt weergegeven en wordt de bestandsnaam van het bestand waar de batch-commando's in staan geven.

De batch-commando's zijn uitgevoerd op de Linux-computers voor studenten op de faculteit, met behulp van het scriptje do-it.sh dat met de batch-commando's zelf bij dit project zou moeten zitten.

Helaas zijn bij het overbrengen van de afbeeldingen van de PDF naar Mathematica de assen verloren gegaan. Bij behoefte hieraan zijn de oorspronkeljike PDF-bestanden beschikbaar.

Afbeelding 120x120 (-8 < x < 8, -8 < y 8, t=100)

Dit is de enige afbeelding die ik bij wijze van voorbeeld wel binnen dit notebook heb staan. Zowel de x als de y lopen tussen -8 en 8 en er zijn per as 120 punten genomen. Hierbij is de resolutie dus 120x120.
Om de hieronder staande berekening over te doen moet de bovenwaarde van t bij het oplossen van de differentialavergelijking bij de definitie van func verandert worden.

Timing[DensityPlot[
          diffunc[100, x, y], {x, -8, 8}, {y, -8, 8},  
          Mesh  ->  False, 
          PlotPoints -> 
120];]

Afbeelding 400x400 (-10 < x < 10, -10 < y < 10, t=100)

Zoals te zien is de mate van detail in deze afbeelding alles behalve afgenomen. Dit geldt ook voor vergrotingen hiervan, zoals in volgende afbeeldingen te zien is.

Bestandsnaam batch: math-400.m
Bestandsnaam PDF: math-400-ow150.pdf

Afbeelding 600x600 (0 < x < 2.5, -0.6 < y < 1, t=100)

Batch-bestand: frac-600.m
PDF-bestand: frac-600-ow143.pdf

Afbeelding 600x600 (1 < x < 2, 0 < y < 0.8, t=100)

Batch-bestand: frac2-600.m
PDF-bestand: frac2-600-ow162.pdf

Afbeelding 200x200 (-10 < x < 10, -10 < y < 10, t={100,300,400,500})

Hier is in een serie plaatjes duidelijk te zien dat er zich in de mate van voorspelbaarheid nog interessante ontwikkelingen voordoen naarmate de tijd voordert. Overigens is het renderen van de afbeelding bij t=200 om misterieuze reden mislukt en vanwege de rekenintensiviteit niet opnieuw berekend.

Discussie

Dankzij een goede inval (en wat extra tijd thuis, wegens ziekte) ben ik er bij deze versie van het project (wel) in geslaagd om een grondiger onderzoek te doen naar de voorwaarden voor chaotische effecten van het hier onderzochte model, met een niet te verwaarlozen resultaat.

Met behulp van de ontworpen 'voorspelbaarheidsfunctie' heb ik een redelijk betrouwbaar lijkende weergave kunnen maken van de mate van voorspelbaarheid van het fysische systeem. Naar aanleiding van de berekende visualisaties hiervan is bij mij sterk het vermoeden onstaan dat het hier een figuur betreft met een fractale dimensie. Dit met name omdat de hoeveelheid weergegeven informatie bij het 'inzoomen' niet af lijkt te nemen.
Als gevolg hiervan stel ik mij dan weer de vraag wat de fractale dimensie hiervan is en aansluitend hoe ik dit zou kunnen berekenen.

Ook ben ik benieuwd wat voor patronen er zouden (kunnen) ontstaan als ik de startpositie-afhankelijkheid afzonderlijk zou bekijken naar de x- danwel de y-as van het systeem.

Naast dat de met deze manier van onderzoeken gemaakte plaatjes van een onmiskenbare esthetische waarde zijn heeft het kleine onderzoekje verder vooral veel vragen opgeroepen omtrend de tijdsafhankelijkheid en de hoeveelheid informatie die langs deze weg verkrijgbaar is. Zijn er patronen volgens welke het hier onderzochte systeem zich in de tijd ontwikkeld, en in hoeverre vertonen deze fractale eigenschappen? Wat voor patronen ontstaan er als er voor een nóg kleiner stukje van nóg 'dichterbij' een plaatje wordt gemaakt?

Verder zijn er nog een aantal andere vragen die bij mij hunkeren naar beantwoording, bijvoorbeeld in hoeverre de visualisaties afhangen van het gebruikte (fysische) model. Wat zou er bijvoorbeeld gebeuren als ik voor de gravitatieconstante niet 1 maar het 'aardse' 9,81 m/(s^2) had ingevuld?

Helaas heb ik vooral uit (reken)tijdsgebrek de hierboven geformuleerde vragen niet kunnen beantwoorden.

Tenslotte heb ik de hier gebruikte werkwijze wel voor een aantal willekeurige waardes van de 'voorspelbaarheidsfunctie' handmatig getest op de mate van correspondentie met de waargenomen afhankelijkheid van startpositie door middel van vergelijking met de animaties (met een overigens positief resultaat) maar zou dit natuurlijk eigenlijk structureel moeten gebeuren om mezelf ervan te verzekeren dat wat er feitelijk wordt weergegeven niet gewoon (deels) een weergave is van afrondingsfouten van de numerieke differentiaalvergelijkings-oplosser van Mathematica. Bij deze handmatige tests (welke vanwege hun willekeur, en het feit dat ik ze niet heb opgeslagen, niet bij dit verslagje vermeld zijn) viel het mij overigens wel op dat er heel grote verschillen zitten in de waarde die de functie teruggeven, vaak in de orde van (10^4).