Procedures voor geïtereerde afbeeldingen in 2D

<< Utilities`FilterOptions`

iteraties[f_, x0_, aantalinitialiseren_, aantalhouden_] := 
      Drop[NestList[f, x0, aantalinitialiseren + aantalhouden], 
        aantalinitialiseren]

\(\(iteraties::usage = "iteraties[f,x0,aantalinitialiseren,aantalhouden] past f eerst aantalinitialiseren keer op x0 toe. Daarna wordt f nog aantalhouden keer toegepast: de resultaten hiervan zijn de output van deze functie.";\)\)

f[α_][{t_, v_}] = {t + v, α*v};

iteraties[f[0.9], {0, 1.45}, 0, 50]

faseplaatje[f_, x0_, aantalinitialiseren_, aantaltonen_, opts___] := 
        ListPlot[
Table[
            iteraties[f, x0, aantalinitialiseren, 
              aantaltonen]], 
AxesLabel -> {t, v}, 
          PlotStyle -> AbsolutePointSize[5], 
          opts
];

\(\(faseplaatje::usage = \*""faseplaatje[f,x0,aantalinitaliseren,aantaltonen] geeft aantaltonen iteranden van f weer, waarbij f een functie van (ℝ^2) naar (ℝ^2) moet zijn."";\)\)

faseplaatje[f[0.9], {0, 1.45}, 0, 50];

Nu naar de aangedreven stuiterbal:

g[α_, γ_][{t_, v_}] = {t + 
            v, α v - γ Cos[t + v]};

Even kijken of hij voor γ = 0 hetzelfde doet als de niet aangedreven functie f:

faseplaatje[g[0.9, 0], {0, 1.45}, 0, 50];

Deze g gaat van (ℝ^2) naar (ℝ^2), we gaan hem nu herdefiniëren om van de cirkel naar de cirkel te gaan (door te tijd modulo 2  π te nemen):

g[α_, γ_][{t_, v_}] = {Mod[t + v, 
            2  π], α v - γ Cos[t + v]};

faseplaatje[g[0.9, 0], {0, 1.45}, 0, 50];

Als we de tafel nu snelheid geven, bijvoorbeeld γ = 1, dan ziet het faseplaatje er iets interessanter uit:

faseplaatje[g[0.9,  1],  {4,  6},  100, 150];

De snelheid kan soms negatief worden. Dit zou je natuurlijk kunnen oplossen door te definiëren

g2[α_, γ_][{t_, v_}] = {Mod[t + v, 2  π], 
          Max[α v - γ Cos[t + v], 0]};

Dit doen we echter niet. De negatieve snelheden zijn een bijkomstigheid van een modelleerfout die we hebben gemaakt. We doen voorlopig gewoon of onze neus bloedt.

We bekijken nu wat g met verschillende beginwaarden doet (met vaste α en γ). Hiertoe kijken we eens wat hij doet met een stel punten. Eerst een functie om zo'n stel punten te tekenen.

wolkplot[f_, input_, 
        opts___] := (
Block[{$DisplayFunction = 
              Identity}, 
output = 
            Map[f, input]; 
inputplot = 
            ListPlot[
              input, 
PlotStyle -> (inputstyle /. 
                    Append[{opts}, 
                      inputstyle -> 
                        Red])]; 
outputplot = 
            ListPlot[
              output, 
PlotStyle -> (outputstyle /. 
                    Append[{opts}, 
                      outputstyle -> 
                        Blue])];
]; 
Show[inputplot, outputplot, 
          AxesLabel -> {t, v}, 
FilterOptions[
            ListPlot, opts], 
          PlotRange -> {{0, 2  π}, Automatic}])

\(\(wolkplot::usage = \*""wolkplot[f,input] geeft voor een lijst punten in (ℝ^2) een plaatje van die punten en de functie f toegepast op die punten. Met de opties inputstyle en outputstyle kan worden opgegeven hoe die punten er uit moeten zien."";\)\)

wolk = Table[{t, 20}, {t, 4, 5, 0.02}]

wolkplot[g[0.9, 16], wolk, inputstyle -> Red, 
        outputstyle -> Blue];

We gaan dit wat uitgebreider bekijken, door te kijken wat er met een bepaalde vierhoek gebeurt, en wel alleen naar de randen (je mag natuurlijk zelf de hoekpunten wat aanpassen).

linksboven = {0, 25};
rechtsboven = {2  π, 17.9};
linksonder = {0, 19.7};
rechtsonder = {2  π, 13.5};
stapjes = 300;
links = 
        Table[linksonder + (i/stapjes)*(linksboven - linksonder), {i, 0, 
            stapjes}];
rechts = 
        Table[rechtsonder + (i/stapjes)*(rechtsboven - rechtsonder), {i, 
            0, stapjes}];
