Logistische groei

We gaan eens kijken naar de groei van een populatie die zich gedraagt volgens het Verhulst model. Zo'n populatie groeit volgens een bepaalde formule, de populatie op tijdstip i + 1 volgt uit de populatie op een tijdstip i, en wel volgens

(Subscript[x, i + 1]) = a (Subscript[x, i]) (1 - (Subscript[x, i])).

We noemen deze functie (Subscript[f, a]) (het subscript a geeft aan dat er een parameter a in zit), zodat we hebben

(Subscript[x, i + 1]) = (Subscript[f, a])((Subscript[x, i])).

We definiëren voor Mathematica hoe de functie (Subscript[f, a]) er uitziet:

(Subscript[f, a])_[x_] = a x (1 - x);

(Om een subscript te maken, gebruik je de toetscombinatie Control-minteken.)

Even kijken of hij het doet:

(Subscript[f, a])[x]

En als we wat getallen invullen:

(Subscript[f, 3.1])[0.5]

En inderdaad,

3.1 * 0.5 (1 - 0.5)

Als we nu een beginpopulatie nemen van (geschaalde) grootte (Subscript[x, 0]) = 0.5, dan is het interessant om te kijken wat er op langere termijn gebeurt:

(Subscript[f, 3.1])[0.5]

Hier is het makkelijk het Mathematica-symbool % te gebruiken, wat staat voor 'het laatst geproduceerde antwoord'.

(Subscript[f, 3.1])[%]

En nog een generatie verder:

(Subscript[f, 3.1])[%]

(Subscript[f, 3.1])[%]

(Subscript[f, 3.1])[%]

(Subscript[f, 3.1])[%]

(Subscript[f, 3.1])[%]

Dit schiet natuurlijk niet op, we moeten veel te veel met de hand typen. (Is er trouwens iets bijzonders met bovenstaande antwoorden?) We kunnen het Mathematica-commando NestList gebruiken. Dit maakt een lijst met alle tussenantwoorden. Als we nog geen getallen invullen zie je wat er gebeurt:

NestList[f, (Subscript[x, 0]), 4]

Vullen we dezelfde waarden als net in, dan zien we ook dezelfde resultaten:

NestList[(Subscript[f, 3.1]), 0.5, 7]

Nu is het ook makkelijk om wat langer door te gaan (we hoeven niet meer voor elke generatie een regel te typen):

NestList[(Subscript[f, 3.1]), 0.5, 50]

Je ziet dat er een patroon ontstaat. In dit patroon zijn we geïnteresseerd, het maakt niet zoveel uit wat er in de eerste zoveel iteraties gebeurt. Die gooien we daarom weg, met het commando Drop.

Drop[%, 40]

Vegen we nu wat van de vorige commando's bij elkaar (zodat we % niet meer nodig hebben), dan is het totale commando om te kijken wat er na een bepaald aantal generaties gebeurt:

Drop[
NestList[(Subscript[f, 3.1]), 0.5, 
        50], 
40]

Nu hebben we gekeken wat er voor het geval a = 3.1 gebeurt. Kijk zelf eens wat er bij a = 3.5 gebeurt. Ontstaat er weer een 2-baan?

We breiden het wat uit, en gaan kijken wat er onstaat bij algemene a. Daarvoor gebruiken we het commando Table, dat een tabel maakt van de oplossing voor elke a in de lijst. Het gedeelte {a, #} &/@ zorgt ervoor dat voor elke oplossing (Subscript[x, i]) ook nog staat bij welke a hij hoorde.

populaties = 
Table[
Drop[
{a, #} & /@ 
            NestList[(Subscript[f, a]), 0.5, 
              50], 
40], 
{a, 2, 
          4,  .1}
]

We kunnen dit opvatten als een lange lijst van (x, y) coördinaten, en we kunnen ze ook in een plaatje zetten. Hiervoor moet eerst de lijst wat worden opgeschoond:

populaties = Flatten[populaties, 1];

ListPlot[populaties];

We kunnen dit allemaal samenvatten in één functie, die het hele plaatje voor je maakt:

bifurcatie[aantalprinten_, aantalinitialiseren_, abegin_, aeind_, 
          astapgrootte_] := 
ListPlot[
Flatten[
Table[
Drop[
{a, #} & /@ 
                  NestList[(Subscript[f, a]),  .5, 
                    aantalprinten + 
                      aantalinitialiseren], 
aantalinitialiseren], 
{a, abegin, 
                aeind, astapgrootte}
], 
1
]
];

\(\(bifurcatie::usage = "bifurcatie[aantalprinten,aantalinitialiseren,abegin,aeind,astapgrootte] maakt een bifurcatiediagram van het Verhulst-model, waarbij eerst aantalinitialiseren generaties worden berekend en weggegooid en daarna aantalprinten generaties worden opgenomen in het diagram. De parameter a loopt van abegin tot aeind in stapjes ter grootte astapgrootte.";\)\)

bifurcatie[10, 20, 2, 4, 0.01];