4 Eindimensionale Differentialgleichungsmodelle

4.1 Exponentielles Wachstum

Gehen wir noch einmal zurück zu dem diskreten Modell des exponentiellen Wachstums (Gleichung 3.1) aus Abschnitt 3.2:

\[x_{n+1}=x_{n} + r \cdot x_{n}\]

Dieses Modell berechnet zu einem Anfangswert \(x_0\) die Bakterienzahl von Zeitschritt zu Zeitschritt (\(n \rightarrow n+1\)).

Betrachtet man nun die Bakterienzahl zu einer Zeit t und wählt den Zeitschritt \(\Delta t\), so bleibt mit

\[x(t+\Delta t)=x(t) + r \cdot x(t) \cdot \Delta t\; ,\]

für \(\Delta t=1\) alles beim Alten. Nun stellen wir die Gleichung um und lassen \(\Delta t\) immer kleiner werden. Der Differenzenquotient geht über in den Differentialquotienten

\[\begin{array}{ccc} \Delta x = \frac{x(t+\Delta t)-x(t)}{\Delta t} &=& r \cdot x(t)\\ \downarrow &&\\ \frac{dx(t)}{dt} &=& r \cdot x(t) \end{array}\]

Man erhält die gewöhnliche Differentialgleichung (Differentialgleichung) des exponentiellen Wachstums:

