Füchse und Kaninchen

Die Entwicklung einer Fuchs‑Kaninchen‑Population:

Ein Paar gewöhnlicher nichtlinearer Differentialgleichungen 1. Ordnung.

von Dr. Rolf Schröder

Inhalt

0. Prolog
1. Die Änderung der Kaninchenzahl K mit der Zeit
2. Die Änderung der Fuchszahl F mit der Zeit
3. Das DGL‑System für Füchse und Kaninchen
4. Stationäre Lösung bzw. Ruhepunkt
5. Die Charakteristik Φ(K,F)
6. Diskussion der Charakteristik Φ(K,F)
7. Analytische Lösung im Grenzfall F → F_s, K → K_s
8. Interaktive graphische Lösung des DGL‑Systems

zurück

0. Prolog

Wer hat in freier Natur heutzutage schon mal einen Fuchs gesehen! Dagegen sehen wir gelegentlich Kaninchen, und sei es nur bei einer nächtlichen Autofahrt. Werden es von Jahr zu Jahr mehr, weil wir keine Füchse mehr haben? Nein, auch ohne Füchse nimmt die Zahl der Kaninchen bei uns offensichtlich nicht unbegrenzt zu. Das, was wir von einer angemessenen Fuchspopulation erwarten könnten, nämlich die natürliche Beschränkung der Kaninchenpopulation, wurde mit anderen Mitteln erreicht, wie etwa mit der bewußt gegen Kaninchen eingesetzten Viruserkrankung Myxomatose.

Wir wollen uns nun in ein „Schlaraffenland für Kaninchen“ versetzt denken, ohne Krankheiten und mit unbegrenztem Nahrungsangebot. Irgendjemand, ich weiß nicht mehr wer, ist dann auf die Idee gekommen, daß auch Füchse in diesem Land ideale Lebensbedingungen finden könnten. Und so kamen gegen den heftigen Protest der Kaninchen auch einige Füchse in das Schlaraffenland der Kaninchen. Schon bald zeigte sich die Auswirkung dieser Maßnahme: unter den zunächst zahlreichen Kaninchen fand die anfangs kleine Schar der Füchse leichte Beute, und sie vergrößerte sich schnell. Wie aber geht die Entwicklung der beiden Gruppen wohl weiter? Wächst die Zahl der Kaninchen und Füchse unbegrenzt? Oder vernichten die Füchse schließlich alle Kaninchen und sterben dann auch aus? Oder stellt sich ein Gleichgewicht ein?

Wir wollen ein einfaches Modell aufstellen, um die Entwicklung der Fuchs‑Kaninchen‑Population zu analysieren.

zurück

1. Die Änderung der Kaninchenzahl K mit der Zeit

Da wir ein unbegrenztes Nahrungsangebot für die Kaninchen annehmen, würde die Zahl der Kaninchen in einer kleinen Zeitspanne dt um dK zunehmen, und zwar proportional zur gerade vorhandenen Kaninchenzahl K, in Formeln:

dK = A_K∙K∙dt ,

dK/dt = A_K∙K ,

mit dem positiven Koeffizienten Ak, der die Stärke der Zunahme beschreibt. Diese Gleichung läßt sich durch die sog. Trennung der Variablen lösen, d. h., wir schreiben die Gleichung so um, daß sich die beiden Variablen K und t auf jeweils einer Seite befinden:

dK/K = A_K∙dt ,

dLn(K) = A_K∙dt ,

und dies ergibt sich offensichtlich aus der Ableitung der Funktion

Ln(K) = A_K∙t + C ,

C eine Konstante.

Wendet man auf beide Seiten der Gleichung die Exponentialfunktion an, so folgt:

K = e^{A_K∙t + C} = e^C∙e^A_K∙t = K₀∙e^A_K∙t ,

mit K₀ als Zahl der Kaninchen zum Zeitpunkt t = 0. Da A_K > 0 ist, nimmt die Kaninchenpopulation also ständig exponentiell zu. Der reziproke Wert der Konstanten A_K hat offenbar die Dimension einer Zeit. Wir können daher statt A_K auch die Zeitkonstante

τ_K ≡ 1/A_K

benutzen. Sie gibt die Zeit an, in der sich die Anzahl der Kaninchen um den Faktor e vergrößert. Damit ist die Entwicklung der Kaninchenpopulation ohne Beeinträchtigung durch Füchse beschrieben.

