Anhang I: Stabilitätsanalyse autonomer Systeme

In diesem Anhang wird die Stabilitätsanalyse zweidimensionaler Systeme vertieft. Der Anhang beinhaltet eine Weiterführung von Kapitel 5.1 und 13.1.1. Für die Analyse von Systemen mit mehr als zwei gekoppelten Differentialgleichung sei auf die Literatur zu dynamischen Systemen, z.B. (Strogatz 2015) oder (Jetschke 1989) verwiesen.

I.1 Fixpunkte

Die Fixpunkte autonomer Systeme mit zwei Zustandsvariablen der Form

\[\begin{aligned} \dot X &= f(X,Y)\\[-4ex]\nonumber \dot Y &= g(X,Y)\nonumber \end{aligned} \tag{I.1}\]

lassen sich bestimmen, indem \(f(X,Y) =0\) und \(g(X,Y)=0\) gesetzt wird, und die Lösungen \((X^\ast,Y^\ast)\) dieses Gleichungssystems bestimmt werden.

I.2 Geschlossene Bahnen

Im Gegensatz zu eindimensionalen Systemen können im zweidimensionalen Fall geschlossene Lösungsbahnen auftreten. Im Normalfall ist die geschlossene Bahn das Bild einer periodischen Lösung. Man nennt eine solche periodische Bahn stabil, wenn alle Bahnen, die in ihrer Nähe starten, zu dieser hin laufen. Man nennt solche stabilen Bahnen auch Grenzkreise. Periodische Bahnen können auch instabil sein. Es gilt der

Satz von Poincare

Bei autonomen Systemen mit zwei Zustandsvariablen läuft jede Lösung entweder

auf einen stationären Zustand zu, oder
nähert sich einer geschlossenen Kurve (Grenzkreis) an, oder
läuft nach unendlich

I.2.1 Teilgleichgewichte (Isoklinen)

Für das System I.1 sind die Teilgleichgewichte durch die Mengen aller Punkte, für die \(\dot X=0\) bzw. \(\dot Y=0\) gilt, gegeben.

\[\mbox{Teilgleichgewicht von X} = \lbrace \dot X=0\rbrace= \lbrace(X,Y) \vert f(X,Y)=0\rbrace\]

\[\mbox{Teilgleichgewicht von Y} = \lbrace \dot Y=0 \rbrace= \lbrace(X,Y) \vert \; g(X,Y)=0\rbrace\]

I.3 Stabilität der Fixpunkte

Ein Fixpunkt heißt stabil, wenn er die Lösungen aus seiner Nachbarschaft anzieht. Man nennt ihn aber auch dann noch stabil (marginal stabil), wenn sich die Lösungen zu allen Anfangspunkten in seiner Nachbarschaft nicht entfernen, also in der Nähe bleiben, ohne zu ihm hinzustreben.

Die Stabilität eines Fixpunkts kann wie im eindimensionalen Fall untersucht werden. Man betrachtet hierzu das System als Vektor

\[ {\dot X \choose \dot Y} ={f(X,Y) \choose g(X,Y)} \]

Analog zum eindimensionalen Fall entwickelt man die rechte Seite des Differentialgleichungssystems am Fixpunkt \((X^\ast,Y^\ast)\) nach Taylor bis zum ersten Glied und erhält so das linearisierte System (ist für jeden Fixpunkt unterschiedlich!!):

\[ {\dot x \choose \dot y} = J(X^\ast,Y^\ast) \cdot {x \choose y} \tag{I.2}\]

mit

\[ J(X^\ast,Y^\ast) =\left ( \begin{array}{cc} \frac{\displaystyle\partial f(X^\ast,Y^\ast)}{\displaystyle\partial X} & \frac{\displaystyle\partial f(X^\ast,Y^\ast)}{\displaystyle\partial Y}\\ &\\ \frac{\displaystyle\partial g(X^\ast,Y^\ast)}{\displaystyle\partial X} & \frac{\displaystyle\partial g(X^\ast,Y^\ast)}{\displaystyle\partial Y} \end{array} \right ) \tag{I.3}\]

