In diesem Modul werden gewöhnliche Differentialgleichungen sowie die Numerik gewöhnlicher Differentialgleichungen vorgestellt, die benötigt wird, um die Bewegungsgleichung des Feder-Masse-Modells numerisch zu lösen. Das Feld der Numerik ist recht umfangreich, so dass hier nur die Grundidee sowie einige Fachbegriffe erklärt werden. Des weiteren wird Matlab mit der grafischen Benutzeroberfläche Simulink als Simulationstool zur Lösung von gewöhnlichen Differentialgleichungen vorgestellt. Nach Durcharbeiten dieses Moduls soll der grundsätzliche Algorithmus zur numerischen Lösung gewöhnlicher Differentialgleichungen sowie dessen Anwendung in Matlab verstanden sein. (Eine eigenständige Programmierung von Algorithmen wird hier nicht angestrebt.)

Gewöhnliche Differentialgleichungen sind Gleichungen, die eine Funktion mit einer unabhängigen Variable und deren Ableitungen beinhaltet.

Anmerkung: Eine andere Klasse von Differentialgleichungen sind partielle Differentialgleichungen. Wie der Name beschreibt, beinhalten solche Differentialgleichungen partielle Ableitungen. Solche Gleichungen werden genutzt, um z.B. die Ausbreitung von Wellen zu beschreiben. Partielle Differentialgleichungen sollen hier aber nicht behandelt werden.

Im Folgenden werden gewöhnliche Differentialgleichungen und deren Lösungen als Beispiel vorgestellt. Der analytische Lösungsweg dieser Differentialgleichungen ist nicht Thema dieses Moduls und wird daher auch nicht präsentiert.

Die einfachste Differentialgleichung lautet:

;#; $ \dot{x} = Cx $ ;#;

In dieser Gleichung ist die unabhängige Variable die Zeit $t$. Die Variable $x$ ist von der Zeit $t$ abhängig und wird daher abhängige Variable genannt.

Die Lösung dieser Gleichung lautet:

;#; $ x = x_0 \:e^{C t} $ ;#;

Für positive Werte von $C$ beschreibt diese Gleichung ein exponentielles Wachstum (z.B. Wachstum von Bakterien), für negative Werte einen exponentiellen Zerfall (z.B. radioaktiver Zerfall). Der Parameter $x_0$ beschreibt den Zustand zum Zeitpunkt $t = 0$. Da in dieser Gleichung die Variable $x$ nur von der ersten Ableitung $\dot{ x}$ abhängt, bezeichnet man Gleichungen dieser Art als Differentialgleichungen erster Ordnung. Für eine Differentialgleichung erster Ordnung benötigt man genau eine Anfangsbedingung, um die Differentialgleichung eindeutig zu lösen.

Gewöhnliche Differentialgleichungen können nicht nur einfache Ableitungen sondern auch höhere Ableitungen beinhalten. Ein Beispiel für eine gewöhnliche Differentialgleichung mit höheren Ableitungen ist die Folgende:

;#; $ \ddot{ x} = -\omega^2\: x $ ;#;

Die Lösung dieser Gleichung lautet:

;#; $ x = x_0\: cos(\omega\: t) + \frac{v_0}{\omega} sin(\omega\: t) $ ;#;

Die Gleichung beschreibt eine harmonische Schwingung mit der Kreisfrequenz $\omega$. Diese Differentialgleichung ist eine Differentialgleichung zweiter Ordnung. Um die Lösung eindeutig zu bestimmen, benötigt man zwei Anfangsbedingungen. In diesem Fall den Ort $x_0$ zum Zeitpunkt $t=0$ sowie die Geschwindigkeit $v_0$ zum Zeitpunkt $t=0$.

Kommen in der Differentialgleichung mehr als nur eine Variable oder höhere Ableitungen vor, spricht man im allgemeinen von einem Differentialgleichungssystem. Um eine allgemeine Lösung für Differentialgleichungen anzugeben, werden Differentialgleichungen mit Ableitungen höherer Ordnung zu einem System von Gleichungen erster Ordnung reduziert. Hat man z.B. eine Gleichung mit Ableitungen bis zur 2-ten Ordnung:

;#; $ \ddot {y}= \dot{ y} +y = g(t,\dot{ y},y) $ ;#;

so kann man diese Gleichung in ein System von zwei Gleichungen umwandeln, indem man $\dot{ y}$ mit $x_2$ und $y$ mit $x_1$ substituiert. Somit erhält man folgendes Gleichungssystem:

