Möglichkeiten in MATLABDie numerische Lösung von Differentialgleichungen ist eine der großen Stärken von MATLAB (und Simulink). Hier wird lediglich auf die mehr oder weniger direkt auf dem Runge-Kutta-Verfahren basierenden Methoden zur Lösung von Anfangswertproblemen (ode23 und ode45) eingegangen. Darüber hinaus stehen aber auch Methoden zur Lösung von Randwertproblemen, Differentialgleichungen mit Verzögerungen und gewissen partiellen Differentialgleichungen zur Verfügung. ode23 und ode45 werden zur numerischen Lösung von gewöhnlichen Differentialgleichungen und von Systemen gewöhnlicher Differentialgleichungen verwendet. Diese sind auch Grundlage der Berechnungen von Simulink (vgl. entsprechendes Kapitel).
Form zur Lösung mit ode23: [T,Y] = ode23(@system_von_DGL, [t<sub>0</sub>, t<sub>final</sub>], [y<sub>0</sub>]) Rückgabeparameter der Funktion sind der berechnete Zeitvektor T und die zugehörige Wertematrix der Lösung(en) Y. Der Methode ode23 muss neben den Anfangsbedingungen y<sub>0</sub> und dem Intervall, in dem die Lösung berechnet werden soll ([t<sub>0</sub>, t<sub>final</sub>]), in einem m-File eine bzw. ein System von gewöhnlichen Differentialgleichungen übergeben werden. Der Dateiname wiederum wird dann der Methode ode23 als zu benutzende Funktion übergeben. Die Vorgehensweise ist bei ode45 und den anderen Verfahren analog. |
Anwendung der Lösungsmöglichkeiten - Beispiel mathematisches PendelBeispiel 1: Mathematisches Pendel
Differentialgleichung: \( \ddot \alpha (t) = - \large \frac {g} {l} \normalsize \; \sin \large ( \normalsize \alpha (t) \large ) \) mit \(g = 9,81 m/s^2\) Dabei ist \(α(t)\) der Winkel, den das Pendel der Länge l zur Zeit t bezogen auf die Ruhelage einnimmt (siehe Abb.). Da es sich um eine Differentialgleichung 2. Ordnung handelt, benötigen wir zu ihrer Lösung zwei Anfangsbedingungen \( \alpha (0) = \begin {pmatrix} \alpha (0) \\ \dot \alpha (0) \end {pmatrix} \) welche die Ausgangslage (Auslenkung) des Pendels und seine Anfangsgeschwindigkeit definieren. Diese Anfangsbedingungen werden als Vektor für y<sub>0</sub> an die Integrationsmethode übergeben. Zur Beschreibung der Differentialgleichung müssen wir in Matlab definieren, wie die Ableitungen der gesuchten Funktion α(t) voneinander abhängen. Dies müssen wir mit Hilfe der numerischen Werte der Funktion tun, d.h. mit den Vektoren der Funktion und ihren Ableitungen zu den Diskretisierungszeitpunkten. Zunächst muss bemerkt werden, dass MATLABs ode23 mit der DGL nur umgehen kann, wenn diese DGL in ein DGL-System 1.Ordnung umgewandelt wird. Dazu setzen wir \( \alpha_{1} (t) := \alpha (t) \) \( \alpha_{2} (t) := \dot \alpha (t) \) Die DGL lässt sich dann folgendermaßen umschreiben: \( \dot \alpha_{1} (t) = \alpha_{2} (t) \) \( \dot \alpha_{2} (t) = - \large \frac {g} {l} \normalsize \; \sin \large ( \normalsize \alpha_{1} (t) \large ) \) mit den Anfangsbedingungen: \( \alpha (0) = \begin {pmatrix} \alpha (0) \\ \dot \alpha (0) \end {pmatrix} = \begin {pmatrix} \alpha_{1} (0) \\ \alpha_{2} (0) \end {pmatrix} \) Auf diese Darstellung beziehen sich ode23 und die zugehörige Repräsentation der DGL in einem M-File. Wir tun dies hier zum Beispiel mit Hilfe einer MATLAB-Funktion pendgl, wobei wir die erste und die zweite Ableitung beschreiben und als Vektor zurückliefern. Die Struktur der Funktion muss dabei einer definierten Form genügen. Die Funktion hat folgende Gestalt:
Funktion als m-File herunterladen Der Parameter t muss in pendgl mit übergeben werden, auch wenn er funktionsintern nicht verwendet wird. Er wird von ode23 benötigt. ode23 verarbeitet auch erst die Anfangsbedingungen. Grundsätzlich müssen alle mit ode23 (oder anderen Methoden) zu lösenden Differentialgleichungen auf eine Form „Vektor der Ableitungen 1. Ordnung gleich Funktion des Lösungsvektors" gebracht und in einer MATLAB-Datei so dargestellt werden. Die folgenden Aufrufe liefern dann eine Darstellung der Lösung für die Anfangswerte \(\alpha (0) = \large \frac{\pi} {4} \normalsize \) und \( \dot \alpha (0) = 0 \) im Intervall [0,20] und für eine Pendellänge von 10. >> [t, loesung] = ode23(@pendgl, [0, 20], [pi/4, 0]) ; >> plot(t, loesung(:,1), 'r-', t, loesung(:,2), 'g--')
Anmerkung zur Lösung:
Der Lösungsalgorithmus arbeitet intern mit einer so genannten Schrittweitensteuerung. Die bedeutet, dass dort, wo die Lösung wenig variiert, große, und dort, wo sie stark variiert, kleine Schrittweiten verwendet werden. Dies erhöht die Effektivität der Algorithmen. Manchmal ist es jedoch wünschenswert, die Lösung mit einer einheitlichen Schrittweite weiterzuverarbeiten. In diesem Fall kann der Vektor, an dem die Lösung ausgegeben werde soll, beim Aufruf von ode23 explizit angegeben werden: >> [t, loesung] = ode23(@pendgl, (0:0.1:20), [pi/4,0]) ; Die Schrittweite hier beträgt nun 0,1. |
Anwendung der Lösungsmöglichkeiten - Beispiel RC-TiefpassBeispiel 2: RC-TiefpassDie folgende Abbildung zeigt das Schaltbild einer RC-Kombination mit Spannungsquelle:
Die Spannung am Ausgang des Systems genügt bekanntlich folgender linearer DGL 1. Ordnung: \( \dot u (t) = - \large \frac{1} {RC} \normalsize \; u(t) + \large \frac{1} {RC} \normalsize \; u_{1} (t) \) Auch hier wird die DGL mit Hilfe der Methode ode23 gelöst. Die benötigte MATLAB-Funktion sieht dann wie folgt aus:
Die verwendete Funktion u<sub>1</sub>(t) muss als MATLAB-Funktion definiert sein. Wir haben im vorliegenden Beispiel den Einheitssprung als Funktion u<sub>1</sub>(t) definiert und in folgender Weise im M-File u1.m im Arbeitsverzeichnis abgelegt:
Funktionen als m-Files herunterladen Mit den Befehlen: >> [t, loesung] = ode23(@rckomb, [0,3],0); >> plot(t, loesung) >> grid erhält man die in folgender Abbildung dargestellte Lösung im Zeitintervall [0,3] sec:
Anmerkung zur Lösung:
Im Programm wurden die die Werte R = 10kΩ und C = 4μF verwendet. Diese ergeben eine Zeitkonstante von τ = RC = 0,47s. |