Aber durch die Anwesenheit der Füchse wird nun das Wachstum der Kaninchen verringert, denn in einer kleinen Zeitspanne dt fressen die Füchse eine gewisse Zahl dK von Kaninchen: diese Zahl ist natürlich proportional zur Zahl F der Füchse, sie ist aber auch zur Zahl K der Kaninchen proportional, da mit wachsender Kaninchenzahl auch die Möglichkeit für die Füchse wächst, Kaninchen zu fangen. In der Zeit dt wird also die Zahl der Kaninchen um B_K∙K∙F reduziert, B_K beschreibt, wie stark diese Reduktion pro Zeit ist. Wir berücksichtigen daher in der anfänglichen Gleichung zusätzlich die Abnahme pro Zeitintervall dt:

dK = A_K∙K∙dt + B_K∙K∙F∙dt ,

dK/dt = A_K∙K + B_K∙K∙F ,

A_K > 0, B_K < 0 .

Man beachte, daß B_K negativ ist und daher eine Abnahme bedeutet.

zurück

2. Die Änderung der Fuchszahl F mit der Zeit

Wir betrachten zuerst die Füchse ohne Kaninchen. Da Kaninchen die einzige Nahrung der Füchse sind, hätte ein plötzlicher Ausfall dieser Nahrungsquelle für die Füchse dramatische Folgen, die Fuchspopulation müßte ständig abnehmen. Wir ersparen uns die detaillierte Beschreibung des Fuchssterbens und nehmen an, daß der Verlust dF an Füchsen pro Zeitintervall dt proportional zur noch vorhandenen Fuchszahl F selbst ist, in Formeln:

dF = A_F∙F∙dt ,

dF/dt = A_F∙F ,

mit dem negativen Koeffizienten A_F, der die Stärke der Abnahme beschreibt.

Ganz analog wie im Falle der Kaninchen läßt sich diese Gleichung lösen, die Lösung lautet entsprechend:

F = e^{A_F∙t + C} = e^C∙e^A_F∙t = F₀∙e^A_F∙t ,

mit F₀ als Zahl der Füchse zum Zeitpunkt t = 0. Da A_F < 0 ist, nimmt die Fuchspopulation also ständig exponentiell ab. Auch hier kann man statt der Konstanten A_F die reziproke Konstante

τ_F ≡ 1/A_F

benutzen, sie hat die Dimension einer Zeit. Da mit A_F auch τ_F negativ ist, beschreibt |τ_F| die Zeit, in der die Fuchspopulation auf e⁻¹ abgefallen ist.

Zum Glück für die Füchse gibt es nun die Kaninchen! Je mehr Kaninchen die Füchse fangen können, umso rascher können sie ihre Abnahme durch Zuwachs ausgleichen. Dieser Zuwachs wird im einfachsten Fall zur Zahl der Füchse aber auch zur Zahl der Kaninchen proportional sein. In der Zeit dt wird also die Zahl der Füchse um B_F∙K∙F wachsen, B_F beschreibt, wie stark dieser Zuwachs pro Zeit ist. Wir berücksichtigen daher in der anfänglichen Fuchsgleichung zusätzlich die Zunahme pro Zeitintervall dt:

dF = A_F∙F∙dt + B_F∙K∙F∙dt ,

dF/dt = A_F∙F + B_F∙K∙F ,

A_F < 0, B_F > 0 .

zurück

3. Das DGL‑System für Füchse und Kaninchen

Wir haben die beiden folgenden Differentialgleichungen zur Beschreibung einer Fuchs‑Kaninchen‑Population aufgestellt:

dK/dt = A_K∙K + B_K∙K∙F ,

dF/dt = A_F∙F + B_F∙K∙F ,

A_K > 0, B_K < 0, A_F < 0, B_F > 0 .

Die Anzahlen K und F sind reine Funktionen der Zeit t:

K = K(t), F = F(t) .

Zur Nomenklatur:

Danach haben wir ein gewöhnliches explizites nichtlineares autonomes Differentialgleichungssystem 1. Ordnung vor uns.

Wir wollen das System noch unter einem anderen Gesichtspunkt betrachten. Aus Gewohnheit schreiben wir statt der Funktionen K(t) und F(t) lieber X(t) und Y(t). Allgemeiner können wir dann unser System auch so schreiben:

dX/dt = f(X,Y) ,

dY/dt = g(X,Y) ,

mit den über den ganzen X‑Y‑Raum stetigen Funktionen f und g. Zur Veranschaulichung betrachten wir ein kartesisches Koordinatensystem in der Ebene. X und Y beschreiben den Ort und dX/dt und dY/dt stellen dann eine Geschwindigkeit am Ort X,Y dar. Wenn f(X,Y) und g(X,Y) in der ganzen Ebene erklärt sind, so gibt es also auch für jeden Punkt eine bestimmte von der Zeit unabhängige Geschwindigkeit, das ist das Bild einer stationären Strömung!