;#; $ \dot x_1 =x_2\\\dot{ x_2} =g(t,x_1(t),x_2(t)) $ ;#;

Dieses System kann man nun in als Vektor zusammenfassen: ;#; $\begin{pmatrix}\dot{x_1} \\\dot{x_2}\end{pmatrix} =\dot{x}\: \: \: \: \: \: \: \: \: \: \:\begin{pmatrix}x_2 \\g(t,x_1(t),x_2(t))\end{pmatrix}=f(t,x(t)) $ ;#;

Damit erhält man in Kurzform:

;#; $ \dot x=f(t,x(t)) $ ;#;

mit dem Anfangsvektor:

;#; $ \dot x(t_0)=x_0 $ ;#;

Sowohl $x$ als auch $f(t,x(t))$ beschreiben hierbei einen Vektor bzw. eine Vektorfunktion. In dieser Form kann man die Lösung der Differentialgleichung allgemein angeben:

;#; $ x(t)=x_0+ \int_{t_0}^{t} \! f(t,x(t)) \, \mathrm{d}t $ ;#;

Diese Umwandlung gilt natürlich für beliebig hohe Ableitungen. In der Mechanik sind jedoch höhere als zweite Ableitungen fast nie anzutreffen. Um die Lösung $x(t)$ explizit anzugeben, muss man somit das Integral lösen. Im Falle der Bewegungsgleichungen des Feder-Masse-Modells ist das Integral analytisch nicht lösbar. Das heißt, dass es keine Stammfunktion für dieses Integral gibt und man dieses Integral nur numerisch lösen kann.

Die numerische Lösung von Differentialgleichungen kann nur diskretisiert erfolgen. D.h. man erhält keinen geschlossenen kontinuierlichen Term für die Lösung, sondern kann die Lösung nur numerisch an ausgewählten Stellen bestimmen. Für das Feder-Masse-Modell bedeutet dies, dass man die Lösung der Bewegungsgleichung nur zu bestimmten Zeitpunkten $t_i$ angeben kann. Der Abstand zwischen zwei Zeitpunkten wird Schrittweite genannt und mit $h$ bezeichnet.

Die numerische Integration erfolgt durch eine Abschätzung des Integrals. Das einfachste numerische Verfahren ist das Euler-Verfahren. Dabei wird das Integral durch eine Konstante über die Schrittweite $h$ abgeschätzt.

;#; $ x(t_{i+1})=x(t_i)+f(x(t_i),t_i)\: h $ ;#;

Je geringer dabei die Schrittweite h gewählt wird, desto geringer ist der Fehler der numerischen Lösung gegenüber der exakten Lösung.

<imgcaption image1 | Darstellung des Euler- und Heun-Verfahrens > </imgcaption>

Die Abschätzung des Integrals erfolgt durch näherungsweise Lösung mit Hilfe eines Polynoms. Dabei erfolgt die Abschätzung beim Euler-Verfahren mit einem Polynom nullter Ordnung d.h. einer Konstanten. Beim Heun-Verfahren erfolgt die Lösung mit einem Polynom erster Ordnung (siehe Bild). Die dazugehörige Rechenvorschrift lautet:

;#; $ x(t_{i+1})=x(t_i)+\frac{h}{2}(f(t_i,x_i)+f(t_{i+1},x_i+h \:f(t_i,x_i))) $ ;#;

Beim Runge-Kutta-Verfahren erfolgt die Abschätzung mit einem Polynom dritter Ordnung. Mit Runge-Kutta wird sowohl sowohl das spezifische Verfahren als auch die gesamte Klassen an Algorithmen bezeichnet. Man kann die Abschätzung des Integrals durch Polynome immer höherer Ordnung weiter verbessern, jedoch bedeutet dies auch einen höheren Rechenaufwand.

Es ist möglich, das Integral über einen großen Bereich mit der konstanten Schrittweite h numerisch zu bestimmen. Jedoch ist es sinnvoll, die Schrittweise in den Bereichen, in denen die Funktion f stark gekrümmt ist, zu verringern um den numerischen Fehler zu minimieren und die Schrittweite in den Bereichen, in denen die Funktion f wenig stark gekrümmt ist größer zu wählen um die Rechenzeit zu minimieren.