\[ \frac{dx(t)}{dt} = r \cdot x(t) \quad \mbox{kurz} \quad x'=r\cdot x \quad \mbox{oder} \quad \dot x =r\cdot x \, , \tag{4.1}\]

sie hat die Lösung

\[ x(t)=C \cdot e^{ r \cdot t} \; . \tag{4.2}\]

C ist eine freie Konstante, die durch den Anfangswert bestimmt wird. Der Parameter \(r\) heißt Wachstumsrate.

Ist ein Anfangswert \(x(t_0)=x_0\) vorgegeben, so nennt man die Suche nach der Lösung der Differentialgleichung, die zu diesem Anfangswert gehört, ein Anfangswertproblem (AWP):

Man kann leicht überprüfen, dass 4.2 eine Lösung von Gleichung 4.1 ist. Hierzu leitet man 4.2 nach t ab:

\[ \frac{dx(t)}{dt} = r \cdot C\cdot e^{ r \cdot t} = r \cdot x(t) \; . \]

Für alle, die noch keine Vorlesung zu Differentialgleichungen gehört haben, sei an dieser Stelle der Anhang C empfohlen.

4.1.1 Vergleich diskretes Modell und Differentialgleichung

Betrachten wir noch einmal den diskreten Fall

\[x_{n+1} = (1+r) \cdot x_n = x_0 \cdot (1+r)^{n+1} = x_0 \cdot e^{(n+1)\cdot \ln(1+r)}\]

Für n=0 ergibt sich \(x_1= x_0 \cdot e^{\ln(1+r)}\). Im stetigen Fall gilt \(x(1)= x_0 \cdot e^{a \cdot 1}\). Um also mit dem diskreten und dem stetigen Modell dasselbe Wachstum zu erreichen, muss \(a=ln(1+r)\) gesetzt werden.

\[x(t)=x_0 \cdot e^{(\ln (1+r)) \cdot t} = x_0 \cdot e^{\ln (1+r)^t}= x_0 \cdot (1+r)^t \; ,\]

Die Wachstumsrate \(a\) ist also immer kleiner als \(r\). Warum ist dies so? Betrachtet man das Modell diskret, so wird immer für einen ganzen Zeitschritt \(\Delta t\) die Änderung auf der Basis der Bakterienzahl zu Beginn des Zeitschritts \(t_0\) angenommen:

\[ \Delta x = r \cdot x(t_0) \cdot \Delta t \]

Im Intervall \([t, t+\Delta t]\) wird also \(x\) als konstant angenommen. Im kontinuierlichen Fall nimmt die Zahl der Bakterien aber stetig zu (Abbildung 4.1).

Abbildung 4.1: Vergleich des exponentiellem Wachstums für das diskrete Modell mit der Wachstumsrate \(r=1\) (grün) und das Differentialgleichungs-Modell mit der Wachstumsrate \(a=1\) (blau) und \(a = \ln(1+r)\) (gelb).

Will man also ganz speziell im Differentialgleichungs-Modell nach einem Zeitschritt eine Verdopplung (r=1) erreichen, so gilt \(a =\ln 2\approx 0.69\):

\[x(t)=x_0 \cdot e^{(\ln 2) \cdot t} = x_0 \cdot e^{\ln 2^t}= x_0 \cdot 2^t \; ,\]

Umgekehrt gilt für \(a =1\), dass man nach einem Zeitschritt eine ver-e-fachung der ursprünglichen Größe erhält:

\[x(t)=x_0 \cdot e^{1 \cdot t} \; .\]

4.1.2 Zur numerischen Lösung der Differentialgleichung

Betrachten wir noch einmal das AWP

\[\dot x(t) = r\cdot x(t), \quad x(0)=1\]

dessen Lösung \(x(t)=e^{r\cdot t }\) ist. Für \(r=1\) gilt also \(x(t)=e^t\) und somit \(x(10)=e^{10}\approx 22026\).

Die Lösung kann numerisch durch das Euler-Verfahren (Anhang D.1) approximiert werden:

\[x_{k+1} = x_k + r\cdot x_k\cdot \Delta t\qquad k=0,1,..\]

\(\Delta t\) ist die Schrittweite des Verfahrens. Je kleiner \(\Delta t\) desto genauer das Ergebnis.

Das Euler-Verfahren liefert für \(\Delta t=0.1\) als Lösung x(10)=13781. Dies ist sehr weit von der tatsächlichen Lösung entfernt. Verkleinert man den Zeitschritt, wird das Ergebnis besser. Für \(\Delta t=0.001\) liefert es das Ergebnis x(10)= 21917. Man kann zwar den Zeitschritt weiter verkleinern, dies verlängert die Rechenzeit erheblich und das Ergebnis wird dadurch nicht unbedingt besser. Das Euler-Verfahren ist nicht geeignet, um numerisch exakte Ergebnisse zu bekommen.

Setzen wir nun den Zeitschritt auf \(\Delta t=1\), dann liefert das Verfahren x(10)=1024. Dies ist genau die Lösung, die uns das diskrete Modell

\[x_{n+1} = x_n + r\cdot x_n \qquad n=0,1,..\]

für \(r=1\) und \(y_0=1\) liefert.

Es ist wichtig zu verstehen, dass das diskrete Modell und die Differentialgleichung zwei unterschiedliche Sachverhalte beschreiben. Beim diskreten Modell geht man davon aus, dass die Werte von Schritt zu Schritt berechnet werden. Es wird also keine Lösung approximiert, die diskreten Werte sind die Lösung. Bei der numerischen Lösung der Differentialgleichung wird die Lösungsfunktion durch das numerische Verfahren approximiert. Das einfachste dieser Verfahren, das Euler-Verfahren, ist gerade die Approximation durch das diskrete Modell. Dieses Verfahren macht aber bei der Lösung große Fehler. Die Approximation an die tatsächliche Lösung der Differentialgleichung gelingt umso besser, je besser das Verfahren ist bzw. je kleiner der Zeitschritt ist. Letzteres funktioniert aber nur in begrenztem Maße. Eine Implementierung in Matlab ist in Abschnitt D.1.2 angegeben.

4.2 Logistisches Wachstum

In Abschnitt 3.2 haben wir das logistische Wachstum kennen gelernt. Im diskreten Fall hatte das logistische Wachstum folgende Form

\[x_{n+1} =x_n + r \cdot \left (1-\frac{x_n}{K} \right )\cdot x_n \; .\]

Diese Gleichung geht nun in

\[ \frac{dx(t)}{dt}=r \cdot \left (1-\frac{x}{K} \right ) \cdot x \] über.

Die Lösungen der logistischen Differentialgleichung weisen im Unterschied zur Differenzengleichung unabhängig von der Parametrisierung nun kein Chaos mehr auf.

Abbildung 4.2: Lösungen des logistischen Wachstums (r=1, K=1) für verschiedene Anfangswerte.

Man kann folgende Lösungstypen unterscheiden:

Lösungen, die zwischen 0 und K bleiben und monoton wachsen
Lösungen, die sich monoton fallend K annähern
Lösungen, die von kleiner 0 gegen \(-\infty\) streben
Lösungen, die konstant sind (\(x\equiv0\) und \(x\equiv\) K)

Das logistische Wachstum weist für den 1. Lösungstyp für kleine Anfangswerte \(x_0\) zunächst ein exponentielles Wachstum auf, denn in diesem Fall gilt \[\frac{dx(t)}{dt}=r \cdot \left (1-\frac{x}{K} \right ) \cdot x = r \cdot x-\frac{r}{K}\cdot x^2 \approx r \cdot x\, . \tag{4.3}\]

Für Startwerte in der Nähe von \(K/2\) ist das Wachstum zunächst linear.

Die sigmoiden Lösungen (\(0<x_0<K/2\)) der logistischen Differentialgleichung sind punktsymmetrisch am Wendepunkt. Daher nähern sich Lösungen in der Nähe von \(K\) exponentiell dem Wert \(K\).

Die logistische Differentialgleichung kann man noch analytisch lösen, die beschränkten Lösungen lauten

\[x(t) = K \cdot \frac{e^{r(t-m)}}{1+e^{r(t-m)}} \qquad -\infty<m<\infty\]

Dies ist eine Funktionenschar, der Schar-Parameter \(m\) kann alle reellen Werte annehmen und hängt für das Anfangswertproblem \(x(0)=x_0\) von \(x_0\) ab.

Weitere Informationen zu Differentialgleichungen findet man in Anhang C .

4.2.1 Beispiel Bohnenwachstum

Abbildung 4.3 zeigt, dass man mit dem kontinuierlichen Modell zum logistischen Wachstum nach Gleichung 4.3 die Daten des Bohnenwachstums aus Abschnitt 2.2.1 einigermaßen reproduzieren lassen. Eine gute Anpassung an die Daten der Hauptwachstumsphase erhält man allerdings nur, wenn man einen Anfangswert annimmt, der deutlich größer als der gemessene ist. Die kann daran liegen, dass die Anfangs sehr kleinen Werte einen Messfehler aufweisen. Wahrscheinlicher ist aber, dass sich das reale Wachstum nicht exakt mit einem so einfachen Modell beschreiben lässt, weil es durch viele weitere Faktoren wie z.B. Temperatur und Lichtintensität beeinflusst ist.

Abbildung 4.3: Simulation des Modells zum logistischen Wachstum nach Gleichung 4.3 mit der Wachstumsrate r, der Kapazität K und dem Anfangswert H(0) (schwarz) verglichen mit den Daten aus Abschnitt 2.2.1 (grün).

4.3 Richtungsfeld

Kann man einer Differentialgleichung ihren Lösungsverlauf “ansehen”?

Betrachten wir zum Beispiel \(\dot{x}=a \cdot x\). Die Steigung (Ableitung) der Funktion \(x\) an der Stelle \(t\) ist proportional zum Funktionswert an der Stelle \(t\).

Um dies zu visualisieren betrachtet man ein \(t,x\)-Koordinatensystem. Zu beliebig vielen Paaren \((t,x)\)-Paaren zeichnet man nun die Steigung an diesem Punkt ein, indem man an \((t,x)\) ein Linienelement oder einen Vektor mit der Steigung \(\dot x=f(t,x)\) einzeichnet (siehe Abbildung 4.4). Ohne die Lösungen der Differentialgleichung zu kennen, weiß man, dass eine Lösung, die durch einen Punkt \((t,x)\) geht, das Linienelement als Tangente haben muss. Zeichnet man viele Linienelemente, so kann man am Richtungsfeld den Verlauf der Lösungen “sehen”.

Die kleinen Striche im Richtungsfeld haben keinen bestimmten Namen, wir nennen sie hier mal Linienelemente. Sie werden von den Lösungskurven niemals durchkreuzt. Die Lösungsbahnen verlaufen immer tangential zu den Linienelementen. Verschiedene Lösungen können sich also nicht schneiden! Zu jedem Anfangswert \(x(t_0)=x_0\) gibt es meistens genau eine Lösung¹. Ist die Differentialgleichung autonom, d.h. die rechte Seite hängt nicht explizit von der Zeit ab, ist das Richtungsfeld t-translationsinvariant. Dies bedeutet, dass sich die Steigungen in t-Richtung nicht verändern.

4.4 Fixpunkte (Stationäre Zustände, Gleichgewichtszustände)

Autonome Systeme² der Form \(\dot x=f(x)\) haben häufig spezielle Lösungen der Form

\[x(t)=x^\ast=\text{konstant} \,,\]

d.h. dass die Änderung von \(x\) gleich Null ist, also \(\dot x=0\). Die Konstante \(x^\ast\) erfüllt also die Gleichung

\[f(x^\ast)=0\,.\]

Diese zeitlich unveränderlichen Lösungen heißen stationäre Zustände, stationäre Lösungen, Gleichgewichtszustände, Gleichgewichtslösungen oder kurz Fixpunkte.

Im Fall des logistischen Wachstums gibt es, wie oben erwähnt, die zwei stationären Lösungen \(x^\ast=0\) und \(x^\ast=K\).

4.4.1 Stabilität von Fixpunkten

Es gibt zwei Arten von Fixpunkten

stabil: kleine Abweichungen wachsen nicht, genauer:
- asymptotisch stabil: kleine Abweichungen verschwinden
- marginal stabil: kleine Abweichungen wachsen nicht
instabil: es gibt kleine Abweichungen, die wachsen

Für das logistische Wachstum gilt, dass \(x^\ast=0\) instabil ist. Sobald man vom Fixpunkt auch nur ein wenig abweicht, fängt die Lösung an zu wachsen. \(x^\ast=K\) ist stabil. Sobald man ein wenig von \(K\) abweicht, läuft das System zu \(K\) zurück.

Mathematisch präzise

Ein Fixpunkt heißt stabil, wenn es zu jedem \(\varepsilon > 0\) ein \(\delta > 0\) gibt, so dass alle Lösungen, die anfangs höchstens um \(\delta\) vom Fixpunkt abweichen, für alle späteren Zeiten höchstens um \(\varepsilon\) vom Fixpunkt abweichen.

Diese Definition schließt die “marginal” stabilen Fälle ein. Ein Fixpunkt heißt auch dann “instabil”, wenn nicht alle kleinen “Störungen” zu wachsenden Abweichungen führen. Es muss nur solche wachsenden Störungen geben.

4.4.2 Bestimmung der Stabilität

Bei autonomen Systemen mit einer Zustandsvariablen kann schnell festgestellt werden, ob ein Fixpunkt stabil ist oder nicht. Man schaut sich hierzu die Änderung an (Abbildung 4.5): Betrachtet man den Fixpunkt \(x^\ast=0\), so sieht man, dass die Änderung positiv wird, sobald man nur ein wenig von \(x^\ast\) nach rechts abweicht. Damit wird \(x\) immer weiter erhöht. Man sagt der Fixpunkt \(x^\ast=0\) ist instabil, weil kleine Störungen dazu führen, dass die Lösung vom Fixpunkt weg läuft. Anders ist es bei \(x^\ast=K\). Verringert man \(x\) ein wenig, so ist die Änderung positiv, \(x\) nimmt also wieder zu. Damit kommt das System von \(x^\ast=K\) nicht weg, der Fixpunkt ist stabil.

Abbildung 4.5: Differentialgleichung des logistischen Wachstums. Aufgetragen ist die rechte Seite der Differentialgleichung \(f(x)=r\cdot x\cdot(1-\frac{x}{K})\) über x. Für \(f>0\) nimmt x zu, für \(f<0\) nimmt x ab und für \(f=0\) bleibt x unverändert, da f gerade die Änderung von x angibt (rechte Seite der Differentialgleichung). Die Pfeile geben die Änderung von x an.

Analytisch kann man die Stabilität also bestimmen, indem man die Steigung der rechten Seite der Differentialgleichung bezogen auf x untersucht. Steigt die Funktion f, so ist der Fixpunkt instabil, fällt sie ist er stabil. Man bildet hierzu also die Ableitung \(\frac{df(x^\ast)}{dx}\) der rechten Seite der Differentialgleichung nach der Zustandsvariablen x im stationären Zustand³. Es gilt :

\(\frac{df(x^\ast)}{dx} > 0\),	dann ist \(x^\ast\) ist instabil
\(\frac{df(x^\ast)}{dx}< 0\),	dann ist \(x^\ast\) ist stabil
\(\frac{df(x^\ast)}{dx}= 0\),	es sind weitere Analysen erforderlich

Für das logistische Wachstum erhält man nun folgende Stabilitäten

\[\dot x= f(x)=r \cdot \left (1-\frac{x}{K} \right ) \cdot x \qquad r> 0\]

gilt \[\frac{df(x)}{dx} =r -\frac{2r}{K} \cdot x\] Für \(x^\ast=0\) gilt: \(\frac{df(x^\ast)}{dx}=r>0\), der Zustand ist instabil
Für \(x^\ast=K\) gilt: \(\frac{df(x^\ast)}{dx} =r-2r<0\), der Zustand ist stabil

Was im Detail passiert

Man betrachte \(x\) als Abweichung vom Fixpunkt \(x^\ast\) mit der Störung \(z\):

\[ x(t)=x^\ast+z(t) \]

Für \(x\) gilt

\[\dot x(t) = \frac{dx(t)}{dt}=\frac{d(x^\ast+z(t))}{dt}=\frac{dx^\ast}{dt}+\frac{dz(t)}{dt}=\frac{dz(t)}{dt} =\dot z(t)\]

da \(\frac{dx^\ast}{dt}=0\), weil \(x^\ast\) konstant ist.
Andererseits gilt die Taylorentwicklung (siehe Anhang A.12 )

\[\dot x(t) = f(x(t))=f(x^\ast+z(t))=f(x^\ast)+\frac{df(x^\ast)}{dx}\cdot z(t) + \cdots\]

Betrachtet man nun die Taylorentwicklung nur bis zum ersten Glied⁴, so erhält man eine neue, lineare Differentialgleichung für die Störung \(z=x^\ast-x\).

\[ \dot z(t)=f(x^\ast)+\frac{df(x^\ast)}{dx}\cdot z(t) = \frac{df(x^\ast)}{dx}\cdot z(t) \tag{4.4}\]

mit der Lösung

\[z(t)= z_0\cdot e^{\frac{df(x^\ast)}{dx}\cdot t}\]

Damit wächst die “Störung” für \(\frac{df(x^\ast)}{dx}>0\) exponentiell (instabil) und zerfällt exponentiell für \(\frac{df(x^\ast)}{dx}<0\) (stabil). Dieses Verfahren nennt man auch “Linearisierung am stationären Zustand”. Man sagt, man habe die ursprüngliche Differentialgleichung linearisiert. Nun ist auch klar, warum man für \(\frac{df(x^\ast)}{dx}=0\) keine Aussage über die Stabilität des Fixpunkts der ursprünglichen Gleichung treffen kann. Da die Störung unverändert bleibt, hängt die Stabilität nun womöglich von den höheren Termen der Taylorentwicklung ab, die beim linearisierten System vernachlässigt wurden.

4.4.3 Beispiel: Schafe am Deich

Es soll bestimmt werden, wie viele Schafe einen Deich beweiden können. Es muss sichergestellt werden, dass das Gras nicht vollständig abgefressen wird, um die Stabilität des Deiches nicht zu gefährden. Hierzu wird nun ein sehr einfaches konzeptionelles Modell aufgestellt. Dieses kann ganz sicher noch nicht das konkrete Problem lösen, zeigt aber das Herangehen an eine solche Fragestellung.

Die Änderung der Vegetation wird durch ihr Wachstum und ihre Beweidung bestimmt: Änderung der Vegetation = Wachstum - Beweidung Durch nicht weiter berücksichtigte Umweltfaktoren wird das Gras (Vegetation V) auch ohne Beweidung nur bis zu einer Kapazität K wachsen. Ohne Beweidung ist also die Annahme des logistischen Wachstums sinnvoll:

\[\frac{dV(t)}{dt}=r \cdot \left (1-\frac{V(t)}{K} \right ) \cdot V(t)\]

Nimmt man nun eine Beweidung durch Herbivoren an, so hängt die gefressene Menge Gras sowohl von der Anzahl der Herbivoren H als auch von der zur Verfügung stehenden Grasmenge V ab. Je mehr Schafe auf dem Deich, desto mehr wird gefressen. Nicht jedoch gilt, dass je mehr Gras vorhanden ist umso mehr gefressen wird, da die Schafe nur eine begrenzte Aufnahmekapazität haben – sie werden satt bzw. benötigen eine gewisse Zeit zur Beweidung. Daher ist der Ansatz mit einem Sättigungsterm sinnvoll. Dieser wird später noch ausführlicher diskutiert.

Hier gehen wir davon aus, dass die maximale Fressrate b erst bei sehr großer Nahrungsmenge erreicht wird und die Fressrate bei einer Grasmenge \(V_H\) gerade zur Hälfte erreicht wird. Deshalb heißt dieser Parameter auch Halbsättigungskonstante. Die gesamte Gleichung für die Vegetation V lautet nun (die Zeitabhängigkeit von V wird hier nun nicht mehr mit aufgeschrieben):

\[\frac{dV}{dt}= r \cdot V \cdot (1-\frac{V}{K}) - b \cdot \frac{V}{V+V_H} \cdot H\]

Hierbei ist es nun wichtig zwischen der Zustandsvariablen V genauer V(t) und den Parametern r, K, b, \(V_H\) und H zu unterscheiden. Parameter sind Größen deren Werte wir vorgeben. Die Zustandsvariable ist unsere gesuchte Funktion, deren zeitlichen Verlauf wir suchen. Insbesondere interessiert uns, wann die Zustandsvariable konstant ist, d.h. keine Änderung erfährt. Wir suchen also die Fixpunkte.

Um das Ganze nun zu vereinfachen und so an dieser Stelle viel Rechnerei zu vermeiden setzen wir r=b=\(V_H\)=1 und K=10. Dies dient hier der Übersichtlichkeit der Rechnungen. In einem realistischen Modell ist es natürlich komplizierter. Insbesondere haben die Parameter alle unterschiedliche Einheiten und können zwar denselben Zahlenwert haben, sind damit aber noch lange nicht gleich.

Wir betrachten nun also das vereinfachte Modell

\[ \frac{dV}{dt}= V \cdot (1-\frac{V}{10}) -\frac{V}{V+1} \cdot H \tag{4.5}\]

Gesucht sind nun die Fixpunkte von 4.5 in Abhängigkeit von H. Wir haben zwar bisher von Schafen gesprochen, dies bedeutet aber nicht, das H nur ganzzahlige Werte annehmen kann. Dies versteht man, wenn man H eher als Herbivorenbiomasse auffasst.

Die Fixpunkte von 4.5 erhält man, wenn man die rechte Seite gleich Null setzt:

\[ V \cdot (1-\frac{V}{10}) -\frac{V}{V+1} \cdot H=0 \tag{4.6}\]

Dies ist für V=0 erfüllt. Damit hat man den ersten Fixpunkt \(V^\ast_1 = 0\) gefunden.

Es gilt aber auch für

\[\begin{aligned} (1-\frac{V}{10}) -\frac{1}{V+1} \cdot H&=0 \Longleftrightarrow\\ V+1 -\frac{1}{10} V^2 - \frac{V}{10} - H&=0 \Longleftrightarrow\\ -V^2 +9V+10 -10 H &=0 \Longleftrightarrow\\ V^2-9V-10+10 H & =0 \end{aligned}\]

Löst man diese quadratische Gleichung so erhält man zwei weitere Lösungen

\[ \nonumber V^\ast_{2/3} = \frac{9}{2} \pm \sqrt{ \frac{9^2}{4} - 10 H+ 10} = 4,5 \pm \sqrt{30,25 - 10 H} \]

Ist der Term unter der Wurzel positiv (\(1\leq H \leq 3,025\)), so existieren also zwei weitere Fixpunkte. Diese sind für uns aber nur interessant, wenn sie nicht negativ sind, das es sich ja um Vegetation handelt. Der Fixpunkt \(V^\ast_{2}\) ist immer positiv, wenn er existiert, \(V^\ast_{3}\) nur dann , wenn \(4,5 - \sqrt{30,25 - 10 H} >0\) gilt, also für \(H>1\).

Die Lage der Fixpunkte ist in Abbildung 4.6 dargestellt.

Abbildung 4.6: Lage der Fixpunkte der Vegetation (Gras am Deich) in Abhängigkeit von der Herbivorenbiomasse (Schafe)

Man sieht, dass es für kleine “Schafzahlen” (\(H<1\)) nur zwei nicht-negative Fixpunkte gibt. Entweder die Vegetation hat sehr hohe Werte oder ist Null. Wenn das Modell sinnvoll ist, sollte sich bei geringer Schafzahl natürlich der hohe Wert als stabil erweisen. Für größere Schafzahlen existieren bis zu dem kritischen Wert von \(H_{crit}=3,025\), bei dem die Wurzeläste zusammenlaufen, drei Fixpunkte. Für noch größere Schafzahlen gibt es als stationäre Lösung nur noch “kein Gras mehr” als Lösung (\(V^\ast_1=0\)). Der Deich ist also überweidet, sofern sich der Fixpunkt \(V^\ast_1=0\) als stabil erweist, also alle Lösungen dorthin streben. Dies würde bedeuten, dass sobald man eine größere Menge Schafe auf den Deich bringt, das Grass komplett abgefressen wird. Für kleinere Schafzahlen (\(1<H<H_{crit}\)) könnte sich das System auf einen der beiden Wurzeläste zu bewegen, wenn einer der beiden stabil ist. Dies gilt es zu untersuchen.

Hierzu wird zunächst gemäß 4.4.2 \(\frac{df(V)}{dV}\) bestimmt:

\[ \frac{df(V)}{dV} = 1-\frac{2V}{10} -\frac{1}{(V+1)^2}\cdot H \tag{4.7}\]

Die Bestimmung der Stabilität der Fixpunkte erfolgt nun wieder durch Einsetzen. Zunächst wird \(V^\ast_1=0\) in Gleichung 4.7 eingesetzt:

\[\frac{dfV^\ast_1=0}{dV} = 1-\frac{2\cdot 0 }{10} -\frac{1}{(0+1)^2}\cdot H = 1-H\]

\(V^\ast_1\) ist also stabil für \(1-H <0\) also für \(H>1\) und instabil für \(1-H >0\) also für \(H<1\) genau wie wir es erwartet haben. Die Bestimmung der Stabilität von \(V^\ast_2\) und \(V^\ast_3\) erfolgt genauso, ist aber deutlich rechenintensiver und wird hier mal ausgelassen. Man kann dies durchaus noch von Hand rechnen. Einfacher ist es mit Maple oder der Symbolic-Toolbox von Matlab. Das Ergebnis ist in Abbildung 4.7 dargestellt.

Man erkennt in Abbildung 4.7, dass bei Werten von H unter 1 (wenig Schafe), die Vegetation startend von einem beliebigen Anfangswert zum stabilen Fixpunkt \(V^\ast_2\) strebt. Im Bereich \(1<H<H_{crit}\) wird die Vegetation wachsen, sofern der Startwert oberhalb von \(V^\ast_3\) liegt. Ist der Anfangswert unterhalb von \(V^\ast_3\), stirbt die Vegetation aus und der Deich ist gefährdet. Oberhalb von \(H_{crit}\) sind zu viele Schafe auf dem Deich. Egal, wie viel Vegetation da ist, letztlich wird alles weggefressen, der einzige stabile Zustand \(V^\ast_1=0\) wird angestrebt. Der zeitliche Verlauf wird deutlich, wenn man sich die einzelnen Lösungskurven für verschiedene Werte von H im jeweiligen Richtungsfeld anschaut (siehe Abbildung 4.8).

Dies gilt nicht für alle Differentialgleichungen. Weiter hilft hier der Satz von Picard-Lindelöf: f muss Lipschitz-Bedingung für y erfüllen.↩︎
Systeme, die nicht explizit von der Zeit abhängen.↩︎
Achtung: Hier wird nun f(x) nach x abgeleitet. Die Funktion f selbst beschreibt die Änderung von x nach der Zeit.↩︎
\[f(x_0+\Delta t)\approx f(x_0)+f'(x_0)\cdot \Delta t\]↩︎