Die Matrix \(J:=J(X^\ast,Y^\ast)\) (I.3) ist die Jacobi-Matrix des Systems. Die Komponenten der Jacobi-Matrix sind die partiellen Ableitungen von f und g nach X und Y jeweils an der Stelle\((X^\ast,Y^\ast)\).

Das linearisierte System I.2 ersetzt in der Nähe des Fixpunkts \((X^\ast,Y^\ast)\) das kompliziertere nichtlineare System I.1. Das linearisierte System ist am Fixpunkt \((X^\ast,Y^\ast)\) entwickelt und gibt die Abweichung (Störung) von diesem an. Die Abweichung besteht darin, dass der “Rest”, also die Summanden höherer Ordnung in der Taylorreihe, weggelassen wird. Die Idee dabei ist, dass man davon ausgeht, dass sich das Modell in der unmittelbaren Nähe des Fixpunkts nahezu linear verhält. Man untersucht nun also stellvertretend das lineare Modell, das seinen Fixpunkt bei (0,0) hat. Lineare Differentialgleichungssysteme sind analytisch lösbar. Hat die Matrix \(J\) zwei verschiedene Eigenwerte \(\lambda_1\) und \(\lambda_2\), so kann man eine allgemeine Lösung des linearen Systems I.2 angeben:

\[\vec{z(t)} = C_1 \cdot e^{\lambda_1 t} \cdot \vec{v_1} + C_2 \cdot e^{\lambda_2 t} \cdot \vec{v_2}\]

Hierbei ist \(\vec{v_1}\) Eigenvektor zum Eigenwert \(\lambda_1\) und \(\vec{v_2}\) Eigenvektor zum Eigenwert \(\lambda_2\), \(C_1\) und \(C_2\) sind freie Konstanten.

Um das Lösungsverhalten und damit die Stabilität des Fixpunkts zu bestimmen, reicht es aus, die Eigenwerte \(\lambda\) der Matrix \(J\) zu kennen und zu interpretieren. Sind beide Eigenwerte reell, so handelt es sich bei den Summanden entweder um exponentielles Wachstum oder exponentiellen Zerfall. Ist also einer der Eigenwerte größer Null, wächst die Störung, der Fixpunkt ist instabil. Sind beide Eigenwerte negativ, so verschwindet eine anfängliche Störung, der Fixpunkt ist stabil.

Für einen komplexen Eigenwert \(\lambda=a+ ib\) gilt

\[e^{\lambda t} = e^{(a+ ib) t} = e^{at} \cdot e^{ibt} = e^{at} \cdot (\cos bt + i \sin bt)\]

Das System vollzieht also Schwingungen. Ist der Realteil \(a\) kleiner Null, so verschwindet die Störung, ist der Realteil \(a\) größer Null, so explodieren Störungen. Bei \(a=0\) bleiben Störungen erhalten. Über die Stabilität des Ursprungssystems kann dann nicht entschieden werden.

Bezüglich der Stabilität der Fixpunkte des Systems I.1 gilt der folgende Satz:

Stabilität der Fixpunkte autonomer Systeme

Der Fixpunkt \((X^\ast,Y^\ast)\) des Systems I.1 ist asymptotisch stabil, falls alle Eigenwerte der Jacobi-Matrix I.3 negative Realteile besitzen.
Der Fixpunkt \((X^\ast,Y^\ast)\) ist instabil, falls mindestens ein Eigenwert der Jacobi-Matrix I.3 einen positiven Realteil besitzt.
Sind alle Realteile aller Eigenwerte der Jacobi-Matrix I.3 kleiner oder gleich Null und mindestens ein Realteil gleich Null, so kann nicht über die Stabilität des Fixpunkts \((X^\ast,Y^\ast)\) entschieden werden.

