Lineare Systeme

From bio-physics-wiki

Jump to: navigation, search

Contents

Lineare Systeme

Im Gegensatz zu den bisher behandelten Differenzialgleichungen, können Systeme von Differentialgleichungen nicht mehr einzeln gelöst und dann Zusammengefügt werden, sie müssen simultan gelöst werden, weil im Allgemeinen jede Gleichung mehrere unbekannte enthält. Wir schreiben allgemein für ein System in zwei Dimensionen

\begin{align} x'&=f(x,y,t)\\ y'&=g(x,y,t)\\ \end{align}

wobei $x$ und $y$ abhängige Variablen und $t$ die unabhängige Variable ist. Für lineare Systeme müssen die Funktionen $f$ und $g$ linear sein, d.h. wir können sie immer in der Form

\begin{align} x'&=ax + by + r_1(t)\\ y'&=cx + dy + r_2(t)\\ \tag{1} \end{align}

schreiben. Die Koeffizienten $a,b,c,d$ sind in linearen Systemen Konstanten. Man nennt Systeme in denen die Zeit $t$ nicht explizit vorkommt, bei denen also $r_1=0$ und $r_2=0$ sind, homogen. Wir behandeln zunächst homogene Lineare Systeme.

ODE Systeme als ODE höherer Ordnung

Jede ODE höherer Ordnung kann auf ein System erster Ordnung reduziert werden (siehe ODE höherer Ordnung als Systeme erster Ordnung). Wir wissen bereits, dass die Lösung einer ODE $n$-ter Ordnung $n$ freie parameter (Konstanten) besitzt. Weil beide Formulierungen äquivalent sind, hat auch das lineare System von ODEs erster Ordnung zwei Konstanten, die mit den Anfangsbedingungen $x(t_0)=x_0$ und $y(t_0)=y_0$ gefunden werden können.

Es ist genauso möglich ein lineares System in eine ODE zweiter Ordnung umzuschreiben. Betrachte z.B.

\begin{align} x'&=-2x+2y\\ y'&=2x-5y\\ \end{align}