Wir können offenbar wegen der bekannten Geschwindigkeit den Nachbarort zum Punkt X,Y, der in einer kleinen Zeitspanne Δt erreicht wird, bestimmen. Wir schreiten entsprechend der Geschwindigkeit dX/dt, dY/dt um ein Wegstück ΔX = (dX/dt)∙Δt, ΔY = (dY/dt)∙Δt weiter. Damit haben wir eine praktisch anwendbare Lösungsvorschrift. Setzen wir uns zur Zeit t₀ in einen willkürlich gewählten Ort X₀,Y₀, so finden wir den Ort X₁,Y₁ nach einem Zeitschritt Δt, indem wir das Wegstück ΔX = f(X₀,Y₀)∙Δt, ΔY = g(X₀,Y₀)∙Δt zu X₀,Y₀ addieren, so daß wir den Punkt X₁,Y₁ erreichen:

X₁ = X₀ + (dX/dt)₀∙Δt = X₀ + f(X₀,Y₀)∙Δt ,

Y₁ = Y₀ + (dY/dt)₀∙Δt = Y₀ + g(X₀,Y₀)∙Δt ,

Dann wiederholen wir diese Prozedur, indem wir X₀,Y₀ durch X₁,Y₁ ersetzen, und erhalten X₂,Y₂ usw. In der X‑Y‑Ebene finden wir so eine Kurve, die offenbar durch den „Anfangspunkt“ oder auch „Anfangswert“ X₀,Y₀ bestimmt ist. Da wir auch die Zeitschritte kennen, haben wir auch X(t) und Y(t) bestimmt. Statt vom Punkt X₀,Y₀ in der Zeit voranzuschreiten, können wir aber auch zurückschreiten, indem wir die Zeitschritte negativ nehmen. D. h., zu jedem Zeitpunkt (vor und nach t₀) ist eine Kurve durch X₀,Y₀ bestimmt (jedoch nicht notwendig die wahre Lösungskurve!).

Für die exakte Lösungskurve müssen die Zeitschritte Δt im Grenzfall gegen Null gehen (und damit wegen der Stetigkeit von f und g auch ΔX und ΔY). Da man beim numerischen Rechnen aber stets mit endlichen Größen arbeiten muß, ist es zweckmäßig, f und g im jeweiligen Punkt X,Y mindestens bis zur zweiten Ordnung nach der Zeit in eine Taylorreihe zu entwickeln, um damit dann mit dem (endlichen) Zeitschritt Δt die nächsten Werte ΔX, ΔY zu berechnen. Die Genauigkeit der Rechnung steigt hierdurch ganz beträchtlich und man kommt der wahren Lösungskurve sehr nahe. Ein bekanntes Verfahren zur numerischen Lösung solcher „Anfangswertprobleme“ ist das Runge‑Kutta‑Verfahren (Fehler pro Schritt ∼ (Δt)⁵), das auch wir später einsetzen werden. Zuvor wollen wir noch einige Eigenschaften der Lösungskurven X(t) und Y(t) bzw. der Kurve (X(t),Y(t)) in der X‑Y‑Ebene bestimmen.

Die Kurve, welche der Punkt X(t),Y(t) – ausgehend von X₀,Y₀ – mit der Zeit als Parameter durchläuft, nennt man „Charakteristik“ des Differentialgleichungssystems. In unserem Fall gibt es eine Schar von Charakteristiken.

Sind die Funktionen f und g eindeutig, dann kann eine Charakteristik sich nicht mit sich selbst schneiden. In einem Schnittpunkt müßten andernfalls zwei verschiedene Geschwindigkeiten vorkommen, d. h., f bzw. g wären mehrdeutig.

Kehren wir wieder zu unseren Kaninchen und Füchsen und damit zum DGL‑System zurück! Obwohl unser DGL‑System analytisch einfache Funktionen f und g hat, lassen sich die Lösungen K(t) und F(t) nicht durch bekannte Funktionen beschreiben, jedenfalls sind keine Lösungen K(t),F(t) in geschlossener Form bekannt. Jedoch können wir gewisse Eigenschaften feststellen.

zurück

4. Stationäre Lösung bzw. Ruhepunkt

Wenn die Zeitableitungen dK/dt und dF/dt bei einer bestimmten Population K_s,F_s (der Index s stehe für stationär) gleich Null sind, so bleiben K(t) und F(t) für alle Zeiten konstant. Eine (triviale) Lösung ist offenbar:

K_s = 0 , F_s = 0 .