Sind beide Eigenwerte des linearisierten Systems reell und verschieden, so gibt es zwei Eigenrichtungen. Diese Richtungen sind selbst Lösungen und können von anderen Lösungen nicht überquert werden. Entlang der Eigenrichtungen bewegen sich die Lösungen in gerader Linie auf den Fixpunkt hin (\(\lambda<0\)) oder weg (\(\lambda>0\)). Die absolute Größe des Eigenwerts bestimmt dabei die Geschwindigkeit. Für Lösungen, die nicht auf einer Eigenrichtung starten, gibt es Einflugschneisen. Für z.B. \(\lambda_1=-0.1\) und \(\lambda_2=-1\) bewegt sich das System zuerst schnell in Richtung von \(v_2\), und dann langsam in Richtung von \(v_1\).

Sind die Eigenwerte komplex (genauer: nicht reell) gibt es keine reellen Eigenrichtungen. Das System bewegt sich an keiner Stelle radial auf den Fixpunkt zu.

Hat das linearisierte System zwei verschiedene negative, reelle Eigenwerte besitzt, so heißt der Fixpunkt stabiler Knoten, sind die Eigenwerte beide positiv heißt der Fixpunkt instabiler Knoten. Sind die Eigenwerte komplex und haben beide negative Realteile, so heißt der Fixpunkt stabiler Fokus, sind die Realteile beide positiv, so heißt der Fixpunkt instabiler Fokus. Ein Fixpunkt, an dem das linearisierte System einen positiven und einen negativen Eigenwert, bzw. einen positiven und einen negativen Realteil besitzt, heißt Sattel. Die sind die 5 Hauptfälle. Bei zusammenfallenden Eigenwerten und bei Eigenwerten, die Null werden, gibt es weitere Formen des Systemverhaltens (siehe Abb. I.1 und Tab. I.1), auf die hier nicht näher eingegangen wird.

Bemerkung: Die Bezeichnungen der Fixpunkte sind in der Literatur leider nicht eindeutig.

Um das Systemverhalten in der Nähe des betrachteten Fixpunkts zu bestimmen, müssen also nur die Eigenwerte der Matrix \(J\) berechnet werden. Die Eigenwerte \(\lambda\) sind die Lösungen der quadratischen Gleichung \(\det\left(J-{\lambda I}\right)=0\) (siehe Anhang F). Man kann leicht ausrechnen, dass die Lösungen nur von der Spur

\[ \tau = \frac{\displaystyle\partial f(X^\ast,Y^\ast)}{\displaystyle\partial X} + \frac{\displaystyle\partial g(X^\ast,Y^\ast)}{\displaystyle\partial Y} \]

und Determinante

\[ \Delta = \frac{ \displaystyle\partial f(X^\ast,Y^\ast)}{\displaystyle\partial X} \cdot \frac{\displaystyle\partial g(X^\ast,Y^\ast)}{\displaystyle\partial Y} - \frac{\displaystyle\partial f(X^\ast,Y^\ast)}{\displaystyle\partial Y} \cdot \frac{\displaystyle\partial g(X^\ast,Y^\ast)}{\displaystyle\partial X} \]

der Matrix J abhängen. Die Lösungen (Eigenwerte) lauten:

\[ \lambda_{1,2}=\frac{\tau}{2}\pm \sqrt{\frac{\tau ^2}{4}-\Delta } \tag{I.4}\]

Aus Spur \(\tau\) und Determinante \(\Delta\) kann der Typ des Fixpunkts ermittelt, ohne die Eigenwerte (I.4) zu berechnen. Darauf basiert auch die Darstellung in der \(\tau\)-\(\Delta\)-Ebene (Abbildung I.1). Die Darstellungen in Abbildung I.1 gelten nur in der Nähe der Fixpunkte. Weiter weg können sich die Bahnen infolge der Nichtlinearitäten des Systems “irgendwie” verhalten. Es ist also, als hätte man in Abbildung I.1 eine Lupe benutzt, um die Umgebung des Fixpunkts zu untersuchen.

Abbildung I.1: Stabilität der Fixpunkte autonomer Systemen mit zwei Variablen in Abhängigkeit von der Spur \(\tau\) und der Determinante \(\Delta\) der Jacobi-Matrix. Benannt sind die 5 Hauptfälle. Die Grenzfälle sind der stabile und der instabile Stern (rot), das Zentrum (orange) , der stabile und instabile degenerierte Knoten (hellblau) und die Grenzfälle, bei denen die Determinante verschwindet (braun).