Dies erfolgt mittels einer Abschätzung des numerischen Fehlers, d.h. der Differenz zwischen exakter Lösung und numerischer Lösung. Da die exakte Lösung nicht vorliegt, wird der numerische Fehler geschätzt, in dem das Integral einmal durch ein Polynom bestimmter Ordnung numerisch bestimmt wird und dann nochmals durch ein Polynom höherer Ordnung. Die Differenz dieser beiden ergibt eine Abschätzung des numerischen Fehlers. Liegt der Fehler oberhalb einer gewissen Schwelle wird die Schrittweite für den nächsten Rechenschritt reduziert, bzw. der vorherige Rechenschritt wird mit geringerer Schrittweite wiederholt. Liegt der Fehler unterhalb einer gewissen Schwelle, wird der nächste Rechenschritt mit einer größeren Schrittweite durchgeführt. Ein Beispiel für so einen Algorithmus ist das Verfahren von Dormand-Prince. Der Matlab- Algorithmus ode45 nutzt dieses Verfahren.

Die Software Matlab wurde Ende der 1970er Jahre an der Universität New Mexico entwickelt. Das Ziel war es, den Studenten ein einfaches Programm zur Lösung von Problemen aus der linearen Algebra zur Verfügung zu stellen. Da in der linearen Algebra häufig mit Matrizen gearbeitet wird, leitet sich daher auch der Name ab: von MATrix LABoratory. Das Programm hat sich schnell auch an anderen Universitäten verbreitet, bis es ab 1984 zu einem kommerziellen Produkt wurde und durch viele Toolboxen wie z.B. Simulink erweitert wurde.

Eine erste Einführung in Matlab und dessen Benutzeroberfläche ist im Modul Matlab Basics beschrieben.

Matlab/Simulink besitzt mehrere Algorithmen zur Lösung von gewöhnlichen Differentialgleichungen (sowohl mit fester Schrittweite als auch variable Schrittweite, siehe Bild). Jeder dieser Algorithmen ist für spezifische Differentialgleichungen konzipiert. Für die Lösung der Bewegungsgleichungen des Feder-Masse-Modells hat sich der Algorithmus ode45 bewährt. Dieser Algorithmus ist für die meisten gewöhnlichen Differentialgleichungen geeignet und sollte der erste Algorithmus sein, der ausprobiert wird.

<imgcaption image2 | Implementierte Runge-Kutta-Verfahren mit fester und variabler Schrittweite in Matlab/Simulink > </imgcaption>

Für die Algorithmen können mehrere Einstellungen gewählt werden. Viele der Optionen können bei den Standard-Einstellungen belassen werden. Wichtige Optionen, die manuell gewählt werden sollten, sind die absolute und relative Genauigkeit. Diese beiden Werte sind standardmäßig auf $10^{-3}$ (d.h. einer Fehlertoleranz von 0,1%) gesetzt. Die Fehlertoleranz sollte dabei nicht größer sondern eher kleiner gewählt werden. Je nach Rechenleistung kann durchaus eine absolute und relative Fehlertoleranz von $10^{-12}$ gewählt werden. In der Wahl der Einstellungen ist jedoch auch häufig Heuristik gefragt.

Gewöhnliche Differentialgleichungen beinhalten eine oder mehrere abhängige Variablen sowie eine unabhängige Variable sowie die Ableitungen nach der unabhängigen Variablen. Diese Art der Gleichungen beschreiben vielfältige Effekte, wie z.B. Wachstums- und Zerfallprozesse, Schwingungen oder auch in der Biologie Räuber-Beute-Systeme. Die Lösung dieser Gleichung ist immer ein Integral. Dieses Integral muss jedoch keine Stammfunktion besitzen, so dass die Lösung nur näherungsweise durch numerische Verfahren abgeschätzt werden kann. Die Genauigkeit der Abschätzung ist dabei abhängig von dem gewählten numerischen Verfahren. Höhere Genauigkeit wird durch einen höheren Rechenaufwand erkauft. Matlab bietet mehrere Algorithmen, um gewöhnliche Differentialgleichungen numerisch zu lösen.

• Stelle die Funktion $x(t)=e^{Ct}$ für $C=1$und $C=-1$ dar. Welche Werte von $C$ beschreiben Wachstum und welche Werte von $C$ beschreiben Zerfall?
• Zeige durch Ableiten von $x(t)=x_0 \: e^{Ct}$ dass diese Funktion die Lösung der Differentialgleichung $\dot{ x}=Cx$ ist.

• Hermann, M. (2004). Numerik gewöhnlicher Differentialgleichungen: Anfangs-und Randwertprobleme. Oldenbourg Wissenschaftsverlag.

indexmenu_n_5