Da aber eine Welt ohne Füchse und Kaninchen recht langweilig ist, nehmen wir nun stets an, daß gilt:

K(t) ≠ 0 , F(t) ≠ 0 .

Wir finden dann eine Lösung ungleich Null, indem wir die linke Seite des DGL‑Systems gleich Null setzen und nach K und F auflösen:

0 = A_K∙K + B_K∙K∙F ,

0 = A_F∙F + B_F∙K∙F ,

A_K ≠ 0, A_F ≠ 0, B_K ≠ 0 , B_F ≠ 0 .

Als einzige stationäre Lösung, bei der also die Anzahl der Füchse und Kaninchen konstant bleibt, erhalten wir:

K_s = −A_F/B_F , F_s = −A_K/B_K ,

A_K ≠ 0 , A_F ≠ 0 , B_K ≠ 0 , B_F ≠ 0 .

zurück

5. Die Charakteristik Φ(K,F)

Zum Auffinden der Charakteristik Φ(K,F)  bedienen wir uns eines kleinen Tricks: wir multiplizieren die erste Gleichung unseres DGL‑Systems, die Kaninchengleichung, mit dF/dt und die zweite Gleichung, die Fuchsgleichung, mit dK/dt und setzen gleich:

(dF/dt)∙(dK/dt) = A_K∙K∙(dF/dt) + B_K∙K∙F∙(dF/dt) ,

(dF/dt)∙(dK/dt) = A_F∙F∙(dK/dt) + B_F∙K∙F∙(dK/dt) ,

es folgt:

A_K∙K∙(dF/dt) + B_K∙K∙F∙(dF/dt) = A_F∙F∙(dK/dt) + B_F∙K∙F∙(dK/dt) .

Wir dividieren mit K∙F und erhalten:

A_K∙(dF/dt)/F + B_K∙(dF/dt) = A_F∙(dK/dt)/K + B_F∙(dK/dt) .

Da nun d(ln(K))/dt = (dK/dt)/K bzw. d(ln(F))/dt = (dF/dt)/F ist, haben wir schließlich:

A_K∙d(ln(F))/dt + B_K∙dF/dt = A_F∙d(ln(K))/dt + B_F∙dK/dt .

Diese Gleichung läßt sich integrieren, man erhält:

A_K∙ln(F) + B_K∙F = A_F∙ln(K) + B_F∙K + C ,

C eine Konstante. Damit haben wir die Gleichung der Charakteristik gefunden:

Φ(K,F) = A_K∙ln(F) + B_K∙F − A_F∙ln(K) − B_F∙K − C .

zurück

6. Diskussion der Charakteristik Φ(K,F)

Bestimmen wir die Extrema von Φ(K,F): notwendige Bedingung ist, daß die partiellen Ableitungen nach K und F beide gleich Null sind:

∂Φ/∂K = −A_F/K − B_F = 0 ,

∂Φ/∂F =  A_K/F + B_K = 0 .

Es zeigt sich, daß die einzige Lösung für ein Extremum für alle K,F > 0 gerade die stationäre Lösung des DGL‑Systems ist:

K_extrem = −A_F/B_F = K_s ,

F_extrem = −A_K/B_K = F_s ,

Φ_extrem = Φ(K_s,F_s) = A_K∙{ln(−A_K/B_K) − 1} − A_F∙{ln(−A_F/B_F) − 1} − C .

Die Bedingung für das Vorliegen einem Maximums lautet:

(∂²Φ/∂K²)∙(∂²Φ/∂F²) − ∂²Φ/(∂K∙∂F) > 0 .

Die partiellen Ableitungen zweiter Ordnung lauten:

∂²Φ/∂K² = A_F/K² ,

∂²Φ/∂F² = −A_K/F² ,

∂²Φ/(∂K∙∂F) = 0 ,

und man stellt fest, daß die Bedingung für ein Maximum für alle K,F > 0 erfüllt ist:

(∂²Φ/∂K²)∙(∂²Φ/∂F²) − ∂²Φ/(∂K∙∂F) = −A_K∙A_F/(K²∙F²) > 0 .

Der Extremwert Φ_extrem ist also ein absolutes Maximum Φ_max:

Φ_max = A_K∙{ln(−A_K/B_K) − 1} − A_F∙{ln(−A_F/B_F) − 1} − C .

Da ferner an den Rändern des K‑F‑Bereiches die Funktion Φ(K,F) gegen −∞ geht, stellen die Funktionen Φ(K,F) = C' < Φ_max geschlossene Kurven um den Ruhepunkt K_s,F_s dar, die Zahl der Füchse und Kaninchen oszilliert um den Ruhepunkt K_s,F_s.