In Tabelle I.1 sind alle Vorzeichenmöglichkeiten der Eigenwerte aufgelistet. Bei komplexen Eigenwerten unterscheidet man nach dem Vorzeichen des Realteils. Außerdem sind in der Tabelle die Sonderfälle Zentrum, Stern und Grenzfälle berücksichtigt. Man muss unterscheiden, ob die Matrix \(J\) ein Vielfaches der Einheitsmatrix \(I\) ist oder nicht.

Der jeweils größere von zwei reellen Eigenwerte wird als \(\lambda_1\) bezeichnet, der kleinere als \(\lambda_2\) . Es ist immer der größte Eigenwert, bei komplexen der mit dem größtem Realteil, der das Verhalten in der Nähe des stationären Zustandes dominant bestimmt!

Im nicht linearisierten Ausgangssystem (I.1) kann die Stabilität der marginal stabilen Fälle durch die bei der Linearisierung weggelassenen Terme höherer Ordnung beeinflusst werden. Besonders kritisch ist das, wenn die Jacobi-Matrix die Nullmatrix ist. Dann findet im linearisierten System keine Bewegung statt, das Systemverhalten des nichtlinearen Ausgangssystems wird nur durch die Nichtlinearitäten bestimmt.

\(A\neq aI\)	Determinante	Spur	Bezeichnung	Stabilität
\(\lambda_{1,2}\) komplex, \(Re(\lambda)<0\)	\(\Delta > \frac{\tau^2}{4}\)	\(\tau< 0\)	Fokus	stabil
\(\lambda_{1,2}\) imaginär, \(Re(\lambda)=0\)	\(\Delta > \frac{\tau^2}{4}\)	\(\tau=0\)	Zentrum	marginal stabil
\(\lambda_{1,2}\) komplex, \(Re(\lambda)>0\)	\(\Delta > \frac{\tau^2}{4}\)	\(\tau> 0\)	Fokus	instabil
\(\lambda_2<\lambda_1< 0\)	\(0<\Delta<\frac{\tau^2}{4}\)	\(\tau< 0\)	Knoten	stabil
\(\lambda_2<\lambda_1= 0\)	\(\Delta=0\)	\(\tau<0\)	Grenzfall Null	marginal stabil
\(\lambda_2< 0<\lambda_1\)	\(\Delta<0\)		Sattel	instabil
\(0=\lambda_2<\lambda_1\)	\(\Delta = 0\)	\(\tau>0\)	Grenzfall Null	instabil
\(0<\lambda_2<\lambda_1\)	\(0 <\Delta<\frac{\tau^2}{4}\)	\(\tau>0\)	Knoten	instabil
\(\lambda_2 = \lambda_1< 0\)	\(\Delta = \frac{\tau^2}{4}\)	\(\tau<0\)	degenerierter Knoten	stabil
\(\lambda_2 = \lambda_1=0\)	\(\Delta = 0\)	\(\tau = 0\)	lineare Bewegung	instabil
\(\lambda_2 = \lambda_1> 0\)	\(\Delta = \frac{\tau^2}{4}\)	\(\tau> 0\)	degenerierter Knoten	instabil
\(A=aI\)
\(\lambda_2=\lambda_1< 0\)	\(a < 0\)		Stern	stabil
\(\lambda_2=\lambda_1= 0\)	\(a = 0\)		ohne Bewegung	marginal stabil
\(\lambda_2=\lambda_1> 0\)	\(a > 0\)		Stern	instabil

Tabelle I.1: Phänomenologie der Fixpunkte autonomer Systemen mit zwei Variablen in Abhängigkeit von den Eigenwerten \(\lambda\) der Jacobi-Matrix \(J\).

I.4 Analyse des skalierten Rosenzweig MacArthur-Modell

Als Beispielsystem wählen wir das allgemeine Räuber-Beute-Modell (5.3) nach Rosenzweig und MacArthur aus Abschnitt 5.1 in der skalierten Form (H.1):

