Preamble (wordt automatisch uitgevoerd):

Iterations[f_,expr_,{m_Integer,n_Integer}]:=NestList[f,Nest[f,N[expr],m],n-m]/;0≤m≤n
Iterations[f_,expr_,n_Integer]:=NestList[f,N[expr],n]
Iterations::usage="Iterations[f,expr,n] gives a list of the result of applying f to expr 0 through n times, applying N after each iteration to reduce computation time when exact results are possible.
Iterations[f,expr,{m,n}] gives a list of the result of applying f to m through n times.";

Unprotect[BifurcationPlot];
Options[BifurcationPlot]=Join[Options[ListPlot],{ParameterSteps->250}];
SyntaxInformation[BifurcationPlot]={"ArgumentsPattern"->{_,_,{_,_},{_,_,_,_.},OptionsPattern[BifurcationPlot]},"LocalVariables"->{"Table",{4,4}}};
BifurcationPlot[f_,x0_,steps:{_,_},{par_,parmin_,parmax_,parstep_:Null},options:OptionsPattern[BifurcationPlot]]:=
ListPlot[Flatten[Table[{par,#}&/@Drop[NestList[f,N[x0],Max@steps],Switch[steps,_List,Min@steps,_,0]],{par,parmin,parmax,Replace[parstep,Null->(parmax-parmin)/OptionValue[ParameterSteps]]}],1],Sequence@@FilterRules[{options},Options[ListPlot]]]
SetAttributes[BifurcationPlot,{HoldAll,Protected}];
BifurcationPlot::usage="BifurcationPlot[f,(Subscript[x, 0]),{m,n},{a,(Subscript[a, min]),(Subscript[a, max]),(Subscript[a, step])}] produces a bifurcation plot for the recursion (Subscript[x, n+1])=f((Subscript[x, n])) .
The argument (Subscript[a, step]) is optional, it is ((Subscript[a, max])-(Subscript[a, min]))/N by default, where N is the value of the option ParameterSteps (default value: 250).
BifurcationPlot admits all options for ListPlot; these are passed on.";

Unprotect[Recursion2D];
Options[Recursion2D]=Join[Options[ListPlot],Options[Manipulate],{Locator->True,PointStyle->{Red,PointSize[Medium]},InitialPointName->"(Subscript[x, 0])"}];
SyntaxInformation[Recursion2D]={"ArgumentsPattern"->{_,_.,{{_,_},{_,_}},_,___,OptionsPattern[Recursion2D]}};
Recursion2D[f_,x0:{_,_},ranges:{{_,_},{_,_}},fixedpoints_:{},steps:(_Integer|{_Integer,_Integer}),pars___,options:OptionsPattern[Recursion2D]]:=Module[{x,
listplotoptions=Sequence@@FilterRules[{options},Options[ListPlot]],
manipulateoptions=Sequence@@FilterRules[{options},Options[Manipulate]],
pointstyle=Flatten[{OptionValue[PointStyle]}],
initialpointname=OptionValue[InitialPointName]},
Manipulate[Show[
ListPlot[Quiet[If[MatchQ[f[N[x]],{_Real,_Real}],Drop[NestList[f,N[x],Max@steps],Switch[steps,_List,Min@steps,_,0]],Message[Recursion2D::nan,f,x];{x}],{General::ovfl,General::unfl,Power::infy}],
listplotoptions//Evaluate,PlotRange->ranges],
Graphics[Append[pointstyle,Point[{fixedpoints}]]]],
pars,If[OptionValue[Locator]//Evaluate,Sequence@@{{{x,x0,initialpointname}},{{x,x0},Locator}},{{x,x0,initialpointname}}]//Evaluate,manipulateoptions//Evaluate]
]
Recursion2D[f_,Automatic|Null,ranges:{{_,_},{_,_}},fixedpoints_:{},steps:(_Integer|{_Integer,_Integer}),pars___,options:OptionsPattern[Recursion2D]]:=Recursion2D[f,Evaluate[Mean/@ranges],ranges,fixedpoints,steps,pars,options]
SetAttributes[Recursion2D,{HoldAll,Protected}];
Recursion2D::usage="Recursion2D[f,{(Subscript[x, 0]),(Subscript[y, 0])},{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},iterations,pars] creates a listplot within a Manipulate for the points obtained by repeatedly applying f to {(Subscript[x, 0]),(Subscript[y, 0])}.
Recursion2D[f,{(Subscript[x, 0]),(Subscript[y, 0])},{{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}},fixedpoints,iterations,pars] does the same and also draws the points specified by fixedpoints, which should be of the form {{(Subscript[X, 1]),(Subscript[Y, 1])},{(Subscript[X, 2]),(Subscript[Y, 2])}…}.
The range {{(Subscript[x, min]),(Subscript[x, max])},{(Subscript[y, min]),(Subscript[y, max])}} is used as the plot range for ListPlot unless something else (such as PlotRange->Automatic) is specified using the Plotrange option. Providing the starting point {(Subscript[x, 0]),(Subscript[y, 0])} is optional, if it is left out {Mean[(Subscript[x, min]),(Subscript[x, max])],Mean[(Subscript[y, min]),(Subscript[y, max])]} is used.
Controls are created for the paramaters specified by pars, as well as the starting point {(Subscript[x, 0]),(Subscript[y, 0])}, which may be left out (in which case it starts at the centre of the specified range). The syntax for pars is the same as for Manipulate.
Which iterations are shown is determined by iterations, which should either be an integer n or of the form {m,n}. The iterations m through n are shown (m=0 if it is left unspecified).
Recursion2D has three options that are unique to it: InitialPointName (default value: \"(Subscript[x, 0])\"), PointStyle (default value: {Red,PointSize[Medium]}) and Locator (default value: True). Other options are passed on to ListPlot or Manipulate where applicable.";
Recursion2D::nan="`1`[`2`] is not a pair of real numbers.";

Unprotect[IterationDiagram];
Options[IterationDiagram]=Join[Options[Plot],{LineStyle->Blue,PointStyle->{Blue,PointSize[Medium]}}];
SyntaxInformation[IterationDiagram]={"ArgumentsPattern"->{_,{_.,_,_},_,_,OptionsPattern[IterationDiagram]},"LocalVariables"->{"Plot",{2,2}}};
IterationDiagram[f_,{xmin_,xmax_},x0_,n_,options:OptionsPattern[IterationDiagram]]:=Module[{data=Quiet[Replace[Cases[Iterations[Composition[N,f],x0,n],_?MachineNumberQ],{a_:>Sequence[a,a]},{1}],{General::ovfl,General::unfl,Power::infy}],
pointstyle=Flatten[{OptionValue[PointStyle]}],linestyle=Flatten[{OptionValue[LineStyle]}]},
data=Drop[Transpose[{data,RotateLeft[data,1]}],-1];
Show[Plot[{f[x],x},{x,0,1},Sequence@@FilterRules[{options},Options[Plot]]//Evaluate,PlotStyle->{Automatic,Black}],Graphics[{Append[linestyle,Line[data]],Append[pointstyle,Point[{x0,x0}]]}]]
]
IterationDiagram[f_,{x_,xmin_,xmax_},x0_,n_,options:OptionsPattern[IterationDiagram]]:=IterationDiagram[Function[f/.x->#//Evaluate]//Evaluate,{xmin,xmax},x0,n,options]
SetAttributes[IterationDiagram,{HoldAll,Protected}]
IterationDiagram::usage="IterationDiagram[f,{x,(Subscript[x, min]),(Subscript[x, max])},(Subscript[x, 0]),n] shows (Subscript[x, 0]),(Subscript[x, 1]),…,(Subscript[x, n]) defined by (Subscript[x, n+1])=f((Subscript[x, n])) in the plot of f, connected by lines.
IterationDiagram[g,{(Subscript[x, min]),(Subscript[x, max])},(Subscript[x, 0]),n] gives the same output if g is a function for which g[x]=f.
IterationDiagram admits all options of Plot and the additional options LineStyle (default value: Blue) and PointStyle (default value: {Blue,PointSize[Medium]}.";

Unprotect[Jacobian];
Jacobian[f_][{var__}]:=Through[Through[(Derivative/@IdentityMatrix[Length[{var}]])[f]][{var}]]\[Transpose]
Jacobian[f_][var__]:=Through[Through[(Derivative@@@IdentityMatrix[Length[{var}]])[f]][var]]\[Transpose]
Protect[Jacobian];
Jacobian::usage="Jacobian[f][x,y,…] gives the Jacobian matrix {(∂_x)f[x,y,…],(∂_y)f[x,y,…],…}\[Transpose] for a function f:(x,y)↦f[x,y] of any number of variables.
Jacobian[f][{x,y,…}] gives the Jacobian matrix {(∂_x)f[{x,y,…}],(∂_y)f[{x,y,…}],…}\[Transpose] for a function f:{x,y}↦f[{x,y}] of any number of variables.
The parameters x,y,… need not be symbolic.";

Recursies in één dimensie:

Dit notebook kan als hulpmiddel worden gebruikt bij het bestuderen van recursierelaties.
Allereerst zetten we een paar nuttige standaardprocedures op een rijtje:

?N
?D
?Nest
?NestList
?Table
?Part
?Drop
?Take
?Plot
?ListPlot
?Simplify
?FullSimplify
?Solve
?FindRoot
?ReplaceAll
?Manipulate

Het model van Verhulst

Als voorbeeld beschouwen we het model van Verhulst dat in detail in sectie 1.3 van het dictaat is beschreven.
De bijbehordend recursierelatie luidt (Subscript[x, n+1])=(Subscript[f, a])((Subscript[x, n])), waarbij (Subscript[f, a]) gegeven wordt door (Subscript[f, a])(x)=a x(1-x):

f[a_][x_]:= a x(1-x)

f[α]
f[α][ξ]

Met NestList kunnen we makkelijk snel voor a=2 en (Subscript[x, 0])=3/4 de eerste twintig iteraties (Subscript[x, 0]),(Subscript[x, 1])…,(Subscript[x, 20]) verkrijgen en deze vervolgens in een grafiek uitzetten:
Let op: het punt (Subscript[x, 0])wordt automatisch geplaatst bij n=1 geplaatst. Dit kan worden opgelost door Mathematica expliciet te vertellen dat n van 0 tot 20 loopt.

NestList[f[2.5],.75,20]
ListPlot[%,AxesLabel->{"n","(Subscript[x, n])"},PlotRange->{0,1},PlotStyle->{PointSize[Medium]}]
ListPlot[%%,AxesLabel->{"n","(Subscript[x, n])"},PlotRange->{0,1},PlotStyle->{PointSize[Medium]},DataRange->{0,Length[%]-1}]

Wees voor zichtig wanneer je Nest of NestList gebruikt: Mathematica probeert wanneer mogelijk uitkomsten exact te bepalen, wat niet altijd wenselijk is.
Een exacte invoer leidt tot een exact resultaat:

Nest[f[2],3/4,10]

Dit zorgt ervoor dat berekeningen vaak onnodig lang duren:

NestList[f[2],3/4,25]//Timing//First

Deze berekening gaat een stuk sneller als minstens één van de argumenten niet exact is.

NestList[f[2.],3/4,25]//Timing//First
NestList[f[2],.75,25]//Timing//First
NestList[f[2],N[3/4],25]//Timing//First
NestList[Composition[N,f[2]],3/4,25]//Timing//First

Het commando Iterations dat in de preamble is gedefinieerd, heeft dezelfde functionaliteit als NestList, maar past eerst voor de zekerheid N toe op (Subscript[x, 0]).

Iterations[f[2],3/4,10000]//Timing//First

??Iterations

Dekpunten

De dekpunten voor ons model kunnen we vinden met behulp van Solve:

Dekpunten[a_]=Solve[f[a][x]==x,x]

Ter controle kunnen we (Subscript[f, a]) evalueren in deze punten door /. (ReplaceAll) te gebruiken:

f[a][x]/.Dekpunten[a]//Simplify

De afgeleide van (Subscript[f, a]) in de dekpunten is vervolgens ook eenvoudig bepalen.

D[f[a][x],x]/.Dekpunten[a]//Simplify

Een goede manier om een overzicht te krijgen van de posities van de dekpunten voor een recursie is het plotten van een grafiek:

Plot[{x,f[2.5][x]},{x,0,1},AxesLabel->{"x","(Subscript[f, 5/2])(x)"},PlotRange->{0,1}]

Mocht de vergelijking (Subscript[f, a])(x)=x niet exact kunnen worden opgelost dan kan FindRoot gebruikt worden in plaats van Solve, mits alle parameters expliciet zijn gegeven.
FindRoot zal geen exacte uitvoer geven.
N.B. FindRoot geeft elke keer slechts één dekpunt, er kunnen er meer zijn!

FindRoot[f[2.5][x]==x,{x,1}]
FindRoot[f[2.5][x]==x,{x,.1}]

Iteratiediagrammen

Het command IterationDiagram, dat in de preamble gedefinieerd is toont in de grafiek voor de functie f zien hoe de iteratie verloopt:

IterationDiagram[f[3.2],{0,1},1/4,100,AxesLabel->{"x","f(x)"}]
IterationDiagram[f[3.2][x],{x,0,1},1/4,100,AxesLabel->{"x","f(x)"}]

?IterationDiagram

Bifurcatiediagrammen

Het commando BifurcationPlot, dat in de preamble gedefinieerd is, kan gebruikt worden om bifurcatiediagrammen te produceren.

BifurcationPlot[f[a],.5,{200,300},{a,0,4,.01},AxesLabel->{"a","x"},PlotStyle->{PointSize[Small]}]

?BifurcationPlot

Manipulate

Voor het bestuderen van de afhankelijkheid van verschillende eigenschappen van een recursie van parameters is Manipulate erg handig:

Manipulate[ListPlot[NestList[f[a],N[x0],25],AxesLabel->{"n","x"},PlotRange->{0,1}],{{a,2},0,4},{{x0,.75,"(Subscript[x, 0])"},0,1}]

Manipulate[Plot[{x,Nest[f[a],N[x],n]},{x,0,1},AxesLabel->{"x","((Subscript[f, a]))^n(x)"},PlotRange->{0,1}],{{a,2},0,4},{n,Range[1,10],ControlType->SetterBar}]

Manipulate[IterationDiagram[f[a],{0,1},x0,n,AxesLabel->{"x","f(x)"}],{{a,2.},0,4},{{x0,.75},0,1},{{n,100},{1,10,100,1000}}]

Manipulate[StringForm["Voor a=`` zijn de dekpunten `` en zijn de afgeleides ``.",NumberForm[a//N,{∞,2}],NumberForm[f[a][x]/.Dekpunten[a]//N,{∞,2}],NumberForm[f[a]'[x]/.Dekpunten[a],{∞,2}]]//Quiet,{{a,2},0,4}]

Recursies in twee dimensies:

Het commando Recursion2D, dat in de preamble is gedefinieerd, kan nuttig zijn bij het bestuderen van algemene recursies in (ℝ^2).

?Recursion2D

Let op: aangezien een punt x∈(ℝ^2) er in Mathematica uitziet als {x,y}, moet de functie f waarvoor (Subscript[x, n+1])=f((Subscript[x, n])) worden ingevoerd als f[{x_,y_}]:= {...,...}.

Een lineair voorbeeld:

Deze recursie bestaat uit herhaalde vermenigvuldiging van een vector met de matrix ρ (({{cos(θ), -sin(θ)}, {sin(θ), cos(θ)}})).

g[ρ_,θ_][{x_,y_}]:=ρ RotationMatrix[θ].{x,y}

Met NestList kunnen we een baan voor deze recursie (exact) produceren:

NestList[g[(Sqrt[2]),π/4],{1,0},16]

We kunnen eenvoudig controleren dat g[ρ,θ] zijn eigen Jacobiaan is:

D[g[ρ,θ][{x,y}],{{x,y}}]//MatrixForm

ListPlot werkt ook voor twee-dimensionale recursies en is dus bruikbaar om recursies in twee dimensies weer te geven:

ListPlot[Iterations[g[0.9,π/8],{1,0},100],AxesLabel->{"x","y"},PlotStyle->PointSize[Medium],PlotRange->{{-1.1,1.1},{-1.1,1.1}},AspectRatio->1]

Met Recursion2D kunnen we gemakkelijk zien wat er gebeurt als de paramater ρ en θ worden gevarieerd:

Recursion2D[g[ρ,2π ϑ],{2,2},{{-4,4},{-4,4}},100,{{ρ,.95},2/3,3/2},{{ϑ,1/12,"θ/2π"},0,1/8},AspectRatio->1,Joined->True,InterpolationOrder->2,Mesh->Full]

Rupsen en sluipwespen:

Dit is het model voor de gecombineerde populatie van rupsen en sluipwespen uit sectie 3.1

rupswesp[a_,b_][{x_,y_}]:={(a-x-y)x,b x y}

D.m.v. Solve kunnen we hiervoor eenvoudig de dekpunten vinden:

Dekpunten[a_,b_]=Solve[rupswesp[a,b][{x,y}]=={x,y},{x,y}]

Voor de stabiliteit van deze dekpunten beschouwen we de Jacobiaan en de bijbehorende eigenwaarden:

D[rupswesp[a,b][{x,y}],{{x,y}}]//MatrixForm
Eigenvalues[%]

We kunnen /. (ReplaceAll) gebruiken om dit te evalueren in de dekpunten. Omdat er drie dekpunten zijn levert dit een lijst van drie matrices op, één voor ieder dekpunt:

D[rupswesp[a,b][{x,y}],{{x,y}}]/.Dekpunten[a,b]//Simplify
MatrixForm/@%
Eigenvalues/@%%

Deze eigenwaarden kunnen gebruikt worden om de stabiliteit van de betreffende dekpunten te analyseren.
N.B. In de bovenstaande cel wordt Map (/@) gebruikt om MatrixForm resp. Eigenvalues toe te passen op ieder van de matrices afzonderlijk.
Over het algemeen geldt dat f/@{a,b,c,…} de lijst {f[a],f[b],f[c],…} als uitvoer geeft.

?Map

Door Recursion2D (of gewoon ListPlot) te gebruiken krijgen we een goed zicht op hoe de populaties zich daadwerkelijk ontwikkelen:

Recursion2D[rupswesp[a,b],{1,1},{{0,3.5},{0,3.5}},{x,y}/.Dekpunten[a,b],{25,5000},{{a,3.35},0,5},{{b,1.1},0,3},AspectRatio->1,Locator->True,PlotStyle->PointSize[Tiny],
AxesLabel->{"x=r/R","y=w/W"},InitialPointName->"{(Subscript[x, 0]),(Subscript[y, 0])}"]