Formen wir die erste Gleichung nach $y$ um erhalten wir $y=\frac{2x+x'}{2}$. Einsetzen von $y$ in die zweite Gleichung ergiebt dann

\begin{align} \left(\frac{2x+x'}{2} \right)'&=2x-5\left(\frac{2x+x'}{2} \right)\\ 2x'+x'' &=4x-10x-5x'\\ x''+7x'+6x &=0\\ \end{align}

eine ODE zweiter Ordnung die sich wie in Artikel Lineare Differentialgleichungen beschrieben lösen lässt. Mit dem Exponentialansatz erhalten wir für die Charakteristische Gleichung

\begin{align} \lambda^2+7\lambda+6 =(\lambda+1)(\lambda+6)=0\\ \end{align}

und damit die Lösungen

\begin{align} x(t)&=c_1 \cdot e^{-t} + c_2 \cdot e^{-6t} \\ y(t)=\frac{2x+x'}{2}&=c_1/2 \cdot e^{-t} - 2c_2 \cdot e^{-6t} \\ \tag{2} \end{align}

Beachte, dass beide Gleichungen die selben Konstanten $c_1$ und $c_2$ beinhalten.

ODE höherer Ordnung als Systeme erster Ordnung

Wir könne jede Differentialgleichung n-ter Ordnung (n>1) umschreiben in ein System von Differentialgleichungen 1. Ordnung.

Beispiel
\[y'''-3y''-5y'-2y=0\] Für Differentialgleichungen n-ter Ordnung schreiben wir nun einen (n)-dimensionalen Vektor, in unserem Beispiel ist $n=3$

$\begin{bmatrix} y''\\ y'\\ y\\\end{bmatrix}'= \begin{bmatrix} y'''\\ y''\\ y'\\\end{bmatrix}$

Wir versuchen nun unsere ODE in Matrixweise darzustellen. Dazu verwenden wir einen Trick. Wir schreiben nicht nur die ODE an, sonder zwei zusätzliche Gleichungen

\[y'''+3y''+5y'+2y=0\] \[y''=y''\] \[y'=y'\]

Jetzt können wir das Ganze in Matrixschreibweise bringen.

$\begin{bmatrix} y''\\ y'\\ y\\\end{bmatrix}'= \begin{bmatrix} 3 & 5 & 2\\1 & 0 & 0 \\ 0 & 1 &0 \end{bmatrix} \cdot \begin{bmatrix} y''\\ y'\\ y\\\end{bmatrix}$

Die erste Zeile der Matrix reproduziert unsere ODE. In der linken unteren Ecke sitzt die Einheitsmatrix.

Wir können also ganz allgemein Schreiben.

$\begin{bmatrix} y^{(n-1)}\\ y^{(n-2)}\\ \vdots \\ y'\\ y\\ \end{bmatrix}'=\begin{bmatrix} y^{(n)}\\ y^{(n-1)}\\ \vdots \\ y''\\ y'\\ \end{bmatrix}= \begin{bmatrix} -a_{(n-1)} & -a_{(n-2)} & \cdots & -a_1 &-a_{0} \\ 1 &0 & & 0 &0\\ 0 &1 & \vdots &\vdots &\vdots \\ \vdots &\vdots & \ddots & &\\ 0 &0 & \dots & 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} y^{(n-1)}\\ y^{(n-2)}\\ \vdots \\ y'\\ y\\ \end{bmatrix}$

Merke: Für eine ODE n-ter Ordnung bekommt man ein System von ODE's mit (n) Gleichungen.

Homogene Lösung - Eigenwertprobleme

Es gibt eine elegantere Methode lineare Systeme zu lösen als sie auf ODEs höherer Ordnung umzuschreiben. Vor allem, wenn es sich um große Systeme handelt, ist diese zu bevorzugen. Die Methode benutzt die Matrixschreibweise, wir können das System (1) auch schreiben als

\begin{align} \left( \begin{array}{c} \dot x(t) \\ \dot y(t)\\ \end{array} \right) = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \left( \begin{array}{c} x(t) \\ y(t)\\ \end{array}\right) \end{align}

oder kurz \begin{align} \dot{\vec{x}}=\textbf{A} \cdot \vec{x} \end{align}

Und die Lösung (2) in Vektorschreibweise angeben \begin{align} \left( \begin{array}{c} \dot x(t) \\ \dot y(t) \\ \end{array} \right)=c_1 \left( \begin{array}{c} 1 \\ 1/2\\ \end{array} \right) \cdot e^{-t} + c_2 \left( \begin{array}{c} 1 \\ -2\\ \end{array} \right) \cdot e^{-6t} \\ \tag{3} \end{align}

Wir führen nun die neue Methode am selben Beispiel ein, welches wir bereits durch zurückführen des ODE Systems auf eine ODE höherer Ordnung gelöst haben.

\begin{align} \left( \begin{array}{c} \dot x(t) \\ \dot y(t)\\ \end{array} \right) = \begin{bmatrix} -2 & 2 \\ 2 & -5 \end{bmatrix} \left( \begin{array}{c} x(t) \\ y(t)\\ \end{array}\right) \end{align}

Die Lösung (3) legt den Ansatz \begin{align} \vec{x}(t)=\left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right) \cdot e^{\lambda t} \\ \tag{4} \end{align} nahe. Setzen wir also diesen Ansatz ein in unser System ein, so erhalten wir

\begin{align} \lambda \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right) \cdot e^{\lambda t} = \begin{bmatrix} -2 & 2 \\ 2 & -5 \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right) \cdot e^{\lambda t} \end{align}

ein System von Linearen Gleichung. Die Lösung solcher Probleme die kurz Eigenwertprobleme genannt werden fällt in das Themengebiet der Linearen Algebra. Lösen wir also das Eigenwertproblem

\begin{align} \begin{bmatrix} -2 & 2 \\ 2 & -5 \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right) &= \lambda \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right) \\ \begin{bmatrix} -2 & 2 \\ 2 & -5 \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right) - \lambda \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right)&=0 \\ \left| \begin{bmatrix} (-2-\lambda) & 2 \\ 2 & (-5 - \lambda) \end{bmatrix} \right|&=0 \end{align}

Berechnung dieser Determinante ergiebt also

\begin{align} (-2-\lambda)(-5 - \lambda)-4&=0\\ \lambda^2+7\lambda+6 &=0\\ (\lambda+1)(\lambda+6)&=0\\ \end{align}

genau die selbe charakteristische Gleichung wie die ODE zweiter Ordnung. Wir haben also die Lösungen für $\lambda$ gefunden, es bleibt also noch die zu $\lambda_1=-1$ bzw. $\lambda_2=-6$ gehörigen Werte für $a_1$ und $a_2$ des Vektors zu finden. Dazu setzen wir die $\lambda$'s in \begin{align} \begin{bmatrix} -2 & 2 \\ 2 & -5 \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right) - \lambda \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right)&=\begin{bmatrix} (-2-\lambda) & 2 \\ 2 & (-5 - \lambda) \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right) =0 \\ \end{align}
ein, sodass wir für jedes $\lambda$ Gleichungen bekommen, die ein Vielfaches von einander, also linear abhängig sind. Für $\lambda_1$ erhalten wir dann die Gleichungen

\begin{align} \begin{bmatrix} (-2-\lambda_1) & 2 \\ 2 & (-5 - \lambda_1) \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right)= \begin{bmatrix} -1 & 2 \\ 2 & -4 \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right)=0 \\ \end{align}
die von \begin{align} \left( \begin{array}{c} a_{11} \\ a_{12}\\ \end{array}\right)=c_1 \cdot \left( \begin{array}{c} 1 \\ 1/2\\ \end{array}\right) \\ \end{align}
gelöst werden. Man nennt Vektoren $\left( \begin{array}{c} a_{1} \\ a_{2}\\ \end{array}\right) = \vec \alpha$ dessen vielfaches $c$ das linear abhängige Gleichungssystem lösen Eigenvektoren. Auf die selbe Art und Weise erhalten wir für $\lambda_2$ die Gleichungen

\begin{align} \begin{bmatrix} (-2-\lambda_2) & 2 \\ 2 & (-5 - \lambda_2) \end{bmatrix} \left( \begin{array}{c} a_{1} \\ a_{2}\\ \end{array}\right)= \begin{bmatrix} 4 & 2 \\ 2 & 1 \end{bmatrix} \left( \begin{array}{c} a_1 \\ a_2\\ \end{array}\right)=0 \\ \end{align}
und den Eigenvektor \begin{align} \left( \begin{array}{c} a_{21} \\ a_{22}\\ \end{array}\right)=c_2 \cdot \left( \begin{array}{c} 1 \\ -2\\ \end{array}\right) \\ \end{align}
Unser Ansatz (4) ergiebt also die Lösungen.

\begin{align} c_1 \cdot \left( \begin{array}{c} 1 \\ 1/2\\ \end{array}\right) \cdot e^{-t} \\ \end{align}
und \begin{align} c_2 \cdot \left( \begin{array}{c} 1 \\ -2\\ \end{array}\right) \cdot e^{-6t} \\ \end{align}

Die Allgmeine Lösung setzt sich gemäß dem Superpositionsprinzip aus der Linearkombination der beiden Lösungen zusammen. Wir erhalten das

\begin{align} \left( \begin{array}{c} \dot x(t) \\ \dot y(t) \\ \end{array} \right)=c_1 \left( \begin{array}{c} 1 \\ 1/2\\ \end{array} \right) \cdot e^{-t} + c_2 \left( \begin{array}{c} 1 \\ -2\\ \end{array} \right) \cdot e^{-6t} \\ \end{align}


Allgemein können wir Differentialgleichung also mit dem Ansatz
\begin{align}
 \vec{x}(t)=\vec{\alpha} \cdot e^{\lambda t} \\
 \tag{5}
 \end{align}
lösen. Indem man (5) in die Gleichung 
\begin{align} \dot{\vec{x}}=\textbf{A} \cdot \vec{x} 
 \end{align}
einsetzt, erhält man das das Eigenwertproblem.
\begin{align} \vec{\alpha}  \cdot \lambda \cdot e^{\lambda t} &=\textbf{A} \cdot \vec{\alpha}  \cdot e^{\lambda t}\\
 \\
 \textbf{A} \cdot \vec{\alpha} &= \lambda \cdot \vec{\alpha}\\
 \\
 (\textbf{A}-\lambda I) \cdot \vec{\alpha} &=0\\
 \\
 det(\textbf{A}-\lambda I)&=0\\
 \end{align}

wobei im letzten Schritt $\lambda$ mit der Einheitsmatrix $I$ multipliziert wurde, weil die Dimensionen nicht mit $\mathbf{A}$ übereinstimmen würden. Der Ausdruck $\lambda \cdot I \cdot \vec{\alpha}$ ist allerdings äquivalent mit $\lambda \cdot \vec{\alpha}$. Schreibt man dies allerdings nicht mit der Einheitsmatrix an, passiert es schnell dass einem dieser Fehler unterläuft.

Für zweidimensionale Systeme erhalten wir allgemein als Lösung der Determinante \begin{align} det\begin{bmatrix} (a-\lambda) & b \\ c & (d - \lambda) \end{bmatrix}=0 \end{align}
die charakteristische Gleichung \begin{align} \lambda^2 + \tau \lambda + \Delta=0 \end{align}
wobei $\tau$ die Spur (engl. trace) also die Summe der Diagonalelemente und $\Delta=det(\mathbf{A})=ad-bc$ ist, mit der Lösung

\begin{align} \lambda_{1,2} = \frac{\tau \pm \sqrt{\tau^2 - 4 \Delta}}{2} \end{align}
Folgende Eigenschaften sind eine große Hilfe wenn es darum geht schnell eine Lösung für die quadtratische Gleichung zu finden. \begin{align} \lambda_1 \lambda_2 = \Delta \hspace{2cm} \lambda_1+ \lambda_2 = \tau \end{align}

LinDynSysTraceDet.jpg

Phasenportrait

Man kann das Verhalten der Lösungen in einem Phasenportrait darstellen. Im zweidimensionalen Fall spricht man auch von der Phasenebene. Mit hilfe von Phasenportraits bekommt man einen intuitiven Zugang zu Systemen und die Möglichkeit Lösungen zu veranschaulichen. Dazu ordenen wir gemäß

\begin{align} \left( \begin{array}{c} \dot x(t) \\ \dot y(t)\\ \end{array} \right) = \begin{bmatrix} a & b \\ c & d \end{bmatrix} \left( \begin{array}{c} x(t) \\ y(t)\\ \end{array}\right) \end{align}

jedem Punkt $(x,y)$ in der Phasen ebene einen Vektor $\dot{\vec{x}},\dot{\vec{y}}$ zu. Lösungen der Differentialgleichung sind Tangentiallinen an dieses Vektorfeld. In einem Phasenportrait ist oft das Vektorfeld selbst nicht eingezeichnet, sondern nur die Lösungen mit der Richtung des Flusses (engl. flow) der die Richtung angibt, in welche das Vektorfeld zeigen würde. Ist kein Vektorfeld eingezeichnet, sieht man also nur die Richtung des Flusses nicht aber sein Betrag. Sind die Lösungen für $\lambda_{1,2}$ nicht komplex, haben also keinen imaginärteil, gibt es immer eine Richtung, bei der der Fluss eine Gerade bildet. Diese Richtungen werden Eigenrichtungen (siehe Abb. in rot) genannt und werden von den Eigenvektoren des Systems angegeben.

Eine Erklärung wie man Phasenportraits ohne Kenntnis von Vektorfeldern zeichnet findet man hier

Stabilität von Systemen - Klassifizierung in der Phasenebene

Die Lösung eines zweidimensionalen Systems mit der charakteristischen Gleichung \begin{align} \lambda^2 - \tau \lambda + \Delta=0 \end{align}
mit

\begin{align} \lambda_{1,2} = \frac{\tau \pm \sqrt{\tau^2 - 4 \Delta}}{2} \end{align}
kann je nach ihrer Form in Verbindung mit einem spezifischen Verhalten bzw. Flüssen im Phasenportrait gebracht werden. Man unterscheidet die Klassen

  • rein reelle Eigenwerte ($\tau^2-4\Delta>0$)
    • stabile Knoten (=Fixpunkt=engl. fixed point): $\lambda_1 , \lambda_2 < 0$
    • instabile Knoten: $\lambda_1 , \lambda_2 > 0$
    • metastabiler Knoten (Sattelpunkt): $\lambda_1<0$ (zugehöriger Eigenvektor gibt stabile Richtung an) und $\lambda_2 > 0$ (zugehöriger Eigenvektor gibt instabile Richtung an)
    • degenerierter Knoten: $\tau^2-4\Delta=0$, nur einen Eigenvektor und nur eine Eigenrichung.
    • nicht isolierte Fixpunkte: $\Delta=0$
  • rein komplexe Eigenwerte ($\tau=0$)
    • Zentrum (engl. center) Lösung bildet eine Ellipse
  • gemischte Eigenwerte: Eigenwert besitzt Real- und Imaginärteil $\lambda=a \pm ib$.
    • stabile Spirale: $a<0$
    • instabile Spirale: $a>0$
  • gleiche Eigenwerte: $\lambda_1 = \lambda_2 \neq 0$, Lösung bildet einen Stern



Ein nettes Tool zur Erstellung von Phasenprotraits gibts hier.

Anfangswertprobleme für Systeme - Fundamentalmatrix

Das System \begin{align} \dot{\vec{x}}=\textbf{A} \cdot \vec{x} \end{align}

in zwei dimensionen hat die Lösung \begin{align} c_1 \cdot {\vec{x}_1}(t)+ c_2 \cdot \vec{x}_2(t) \end{align}

Die Vektoren $\{ \vec{x}_1,\vec{x}_2 \}$ sind linear unabhängig und bilden die Basis des zweidimendionalen Lösungsraums. Sie werden als Fundamentalsystem bezeichnet. Die Matrix, deren Lösungensvektoren $\vec{x}_1$ und $\vec{x}_2$ die Spalten bilden, heißt Fundamentalmatrix des Systems und wird

\begin{align} X(t)= [ \vec{x}_1 \vec{x}_2 ]=\begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \end{bmatrix} \end{align} geschrieben. Beachte, dass die Fundamentalmatrix nicht eindeutig ist, weil jede Lösung $\mathbf{x}_i$ mit einer Konstanten multipliziert wiederum einer Lösung Die Determinante einer Fundamentalmatrix wird als Wronskideterminate bezeichnet.

\begin{align} |X(t)|= \left| \begin{bmatrix} x_1 & y_1 \\ x_2 & y_2 \end{bmatrix} \right| \end{align}

Mithilfe einer Fundamentalmatrix kann man die allgemeine Lösung in Matrixschreibweise anschreiben.

\begin{align} x= c_1 \begin{pmatrix} x_1 \\ y_1 \end{pmatrix} + c_2 \begin{pmatrix} x_2 \\ y_2 \end{pmatrix} = \begin{bmatrix} x_1 & x_1 \\ x_2 & y_2 \end{bmatrix}\begin{pmatrix} c_1 \\ c_2 \end{pmatrix} \end{align} oder kurz \begin{align} \mathbf{x}= \mathbf{X}(t) \mathbf{c} \tag{6} \end{align} wobei die Notationen $\vec{x}$ und $\mathbf{x}$ die selbe Bedeutung haben. Die Fundamentalmatrix erlaubt es Anfangswertprobleme für große Systeme bequem mit Matrixalgebra zu lösen. Wir gehen vom System

\begin{align} \dot{\mathbf{x}}=\textbf{A} \mathbf{x} \hspace{2cm}\mathbf{x}(t)=\mathbf{x}_0 \end{align}

und den Anfangsbedingungen $\mathbf{x}(t)=\mathbf{x}_0$ mit der allgemeinen Lösung (7) aus und suchen einen Vektor $\mathbf{c}$ der die Anfangsbedingungen erfüllt. Einsetzen der Anfangsbedingungen ergiebt \begin{align}\mathbf{X}(t_0) \mathbf{c}=\mathbf{x}_0 \tag{7} \end{align} Wenn die Lösungen eines $n$ dimensionalen Systems die das Fundamentalsystem $\{\vec{x}_1,\vec{x}_2, \dots ,\vec{x}_n\}$ bilden linear unabhängig sind, ist die Wronskydeterminante $|\mathbf{X}(t_0)| \neq 0$ und daher die inverse Matrix $\mathbf{X}(t_0)^{-1}$ existiert, erhalten wir $\mathbf{c}$ \begin{align}\mathbf{c}=\mathbf{X}(t_0)^{-1} \mathbf{x}_0 \end{align} Durch einsetzten in (6) erhalten wir dann \begin{align} \mathbf{x}(t)= \mathbf{X}(t) \mathbf{X}(t_0)^{-1} \mathbf{x}_0 \end{align}

Außerdem gilt für die Fundamentalmatrix $\mathbf{X}(t)$, dass sie die Matrixdifferentialgleichung \begin{align} \dot{\mathbf{X}}= \mathbf{A}\mathbf{X} \end{align} erfüllt, denn $\mathbf{X}(t)$ enthält die Lösungen $\mathbf{x}_1, \mathbf{x}_2$. Beweis: \begin{align} \mathbf{A}\mathbf{X} = \mathbf{A} (\mathbf{x}_1 \mathbf{x}_2)= (\mathbf{A} \mathbf{x}_1 \mathbf{A} \mathbf{x}_2)=(\dot{\mathbf{x}}_1 \dot{\mathbf{x}}_2)=\dot{\mathbf{X}} \end{align}

Inhomogene Lösung - Variation der Konstanten

Für das homogene System \begin{align} \dot{\mathbf{x}}=\textbf{A} \cdot \mathbf{x} \end{align}

in zwei Dimensionen hatten wir die Lösung \begin{align}\mathbf{x}= c_1 \cdot {\mathbf{x}_1}(t)+ c_2 \cdot \mathbf{x}_2(t) \end{align} gefunden. Nun wollen wir das inhomogene System \begin{align} \dot{\mathbf{x}}=\textbf{A} \cdot \mathbf{x} + \mathbf{r}(t) \tag{8} \end{align} lösen, dazu bedienen wir uns der Methode der Varitation der Konstanten. Diese Methode nutzt die Lösung des homogenen Systems, nimmt aber an, das die Konstanten $c_1,c_2$ Variablen sind, die wir $v(t)$ nennen. So erhält man den Ansatz \begin{align}\mathbf{x}= v_1(t) \cdot {\mathbf{x}_1}(t)+ v_2(t) \cdot \mathbf{x}_2(t) \end{align} den wir mithilfe der Fundamentalmatrix schreiben können als \begin{align} \mathbf{x}=(\mathbf{x}_1 \mathbf{x}_2)\mathbf{v}=\mathbf{X}\mathbf{v} \tag{9} \end{align} Einsetzen in die Differentialgleigung (8) liefert wegen der Produktregel \begin{align} \dot{\mathbf{X}}\mathbf{v}+\mathbf{X} \dot{\mathbf{v}}=\textbf{A} \mathbf{X}\mathbf{v} + \mathbf{r} \end{align} Da die Matrixdifferentialgleichung \begin{align} \dot{\mathbf{X}}= \mathbf{A}\mathbf{X} \end{align} gilt, können wir kürzen und erhalten \begin{align} \mathbf{X} \dot{\mathbf{v}}= \mathbf{r} \end{align} Durch umformen und integrieren erhalten wir \begin{align} \mathbf{v}(t)=\int \mathbf{X}^{-1} \mathbf{r} \hspace{0.1cm} dt \end{align} sodass wir schließlich $\mathbf{v}(t)$ in (9) einsetzen können \begin{align} \mathbf{x}(t)=\mathbf{X} \int \mathbf{X}^{-1} \mathbf{r} \hspace{0.1cm} dt \end{align}

Homogene Lösung - Matrixexponential

Die Berechnung des Fundamentalsystems mithilfe des Matrixexponentials ist die eleganteste Methode. Sie erlaubt es die Lösung eines Systems in einer einzigen Formel darzustellen, während wir bisher die Eigenwerte und Eigenvektoren berechnen mussten.

Sei $\mathbf{A}$ eine reelle oder komplexe $n \times n$-Matrix. Das Exponential von $\mathbf{A}$, welches durch $e^{\mathbf{A}}$ bezeichnet wird, ist die $n \times n$-Matrix, welche durch die folgende Potenzreihe definiert ist.

\[e^A = \sum_{k=0}^\infty\frac{X^k}{k!}. \tag{10}\]

Diese Reihe konvergiert immer. Daher ist das Exponential von $\mathbf{A}$ wohldefiniert. Wenn $\mathbf{A}$ eine 1×1-Matrix ist, entspricht das Matrixexponential von $\mathbf{A}$ der gewöhnlichen Exponentialfunktion.

Für den Fall, dass A eine 1×1-Matrix ist, können wir die Differentialgleichung $\dot{x}=A x$ durch Seperation der Variablen lösen. Die Lösung ist dann $x(t)=c \cdot e^{At}$. Wir werden nun sehen, dass dies ihn ähnlicher Weise auch möglich ist, wenn $\mathbf{A}$ eine Matrix ist, denn gemäß (10) erhalten wir für $e^{\mathbf{A}t}$

\begin{align} e^{\mathbf{A}t}=I+\mathbf{A}t+\mathbf{A}\frac{t^2}{2!}+\mathbf{A}\frac{t^3}{3!}+ \dots \tag{11} \end{align}

und haben damit einen Ausdruck für das Matrixexponential. Angenommen das Matrixexponential $e^{\mathbf{A}t}$ hat die Dimension $n \times n$, dann ist sowohl die quadratische Matrix $\mathbf{A}$, wie auch der Lösungsvektor $\mathbf{x}(t)$ von der Dimension $n$. Dies bedeutet, dass die Konstante (im eindimensionalen Fall $c$) ebenfalls ein Vektor der Dimension $n$ sein muss. Diese Konstante $\mathbf{c}$ muss von rechts mulitpliziert werden, dass die Dimensionen übereinstimmen.

\begin{align} \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{c} \tag{12} \end{align}

Vergleich dieser Lösung mit der im Artikel Anfangswertprobleme für Systeme - Fundamentalmatrix eingeführten Schreibweise für die Lösung des Systems $\mathbf{x}=\mathbf{X}\mathbf{c}$ zeigt, dass ein Zusammenhang zwischen der Matrix $e^{\mathbf{A}t}$ und der Fundamentalmatrix besteht. Wie genau diese Beziehung aussieht erfahrt man am Ende dieses Artikels.

Das Matrixexponential kann durch einsetzen der Matrix $\mathbf{A}$ in die Taylorreihenentwicklung (11) berechnet werden. Mit einem Computeralgebrasystem wie Matlab oder Mathematica kann das Matrixexponential unmittelbar durch Eingabe von $e^{\mathbf{A}t}$ berechnet werden. Das erspart einem den Umweg über die Berechnung der Eigenwerte und Eigenvektoren.

Beispiel
\begin{align} \mathbf{A}= \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \end{align} Für die Taylorreihenentwicklung brauchen wir $\mathbf{A}^2,\mathbf{A}^3,\dots ,\mathbf{A}^n$. Für $\mathbf{A}^2$ erhalten wir \begin{align} \mathbf{A}^2= \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \cdot \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix}=\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}\end{align} die Einheitsmatrix, alle weiteren Lösungen sind daher einfach zu finden und haben die Form \begin{align} \mathbf{A}^{2k}&= \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \hspace{2cm} k=1,2, \dots\\ \mathbf{A}^{2k+1}&= \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \hspace{2cm} k=1,2, \dots \\ \end{align} Einsetzen von $\mathbf{A}^{n}$ in (11) ergibt \begin{align} e^{\mathbf{A}t}&=\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}+\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \cdot t+\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \frac{t^2}{2!}+\begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix} \frac{t^3}{3!}+ \dots \\ e^{\mathbf{A}t}&=\begin{bmatrix} 1+\frac{t^2}{2!}+\frac{t^4}{4!} + \dots & t+\frac{t^3}{3!}+ \frac{t^5}{5!}+ \dots \\ t+\frac{t^3}{3!}+ \frac{t^5}{5!}+ \dots & 1+\frac{t^2}{2!}+\frac{t^4}{4!} + \dots \end{bmatrix}=\begin{bmatrix} e^t + e^{-t} & e^t - e^{-t} \\ e^t - e^{-t}& e^t + e^{-t} \end{bmatrix}=\begin{bmatrix} cosh(t) & sinh(t) \\ sinh(t) & cosh(t) \end{bmatrix} \end{align}

In diesem Beispiel war die Berechnung von $e^{\mathbf{A}t}$ analytisch möglich, in vielen Fällen ist das nicht möglich. Wir untersuchen deshalb mögliche Situationen die bei der händischen Berechnung probleme bereiten können. Zuerst untersuchen wir aber die Eigenschaften des Matrixexponentials. Im Allgemeinen gilt \begin{align} e^{\mathbf{A}+\mathbf{B}}& \neq e^{\mathbf{A}}e^{\mathbf{B}} \end{align} Für den speziellen Fall, dass $\mathbf{A}$ und $\mathbf{B}$ Kommutieren, d.h. $[\mathbf{A},\mathbf{B}]=\mathbf{A}\mathbf{B}-\mathbf{B}\mathbf{A}=0$ ist, gilt die Gleichung \begin{align} e^{\mathbf{A}+\mathbf{B}}& = e^{\mathbf{A}}e^{\mathbf{B}} \end{align} Außerdem ist wie im eindimsionalen Fall $e^0=1$ \begin{align} e^{\mathbf{0}}=\mathbf{I} \end{align} wobei $\mathbf{I}$ die $n$-dimensionale Einheitsmatrix ist. Außerdem gilt für skalare $a,b$ dass \begin{align} e^{a\mathbf{X}}e^{b\mathbf{X}}=e^{(a+b)\mathbf{X}} \end{align} und damit \begin{align} e^{1\mathbf{X}}e^{-1\mathbf{X}}=e^{(1-1)\mathbf{X}}= e^{\mathbf{0}}=\mathbf{I} \end{align} deshalb muss $e^{-1\mathbf{X}}$ die Inverse $\left( e^{1\mathbf{X}} \right)^{-1}$ von $e^{1\mathbf{X}}$ sein. \begin{align} e^{-1\mathbf{X}}&=\left( e^{1\mathbf{X}} \right)^{-1} \end{align}

  • Die Berechnung von $e^{\mathbf{A}t}$ ist im Allgemeinen analytisch nicht möglich, da die Reihenentwicklung (11) unendlich viele Terme hat.
  • Für den Fall das $\mathbf{A}$ die Form

\begin{align} \mathbf{A}= \begin{bmatrix} a & b \\ b & a \end{bmatrix} \end{align}

und die Matrizen

\begin{align} \mathbf{A}= \underbrace{\begin{bmatrix} a & 0 \\ 0 & a \end{bmatrix}}_{\mathbf{X}} + \underbrace{\begin{bmatrix} 0 & b \\ b & 0 \end{bmatrix}}_{\mathbf{Y}} \end{align}

aufgespalten werden kann. Ist es möglich $e^{\mathbf{A}t}$ zu schreiben als

\begin{align} e^{\mathbf{A}}=e^{\mathbf{X}+\mathbf{Y}}=e^{\mathbf{X}}e^{\mathbf{Y}} \end{align}

weil in diesem Fall $\mathbf{X}$ und $\mathbf{Y}$ kommutieren.

\begin{align} \begin{bmatrix} a & 0 \\ 0 & a \end{bmatrix} \begin{bmatrix} 0 & b \\ b & 0 \end{bmatrix} - \begin{bmatrix} 0 & b \\ b & 0 \end{bmatrix} \begin{bmatrix} a & 0 \\ 0 & a \end{bmatrix} = \begin{bmatrix} ab & ab \\ ab & ab \end{bmatrix} -\begin{bmatrix} ab & ab \\ :ab & ab \end{bmatrix}=\begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}\end{align}

  • Für alle anderen Fälle muss $e^{\mathbf{A}t}$ über den mühsamen Weg gefunden werden, indem man die Eigenvektoren berechnet.

\begin{align} e^{\mathbf{A}t}= \mathbf{X}(t)\mathbf{X}(0)^{-1} \end{align}

Die Matrix $\mathbf{X}(t)\mathbf{X}(0)^{-1}$ hat die selben Eigenschaften wie $e^{\mathbf{A}t}$, sodass z.B.

\begin{align} e^{\mathbf{A}0}=\mathbf{X}(0)\mathbf{X}(0)^{-1}=\mathbf{I} \end{align}

die Einheitsmatrix ist.

Das Anfangswertproblem $\mathbf{x}(t=0)=\mathbf{x}_0$ des Systems $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}$ nimmt mit der Methode des Matrixexponentials eine besonders einfache Form an, denn Einsetzen der Angfangswerte in die Lösung (12) ergibt. \begin{align} \mathbf{x}(t)=e^{\mathbf{A}0}\mathbf{c}=\mathbf{c}=\mathbf{x}_0 \tag{13} \end{align} Der Vektor $\mathbf{c}$ der die Konstanten enthält, entspricht also einfach dem Vektor $\mathbf{x}_0$ der Anfangsbedingung. Wir erhalten damit die homogene Lösung \begin{align} \mathbf{x}(t)=e^{\mathbf{A}t}\mathbf{x}_0 \tag{14} \end{align}

Homogene Lösung - Entkoppeln

Im Allgemeinen sind homogene ODE Systeme von der Form \[ \dot{x_1} = \frac{dx_1}{dt} = f_1(x_1, \ldots, x_n) \]

\[ \vdots \]

\[ \dot{x_i} = \frac{dx_i}{dt} = f_i(x_1, \ldots, x_n) \]

\[ \vdots \]

\[ \dot{x_n} = \frac{dx_n}{dt} = f_n(x_1, \ldots, x_n) \]

sodass die Gleichungen gleichzeitig gelöst werden müssen weil $f_i(x_1, \ldots, x_n)$ nicht nur von $x_i$ abhängt. Man spricht davon, dass die Differentialgleichungen gekoppelt sind. Man kann das System entkoppeln und damit in die Form

\[ \dot{x_1} = \frac{dx_1}{dt} = f_1(x_1) \]

\[ \vdots \]

\[ \dot{x_i} = \frac{dx_i}{dt} = f_i(x_i) \]

\[ \vdots \]

\[ \dot{x_n} = \frac{dx_n}{dt} = f_n(x_n) \]

bringen, indem man das System mithilfe der Matrix $\mathbf{E}$ in ein anderes Koordinatensystem transformiert. In diesem neuen Koordinatensystem sind die Gleichungen entkoppelt. Man erhält $n$ Gleichungen erster Ordnung, die einzeln und unabhängig voneinander gelöst werden können. Welche Abbildung $\mathbf{E}$ \begin{align} \textbf{x}&=\textbf{E} {\mathbf{u}}\\ \textbf{u}&=\textbf{E}^{-1} {\mathbf{x}} \end{align}

erfüllt die Bedingung \begin{align*} \dot{\mathbf{x}}&=\mathbf{A} \cdot \mathbf{x} \tag{15}\\ \mathbf{E}\dot{\mathbf{u}}&=\mathbf{AE} \cdot \mathbf{u} \tag{16}\\ \mathbf{E}\dot{\mathbf{u}}&=\mathbf{E}\mathbf{\Lambda} \cdot \mathbf{u} \tag{17}\\ \mathbf{E}^{-1}\mathbf{E}\dot{\mathbf{u}}&=\mathbf{E}^{-1}\mathbf{E}\mathbf{\Lambda} \cdot \mathbf{u} \tag{18}\\ \dot{\mathbf{u}}&=\mathbf{\Lambda} \cdot \mathbf{u} \tag{19} \end{align*} Wobei $\mathbf{\Lambda}$ eine Diagonalmatrix z.B.

\begin{align} \mathbf{\Lambda} = \begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix} \end{align}

ist, sodass $\mathbf{E}$ die Matrix $\mathbf{A}$ entkoppelt. Vergleichen wir Gl. (16) mit Gl. (17) so muss gelten. \begin{align*} \mathbf{AE} &=\mathbf{E}\mathbf{\Lambda} \\ \mathbf{A}(\mathbf{e}_1\mathbf{e}_2) &=(\mathbf{e}_1\mathbf{e}_2)\mathbf{\Lambda} \\ \mathbf{A}\mathbf{e}_1 &=\lambda_1 \mathbf{e}_1 \\ \mathbf{A}\mathbf{e}_2 &=\lambda_2 \mathbf{e}_2 \\ \end{align*} Die Gleichung ist also erfüllt wenn $\mathbf{e}_i$ Eigenvektoren zu den Eigenwerten $\lambda_i$ sind. Wir verwenden im Folgenden die gewohnte Notation in welcher wir die Eigenvektoren von $\mathbf{A}$ mit $\vec{\alpha}_i$ bezeichnet haben.

Die spalten der Matrix $\mathbf{E}$, welche die Matrix $\mathbf{A}$ entkoppelt, bestehen aus den Eigenvektoren $\vec{\alpha}_i$ der Matrix $\mathbf{A}$

Beispiel
\begin{align} \left( \begin{array}{c} \dot x \\ \dot y\\ \end{array} \right) = \underbrace{\begin{bmatrix} -2 & 2 \\ 1 & -1 \end{bmatrix}}_{\mathbf{A}} \tag{20} \left( \begin{array}{c} x \\ y\\ \end{array}\right) \end{align}

Wir entkoppeln nun dieses System, dazu suchen wir die Eigenvektoren der Matrix $\mathbf{A}$ \begin{align} \lambda^2 - trace(\mathbf{A})\lambda+det(\mathbf{A}) = \lambda^2 -3\lambda=0 \hspace{2cm} \lambda_1= 0;\lambda_2=-3 \end{align} \begin{align} \vec{\alpha}_1= \begin{pmatrix} 1 \\ 1 \end{pmatrix} \hspace{2cm} \vec{\alpha}_2= \begin{pmatrix} 2 \\ -1 \end{pmatrix} \end{align} Die Spalten der Entkoppelungsmatrix $\mathbf{E}$ bestehen aus den Eigenvektoren $\mathbf{E}=(\vec{\alpha}_1 \vec{\alpha}_2)$

\begin{align} \mathbf{E}= \begin{bmatrix} 1 & 2 \\ 1 & -1 \end{bmatrix} \end{align} Um die Gleichungen (15) - (19) besser zu verstehen, führen wir dieselben Bechrechnungen am Beispiel (20) dieses Systems durch. Dazu ersetzen wir $\dot{\mathbf{x}}$ in Gl. (20) mit $\mathbf{E} \dot{\mathbf{u}}$ \begin{align} \left( \begin{array}{c} \dot x \\ \dot y\\ \end{array} \right) = \underbrace{\begin{bmatrix} 1 & 2 \\ 1 & -1 \end{bmatrix}}_{\mathbf{E}=\begin{bmatrix} \vec{\alpha}_1 \vec{\alpha}_2 \end{bmatrix}} \underbrace{\left( \begin{array}{c} \dot u \\ \dot v\\ \end{array} \right)}_{\dot{\mathbf{u}}} &= \underbrace{\begin{bmatrix} -2 & 2 \\ 1 & -1 \end{bmatrix}}_{\mathbf{A}} \underbrace{\begin{bmatrix} 1 & 2 \\ 1 & -1 \end{bmatrix}}_{\mathbf{E}=\begin{bmatrix} \vec{\alpha}_1 \vec{\alpha}_2 \end{bmatrix}} \left( \begin{array}{c} u \\ v\\ \end{array}\right)= \underbrace{\begin{bmatrix} 1 & 2 \\ 1 & -1 \end{bmatrix}}_{\mathbf{E}=\begin{bmatrix} \vec{\alpha}_1 \vec{\alpha}_2 \end{bmatrix}} \underbrace{\begin{bmatrix} 0 & 0 \\ 0 & -3 \end{bmatrix}}_{\begin{bmatrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{bmatrix}} \left( \begin{array}{c} u \\ v\\ \end{array}\right)\\ \left( \begin{array}{c} \dot u \\ \dot v\\ \end{array} \right) &= \begin{bmatrix} 0 & 0 \\ 0 & -3 \end{bmatrix} \left( \begin{array}{c} u \\ v\\ \end{array}\right)\\ \dot{u}&=0 \Rightarrow u=const.=c_1\\ \dot{v}&=-3 \Rightarrow v=c_2 \cdot e^{-3t}\\ \end{align}


Um $\mathbf{u}$ zu berechnen benötigen wir die $\textbf{E}^{-1}$ \begin{align} \mathbf{x}=\textbf{E}\mathbf{u} = \begin{bmatrix} 1 & 2 \\ 1 & -1 \end{bmatrix} \left( \begin{array}{c} c_1 \\ c_2 \cdot e^{-3t} \end{array} \right) =\left( \begin{array}{c} c_1+ 2c_2 \cdot e^{-3t} \\ c_1 - c_2 \cdot e^{-3t} \\ \end{array} \right) \end{align}

Allgemeine Lösung

Die Allgemeine Lösung eines linearen Systems von gewöhnlichen Differentialgleichungen mit Anfangsbedingungen $\mathbf{x}(0)=\mathbf{x}_0$ \begin{align} \dot{\mathbf{x}}=\mathbf{A}\mathbf{x} + \mathbf{r}(t) \end{align} besteht aus der Lösung des homogenen Systems $\mathbf{x}_c$ mit $\mathbf{r}(t)=0$ die auch komplementäre Lösung (engl. complementary) genannt wird \begin{align} \dot{\mathbf{x}}_c(t)=e^{\mathbf{A}t}\mathbf{x}_0 \end{align} und der inhomognen Lösung $\mathbf{x}_p$, auch partikuläre Lösung (engl. particular) genannt, die durch Variation der Konstanten gefunden werden kann. \begin{align} \mathbf{x}_p(t)=\mathbf{X}(t) \int \mathbf{X}^{-1}(t) \mathbf{r}(t) \hspace{0.1cm} dt \end{align} Da die Fundamentalmatrix nicht eindeutig bestimmt ist, bildet auch die mit $\mathbf{X}^{-1}(0)$ nomalisierte Matrix $\mathbf{X}(t)\mathbf{X}^{-1}(0)$ eine Fundamentalmatrix. Diese entspricht genau dem Matrixexponential $e^{\mathbf{A}t}$ und wir können $\mathbf{x}_p$ auch schreiben als \begin{align} \mathbf{x}_p(t)=e^{\mathbf{A}t} \int e^{- \mathbf{A}t} \mathbf{r}(t) \hspace{0.1cm} dt \end{align}

Allgemeine Lösung

\begin{align}\mathbf{x}(t)=\mathbf{x}_c(t)+ \mathbf{x}_p(t)=e^{\mathbf{A}t}\mathbf{x}_0 + e^{\mathbf{A}t} \int e^{- \mathbf{A}t} \mathbf{r}(t) \hspace{0.1cm} dt
\end{align}