\[\begin{aligned} \dot B&= \left (1-\varepsilon B \right ) \cdot B\;-\; \frac{R}{1+\kappa B} \cdot B \\ \dot R &= \delta \frac{B}{1+\kappa B} \cdot R \; - \delta R \end{aligned}\]

Das System hat die Parameter \(\varepsilon\), \(\kappa\) und \(\delta\), für die angenommen wird, dass sie alle größer oder gleich Null sind. Für \(\varepsilon=0\) und \(\kappa=0\) hat man das skalierte Lotka-Volterra-Modell nach Gleichung 5.1. Das System hat 3 Fixpunkte

\(B^\ast\)=0 und \(R^\ast=0\) (trivial)
\(B^\ast=\frac{\displaystyle 1}{\displaystyle\varepsilon}\) und \(R^\ast=0\) (Beute ohne Räuber, bei \(\varepsilon>0\))
\(B^\ast=\frac{\displaystyle 1}{\displaystyle 1-\kappa}\) und \(R^\ast=\frac{\displaystyle 1-\kappa-\varepsilon}{\displaystyle(1-\kappa)^2}\) (positiver stationärer Zustand)

Der positive stationäre Zustand existiert im Zustandsraum (\(B,R\ge0\)), wenn \(\kappa+\varepsilon<1\) erfüllt ist. Für \(\kappa+\varepsilon=1\) fällt der Zustand (3) mit Zustand (2) zusammen. Die Jacobi-Matrix des System lautet

\[ J(B,R)= \left[ \begin {array}{cc} 1-2\,\varepsilon\,B-\dfrac {R}{\left( 1+\kappa\,B \right) ^{2}} &-\dfrac {B}{1+\kappa\,B}\\ \delta \dfrac {R}{ \left( 1+\kappa\,B \right) ^{2}} &\delta \left( \dfrac {B}{1+\kappa\,B}-1\right) \end {array} \right] \]

Stabilität von \((B_1^\ast,R_1^\ast)=(0,0)\)

Für \((B_1^\ast,R_1^\ast)=(0,0)\) erhält man

\[ J(B_1^\ast,R_1^\ast) = \left[ \begin {array}{cc} 1 &0\\ 0 &-\delta \end {array} \right] \]

Damit gilt \(\det J<0\), der Zustand ist instabil.

Stabilität von \((B_2^\ast,R_2^\ast)=(\frac{1}{\varepsilon},0)\) (Beute ohne Räuber)

Für \((B_2^\ast,R_2^\ast)=(\frac{1}{\varepsilon},0)\) erhält man

\[ J(B_2^\ast,R_2^\ast) = \left[ \begin {array}{cc} -1 &-\dfrac{1}{\kappa+\varepsilon}\\ 0 &\delta \left ( \dfrac{1}{\kappa+\varepsilon} -1\right ) \end {array} \right] \]

Es gilt \[\det J =\delta \left ( 1-\dfrac{1}{\kappa+\varepsilon} \right )\]

\(\kappa+\varepsilon <1\), Zustand 3 existiert:
\(\det J<0\), der Zustand 2 ist instabil.
\(\kappa+\varepsilon >1\), Zustand 3 existiert nicht:
\(\det J>0\) und \(\mbox{spur }J = -1+ \delta \left ( \dfrac{1}{\kappa+\varepsilon} -1\right ) <0\), der Zustand 2 ist stabil.
\(\kappa+\varepsilon =1\), Zustand 2 und 3 fallen zusammen:
\(\det J=0\), es ist keine Entscheidung über die Stabilität des Zustands 2 möglich.

Stabilität von \((B_3^\ast,R_3^\ast)=\left (\frac{\displaystyle 1}{\displaystyle 1-\kappa},\frac{\displaystyle 1-(\kappa+\varepsilon)}{\displaystyle(1-\kappa)^2} \right )\)

Der Zustand 3 existiert für \(\kappa+\varepsilon <1\).

\[ J(B_3^\ast,R_3^\ast) = \left[ \begin {array}{cc} -\varepsilon \dfrac{1+\kappa}{1-\kappa}+\kappa &-1\\ \delta \left (1- (\kappa+\varepsilon) \right ) &0 \end {array} \right] \]