onder = 
        Table[linksonder + (i/stapjes)*(rechtsonder - linksonder), {i, 0, 
            stapjes}];
boven = 
        Table[linksboven + (i/stapjes)*(rechtsboven - linksboven), {i, 0, 
            stapjes}];
Clear[stapjes];

h = g[0.9, 16];
Block[{$DisplayFunction = Identity}, 
linksplot = 
          wolkplot[h, links, outputstyle -> Green, 
            inputstyle -> Green]; 
rechtsplot = 
          wolkplot[h, rechts, outputstyle -> Black, 
            inputstyle -> Black]; 
onderplot = 
          wolkplot[h, onder, outputstyle -> Blue, 
            inputstyle -> Blue]; 
bovenplot = 
          wolkplot[h, boven, outputstyle -> Magenta, 
            inputstyle -> 
              Magenta];
];
Show[linksplot, rechtsplot, onderplot, bovenplot];

Nog enige losse opmerkingen t.b.v. het verslag (allemaal losse tips die niet al te serieus genomen moeten worden):
 - Wiskundigen hebben de neiging altijd lelijke kleuren te kiezen in hun plaatjes. Probeer eens of je hier een uitzondering op kunt zijn.
 - Als je een plaatje uit Mathematica wilt exporteren voor in je LATEX-verslag, kun je dit het beste doen via een EPS- of PDF-bestand. Exporteer eerst het plaatje als EPS (rechtermuisknop op het plaatje, menu Edit, Save Selection As, EPS) en converteer het dan via de terminal eventueel naar PDF met het commando eps2pdf. Grafieken in JPG zien er erg knullig uit (JPG is een compressietechniek die op gladde plaatjes zoals foto's goed werkt, maar zeker niet voor lijntekeningen).
 - Als je conclusies in je verslag schrijft, laat die dan enigszins zinvol zijn.
 - Probeer niet de opdracht letterlijk over te schrijven, het staat veel leuker als je er zelf een beetje een verhaal van maakt.

Bifurcaties

\(bifurcatiefilmpje[γbegin_, γeind_, γstapjes_, 
        x0_List, 
        opts___] := 
Module[{initialiserenbegin = 400, 
          initialiseren = 50, plotten = 100, 
          lijst, γ}, 
lijst = 
          iteraties[g[0.9, γbegin], x0, initialiserenbegin, 
            plotten]; 
For[γ = γbegin, γ ≤ γeind, γ = γ + \(γeind - γbegin\)\/γstapjes, 
x00 = 
            Last[lijst]; 
lijst = 
            iteraties[g[0.9, γ], x00, initialiseren, 
              plotten]; 
ListPlot[lijst, opts, 
            PlotLabel -> \(("γ = " <> 
                  ToString[PaddedForm[γ, 3]] <> "
x = " <> 
                  ToString[PaddedForm[Last[lijst], 5]])\), 
            PlotStyle -> AbsolutePointSize[5], 
            AxesLabel -> {"t", "v"}]
]
]\)

\(\(bifurcatiefilmpje::usage = "bifurcatiefilmpje[γbegin,γeind,γstapjes,x0] geeft een serie faseplaatjes van de afbeelding g[0.9,γ] met γ in een bepaalde range. Dubbelklik op het eerste plaatje om het als filmpje te bekijken";\)\)

\(bifurcatieplaatje[γbegin_, γeind_, γstapjes_, 
        x0_List, 
        opts___] := 
Module[{initialiserenbegin = 400, 
          initialiseren = 50, plotten = 100, lijst, γ, 
          output = {}}, 
lijst = 
          iteraties[g[0.9, γbegin], x0, initialiserenbegin, 
            plotten]; 
For[γ = γbegin, γ ≤ γeind, γ = γ + \(γeind - γbegin\)\/γstapjes, 
x00 = 
            Last[lijst]; 
lijst = 
            iteraties[g[0.9, γ], x00, initialiseren, 
              plotten]; 
output = 
            Join[output, \((\({γ, #}
 &\) /@ \((First /@ lijst)\))\)]]; 
ListPlot[output, opts, 
          AxesLabel -> {"γ", "t"}]
]\)

\(\(bifurcatieplaatje::usage = "bifurcatieplaatje[γbegin,γeind,γstapjes,x0] geeft een bifurcatieplaatje voor de functie g[0.9,γ] met γ in een bepaalde range. Let op afhankelijkheid van beginvoorwaarden. Alleen de t-waarde wordt geplot (verander in de code de First/@lijst in Last/@lijst om de v-waarde te plotten).";\)\)

bifurcatiefilmpje[3.3, 5, 30, {0, 4.1}, 
        PlotRange -> {{0, 2  π}, {-3, 9}}];

bifurcatieplaatje[3.3, 5, 30, {0, 4.1}];