10.1 “Brute force”Attacke
10.3 Stationäre Lösung und Richtungsfeld
Wennp˙ undf˙gleichzeitig verschwinden in 10.1, ändern sich die Populationen nicht mehr!
0 = a·p −b·p·f Substitution in die obere Gleichung ergibt
0 =bd
c f−bd b ac f ·
f
Dies führt zum trivialen Lösungspaar(p, f) = (0,0)und zur stationären Lösung
(ps, fs) =
Man kann sich leicht vergewissern, wenn man in obigen Programmen als Anfangsvektor d
c, a b
vorgibt, bleiben beide Populationen konstant.
Jede Populationsentwicklung ~x(t) := (p(t), f(t)) lässt sich auch als parametrisierte Kurve im (p, f)-Raum interpretieren. 10.1 gibt dann an, wie die Tangentenvektoren einer solchen Kurve auszuschauen haben.
~x˙ =F~(~x(t)) (10.4)
Zeichnet man in jedem Punkt des(p, f)-Raumes (Phasenraum) den Tangentenvektor ein (Rich-tungsfeld) bekommt man einen guten Überblick wie die Kurve (Trajektorie = Entwicklung der Population) ausschauen wird - sie ist ja jene Kurve, die die Vektoren als Tangente besitzt.
Man wählt also einen Startpunkt und folgt mit einem Stift einfach den Vektoren. Wobei wir durchaus die Vermutung hegen, dass die zyklischen Populationen “irgendwie” um die stationäre Lösung “kreisen” werden.
InGNU-Octave heißt der Befehl für das Zeichnen eines Richtungsfeldsquiver(~x, ~y, Fx, Fy) Dabei ist
n ~x (bzw.~y) ist der Stützstellenvektor inx (bzw.y)-Richtung n Fx (bzw.Fy) sind die Terme von 10.4 (rechte Seite)
Als Output bekommt man an den Stellen (x0, y0),(x1, y1), . . .(xn, yn) die Tangentenvektoren gezeichnet. Wir haben im folgenden GNU-Octave-Skript auch deutlich die stationäre Lösung eingezeichnet:
10.3 Stationäre Lösung und Richtungsfeld
a=1;
2 b =0.2;
c =0.04;
4 d =0.8;
f_s=a/b ; # s t a t i o n a r y s o l u t i o n f_s = 5
6 p_s=d/ c ; # s t a t i o n a r y s o l u t i o n p_s= 20
# l i n s p a c e ( lowerBound , upperBound , number o f i n t e r s e c t i o n −p o i n t s )
8 rangeP = l i n s p a c e (0 , 2∗p_s , 15) ; rangeF = l i n s p a c e (0 , 2∗ f_s , 15) ;
10# get a 2−dim g r i d f o r [ x , y ] −−> here [ p , f ] [ p , f ]=meshgrid( rangeP , rangeF ) ;
12# e x p r e s s i o n f o r change in x−d i r e c t i o n dp=a∗p−b∗p . ∗ f ;
14# e x p r e s s i o n f o r change in y−d i r e c t i o n df=c ∗p . ∗ f−d∗ f ;
16 hold on ;
a x i s ([ −1 45 −1 1 0 . 5 ] )
18 quiver(p , f , dp , df ) ;
p l o t(p_s , f_s , "o", " markersize ", 10 , " m a r k e r f a c e c o l o r "," red ") ;
20 t i t l e ("Tangent−Vectors ") ; hold o f f ;
Abb.73 : Richtungsfeld um die stationäre Lösung mitquiver
Deutlich kann man das Wirbelfeld um die stationäre Lösung erkennen (in der Physik spielen Wirbelfelder eine große Rolle, z.B. erzeugen zeitlich variable magnetische Wirbelfelder elektri-sche Wirbelfelder). Wir erkennen auch, dass die Vektoren nicht symmetrisch um die stationäre Lösung verteilt sind - welche Formen haben die Populationskurven?
10. Räuber-Beute Modell
In wxMaxima heißt der Befehl zur Darstellung eines Richtungsfelds (directionfield or slope-field) plotdf. Eingabeparameter sind: Liste der Ableitungsterme, Variablenlist, Trajektorie durch die Anfangswerte, Intervalle in denen die Variablen gezeichnet werden. Klickt man mit der Maus auf die Zeichnung wird dies als Ausgangspunkt einer Trajektorie gemacht (hier wurde in die Nähe der stationären Lösung geklickt). Das sieht dann bei uns so aus:
(%i1) (a:1,b:0.2,c:0.04,d:0.8)$
(%i2) dpdt:a*p-b*p*f$
(%i3) dfdt:c*p*f-d*f$
(%i6) plotdf([dpdt,dfdt],[p,f],[trajectory_at,50,10],[p,-1,75],[f,-1,15])$
Abb.74 : Richtungsfeld und Trajekorie (50,10) inxmaximamitplotdf
Damit der Befehl plotdf funktioniert, muss xmaxima installiert sein - plotdf kann man zwar vonwxMaxima aufrufen, schickt dann aber die Ausgabe immer an xmaxima.
10.3 Stationäre Lösung und Richtungsfeld
10.3.1 Richtungsfeld in Geogebra Befehl Slopefield
Bei der englischen Version vonGeogebra (die ich benutze) lautet der eingebaute Befehl für das Zeichnen eines Richtungsfeldes
SlopeField(<f(x,y)>,<Number n>,<Length Multiplier a>,
<Min x>,<Min y>,<Max x>,<Max y>)
Folgende Angaben werden (für die größtmögliche Kontrolle) benötigt:
n die Differentialgleichung dy
dx =f(x,y)- wir haben allerdings dx
dt und dy dt also müssen wir umformen dy
dx = dy dt
dt dx = dy
dt 1
dx dt
n das Richtungsfeld wird in einemn×ngrid (Gitter) angegeben
n die Länge der Steigungslinien wird mit dem nächsten Parameter bestimmt n der rechteckige Bereich (left bottom - right top) für das Richtungsfeld
Abb.75 : DerSlopefield - Befehl und seine Ausgabe
Nachteile:
* keine Orientierung des Feldes
* alle Linien gleich lang (Stärke nicht abschätzbar→ stationaäre Lösung kaum sichtbar)
* umschreiben der Gleichung bzw. Gleichungssystems
10. Räuber-Beute Modell
Selbstgebastelt - um zu lernen
n Mit Corner(1) bzw. Corner(3) holen wir uns die Koordinaten der linken unteren Ecke (lb - left bottom) bzw. der rechten oberen Ecke (rt - right top)
n xDensity - wieviele Punkte in x-Richtung wollen wir haben?
n DerscaleFactor verkürzt/verlängert die Vektoren im Vektorfeld n ∆xist die Abstand des Gitters inx-Richtung
n xGrid sind diex-Koordinaten des Gitters
n 7-9 dasselbe für diey-Richtung; 10 erzeugt die Gitterpunkte n 11-16 erzeugen die Änderungsraten für Füchse bzw. Pinguine
n h1 sind Tangentenvektoren im Phasenraum (Pinguin-Füchse) an den Gitterpunkten n Wir setzen das Zeit-Inkrement (so klein als nötig!?) fest
n 19 - Wir versetzen die Gitterpunkte um∆t·h1 → grid2
Abb.76 : Hier das Konstruktionsprotokoll
n VF bzw. VFs sind jetzt das Vektorfeld bzw.
skaliertes Vektorfeld
die restlichen Befehle zeichnen den Zyklus im Phasenraum Pinguine - Füchse
Abb.77 : Fortsetzung
n 23/24 sind wieder unsere Differentialgleichungen
10.3 Stationäre Lösung und Richtungsfeld
n PunktA ist ein frei wählbarer Anfangspunkt(p0, f0) für Pinguine und Füchse
n Wir lösen das Differentialgleichungssystem numerisch und bekommen die Ortskurven für Pinguine und Füchse in Abhängigkeit von t in [0,9] - das ist leider nicht das was wir brauchen und darum müssen wir uns noch kümmern!
n Anfangszeitpunkt t0 und Endzeitpunkt tf werden gesetzt - das wäre vorher günstiger gewesen, aber was soll’s!
n Jetzt “samplen” wir die Ortslinien für Pinguine und Füchse indem wir sie “scannen” - im Protokoll wird das mit “Freehand function” bzw. “freehand(x)” bezeichnet, was die Sache nicht ganz trifft, denn die Befehle lauten
g(x)=Function(Join({{t_0, t_f}, y(First(penguins, Length(penguins)))})) f(x)=Function(Join({{t_0, t_f}, y(First(foxes, Length(foxes)))}))
y(First(penguins, Length(penguins)) → y-Werte der Pinguinpunkte
Function - benötigt t0 undtf und Werte - alles in 1 Liste (wird mitJoin erledigt) Curve(g(t), f(t), t, t_0, t_f) zeichnet die Kurve im Phasenraum!
Durch ziehen vonA können die Anfanggsbedingungen gewählt werden!
Abb.78 : Richtungsfeld und Trajekorie durch PunktAinGeogebra
10. Räuber-Beute Modell
Wie man in Abb. 78 sieht, sind die Richtungsvektoren in der Nähe des stationären Punktes recht kurz, dies könnte man verbessern, indem man diese Vektoren (falls sie eine gewisse Länge unterschreiten) verlängert.
Sowohl dieser “Grenzwert”(sliderthreshold) als auch der Verlängerungsfaktor(sliderampFactor) sollen eingestellt werden können. Natürlich muss man diese Vektoren farblich hervorheben.
Abb.79 : Zusätzliche Befehle für das “Zoomen”
34: grid3 - man trägt von den Gitterpunkten(grid) die skalierten Vektoren(VFs) auf 35: Längenmessung zwischen grid und grid3
36: Wir konstruieren Punkte Pi(di, i) -di ist die Länge vom i-ten Vektorgrid→grid3, falls kleiner alsthreshold - wird der Index als y-Koordinate gespeichert (Null sonst).
37: Welche Indices haben die “kurzen” Vektoren und diese werden bei 38 verlängert!
Abb.80 : Kurze Richtungsvektoren werden verlängert
10.3 Stationäre Lösung und Richtungsfeld
10.3.2 Stützpunkte der Phasenraumkurve
Ein interessantes Zusatzproblem: Greife von der Lösungskurve im Phasenraum äquidistante Stützpunkte heraus (z.B. für eine Splinefunktion). (Anmerkung:günstiger wäre es natürlich bei starker Krümmung die Dichte der Stützpunkte zu erhöhen). Da man dabei viel überGeogebra -Skript lernen kann, möchte ich das hier zeigen.
Ein weiterer kritischer Punkt ist die Periodendauer festzustellen - dies wird normalerweise mit einer Fourieranalyse bzw. Autokorrelation einer Datenreihe (die bei uns hierNSolveODE liefert) bewerkstelligt. Weil dies sehr aufwendig ist, behelfen wir uns mit einem Slider für tf
(Ende des Lösungsintervalls) und vergrößern (ev. verkleinern) tf solange, bis im Phasenraum ein Zyklus beendet ist - dies hängt natürlich von den Parametern ab!
Abb.81 :tf ist noch zu klein - Kurve im Phasenraum noch offen!
Abb.82 : Letzter Stützpunkt stimmt (fast) mit AnfangspunktAüberein
10. Räuber-Beute Modell
Um das Programm nicht noch mehr aufzublähen, lassen wir jetzt das Richtungsfeld weg.
Nun zur Erklärung des Programms und wozu die seltsame Kurve in Graphics2 dient.
Abb.83 : Lösung des ODE-Systems
Bis auf Zeile 8 (ein Slider zum Einstellen der Anzahl der Stützpunkte, die äquidistant vom Anfangspunkt entfernt sind) sind die Zeilen 1-13 sind lauter “alte Bekannte”.
Hier die neuen Befehle:
Abb.84 : Konstruktion der äquidistanten Stützpunkte
14-15 Wir holen uns die Punkte von den Ortslinien - davon diey-Werte.
Es ist wichtig, dass der Slider tf so eingestellt wird, dass genau(!) 1 Zyklus komplett ist -nicht mehr und -nicht weniger!
16 Wir konstruieren die Punkte im Phasenraum(p, f)(penguins and foxes) - diese liegen aber verschieden weit auseinander (zoomen!)
17 di =
(pi, fi)−(pi−1, fi−1)
∀i∈ {2,3, . . . N} - Entfernung der Nachbarpunkte 18 Wir konstruieren eine Folge von Punkten:(D0,0),(D1,1),(D2,2), . . .(DN, N)
wobeiDk=Pk
i=1di - also der Abstand desk-ten Punktes vom Anfangspunkt.
Es gilt natürlich die Rekursion:D0 = 0 und Dk=Dk−1+dk Und obige Punkte (Dk, k)bilden die Kurve in Graphics2:
x-Werte entsprechen dem Abstand vom Anfangspunkt (auf der Kurve gemessen!) diey-Werte sind der Index des Punktes in der(p, f)-Liste.
Dass diese Kurve keine Gerade ist, zeigt die variierende Dichte der Punkte
Punkte mit gleichem Abstand sind in diesem monoton steigendem Graph Punkte, derenx-Werte sich mit gleichem Wert unterscheiden. Weitere Vorgangsweise:
20,21 GesamtlängeL holen, `n=L/n nist die Anzahl der Stützpunkte(T otP nt) 22 Wir schneiden die Senkrechten x=i·`n i∈ {1,2, . . . n} mit
Polyline von L2Index, nehmen deny-Wert und runden diesen auf eine
ganze Zahl (nächster Index) und bekommen eine Liste der Stützpunkt-Indices 23 Diese Punkte holen wir uns aus dem(p, f)-Feld
In unserer Abb. 82 sind offensichtlich die 2 Stützunkte, die sich am weitesten links im Zyklus befinden weiter auseinander als die anderen. Wenn Sie das Achsenverhältnis auf 1:1 schalten, klärt sich dieser Sachverhalt auf!