Für \(\kappa+\varepsilon <1\) gilt \(\Delta=\det J= \delta \left (1- (\kappa+\varepsilon) \right ) > 0\).
Für \(-\varepsilon \dfrac{1+\kappa}{1-\kappa}+\kappa <0\) gilt \(\tau= \mbox{spur }J = -\varepsilon \dfrac{1+\kappa}{1-\kappa}+\kappa<0\).

Damit ist Zustand 3 für

\[\varepsilon> \kappa\dfrac{1-\kappa}{1+\kappa}\]

stabil.

Das Verhalten des Systems für \(\delta=1\), \(\kappa=0.5\) und verschiedene \(\varepsilon\) ist in Abbildung I.2 dargestellt.

Die Stabilität des 3. Fixpunkts in Abhängigkeit von \(\varepsilon\) und \(\kappa\) ist in Abbildung I.3 für \(\delta=0.2\) dargestellt.

Abbildung I.3: Stabilität des 3. Fixpunkts (Räuber und Beute existieren) des verallgemeinerten skalierten Räuber-Beute-Systems nach Rosenzweig-MacArthur in Abhängigkeit von \(\kappa\) und \(\varepsilon\) für \(\delta=0.2\).

I.4.1 Teilgleichgewichte (Isoklinen)

DasTeilgleichgewicht von B trennt die Bereiche des Zustandsraums, in denen die Bewegung nach rechts bzw. links erfolgt (senkrechte Linien in Abb. I.4). Das Teilgleichgewicht von R trennt die Bereiche des Zustandsraums, in denen die Bewegung nach oben bzw. unten erfolgt (nach unten geöffnete Parablen in Abb. I.4). Die Schnittpunkte der Teilgleichgewichte sind gerade die Fixpunkte.

I.4.2 Analyse des SIR-Modells

Zur Untersuchung des SIR-Modells (Abschnitt 13.1) wird vereinfachend eine konstante Bevölkerung \(N=S+I+R\) angenommen. In diesem Fall kann R aus den Zustandsgrößen S und I bestimmt werden und es reicht, das SI-Modell (13.2) zu untersuchen. Die Fixpunkte des Modells

\[\begin{aligned} \dot S &= - \beta \cdot \frac{I}{N} \cdot S \\[1ex] \dot I\; &= \phantom{-} \beta \cdot \frac{I}{N} \cdot(S -\rho) \end{aligned}\]

sind alle Punkte \((S^\ast,I^\ast)=(S,0), S\ge0\). Damit besteht die gesamte S-Achse aus Fixpunkten. Es kann somit keine asymptotische Stabilität erwartet werden, da die Fixpunkte dicht nebeneinander liegen. Die Linearisierung des Modells 13.2 ergibt die Jacobimatrix

\[ J(S^\ast,I^\ast) = \left[ \begin {array}{cc} -\beta\, I^\ast &-\beta \, S^\ast\\ \beta\, I^\ast &\beta \, (S^\ast-\rho) \end {array} \right] \]

Für \((S^\ast,I^\ast)=(S,0)\) gilt \(\tau=\mbox{spur }J =\beta\, (S-\rho)\) und \(\Delta=\det J =0\), und man erhält die Eigenwerte

\[\begin{aligned} \lambda_1&=& \beta\, (S-\rho)\\ \lambda_2&=& 0\\ \end{aligned}\]

Für \(S>\rho\) gilt \(\lambda_1>0\). In diesem Fall ist der Fixpunkt \((S,0)\) instabil, die Epidemie bricht aus. Da \(\lambda_2=0\), gibt es keine asymptotisch stabilen Fixpunkte. Für \(S<\rho\) ist \((S,0)\) (marginal) stabil, die Infizierten verschwinden.

Jetschke, Gottfried. 1989. Mathematik der Selbstorganisation. Wiesbaden: Vieweg+Teubner Verlag. https://doi.org/10.1007/978-3-322-85918-1.

Strogatz, Steven H. 2015. Nonlinear dynamics and chaos: with applications to physics, biology, chemistry, and engineering. Second edition. Boulder, CO: Westview Press, a member of the Perseus Books Group.