Das genaue Zeitverhalten läßt sich durch numerische Berechnung entlang der Charakteristik bestimmen: für jeden Kurvenpunkt kann man ja die zeitliche Änderung berechnen. Eine geschlossen angebbare Lösung ist aber nicht bekannt.

Ändert man relativ zueinander die Vorzeichen der Koeffizienten des DGL‑Systems, so erhält man die (nicht mehr zyklischen) Lösungen anderer Probleme wie z. B. „Symbiose“ oder „Totaler Krieg“.

zurück

7. Analytische Lösung im Grenzfall F → F_s, K → K_s

Um eine Lösung in der Umgebung des stationären Punktes zu bestimmen, führt man eine Variablentransformation durch:

K = K_s + x , F = F_s + y .

Damit geht das DGL‑System über in:

dx/dt = A_K∙(x + K_s) + B_K∙(x + K_s)∙(y + F_s) ,

dy/dt = A_F∙(y + F_s) + B_F∙(x + K_s)∙(y + F_s) ,

und nach wenigen Umstellungen erhalten wir:

dx/dt = (A_K + B_K∙F_s)∙x + B_K∙K_s∙y + B_K∙x∙y + K_s∙(A_K + B_K∙F_s) ,

dy/dt = B_F∙F_s∙x + (A_F + B_F∙K_s)∙y + B_F∙x∙y + F_s∙(A_F + B_F∙K_s) ,

und da die beiden Terme A_K+B_K∙F_s und A_F+B_F∙K_s gleich Null sind (siehe unter 4. Punkt), folgt:

dx/dt = B_K∙K_s∙y + B_K∙x∙y ,

dy/dt = B_F∙F_s∙x + B_F∙x∙y .

Wenn wir uns nun dem stationären Punkt K_s,F_s nähern, dann wird der Term mit dem gemischten Glied x∙y gegenüber den Gliedern mit x und y vernachlässigbar klein und wir erhalten die sog. lineare Näherung:

dx/dt = B_K∙K_s∙y ,

dy/dt = B_F∙F_s∙x ,

|x| ≪ K_s , |y| ≪ F_s .

Zur Lösung leiten wir die erste Gleichung einmal ab und setzen dann die zweite Gleichung ein, wir erhalten:

d²x/dt² = B_K∙B_F∙K_s∙F_s∙x .

Die allgemeine Lösung dieser DGL 2. Ordnung lautet:

x(t) = a∙sin(ω∙t + φ) ,

mit beliebigen Konstanten a und φ. Die zweimalige Ableitung ergibt:

d²x/dt² = −ω²∙x(t) ,

so daß also folgt:

ω = (−B_K∙B_F∙K_s∙F_s)^½ = (2∙π)/P .

Auch das negativ Vorzeichen für ω ist mathematisch eine Lösung, sie stellt eine zeitlich rückwärtslaufende Änderung dar, die uns aber nicht weiter interessiert. Die Periode P der oszillierenden Lösung ist jedenfalls:

P = 2∙π∙(−B_K∙B_F∙K_s∙F_s)^−½ ,

B_K∙B_F < 0 ,

und ersetzt man in diesem Ausdruck K_s und F_s mit Benutzung der Zeitkonstanten τ_K und τ_F, so zeigt sich, daß die Periode P allein durch die Zeitkonstanten beider Populationen beschrieben wird:

P = 2∙π∙(−τ_K∙τ_F)^½ .

Die Lösung für x(t) lautet also:

x(t) = a∙sin(2∙π∙t/P + φ) .

Leitet man diese Gleichung nach t ab und setzt sie in die erste DGL dieses Abschnitts ein, so erhalten wir y(t):

y(t) = (dx/dt)/(B_K∙K_s) = {2∙π/(P∙B_K∙K_s)}∙a∙cos(2∙π∙t/P + φ) ,

y(t) = −{B_F∙F_s/(−B_K∙K_s)}^½∙a∙cos(2∙π∙t/P + φ) .

Die Funktion y(t) ist gegenüber x(t) um π/2 phasenverschoben. Der Punkt x(t),y(t) - und damit auch K(t),F(t) - umläuft den stationären Punkt mit zunehmender Zeit t auf einer Ellipse im mathematisch positiven Sinn. Man beachte, daß die Periode P und die Ellipsenform nur im Grenzfall |K−K_s| ≪ K_s und |F−F_s| ≪ F_s gelten.

zurück

8. Interaktive graphische Lösung des DGL‑Systems

(muß noch mit Leben erfüllt werden)

zurück