Zurück
Vor 

Möglichkeiten in MATLAB

Die 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 Pendel

Beispiel 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:
function [alphadot] = pendgl(t,alpha)
Funktion pendgl
Aufruf über:
  [alphadot]= pendgl(t,alpha)
Beispiel für die Lösung von DGL's unter MATLAB mit ode23. Dieese Funktion definiert die DGL des mathematischen Pendels. Damit MATLAB damit umgehen kann, muss die DGL 2. Ordnung vorher in ein System von DGL 1. Ordnung überführt werden.
l = 10 ;
g = 9.81 ;
Definition der Konstanten:
l = Pendellänge
g = Erdbeschleunigung

Diese Parameter können auch als weitere Parameter an die Definitionsfunktion der Differentialgleichungen übergeben werden.
alphadot=[0; 0];
Vorinitialisierung
%% Darstellung der DGL
alphadot(1) = alpha(2) ;
erste Gleichung 1.Ordnung
alphadot(2) = -(g/l)sin(alpha(1)) ;
zweite Gleichung 1.Ordnung
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-Tiefpass

Beispiel 2: RC-Tiefpass

Die 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:
function [udot] = rckomb(t,u)
Funktion rckomb
Aufruf über:
[udot]= rckomb(t,u)

die vorliegende Funktion definiert die DGL einer RC-Kombination mit Spannungsquelle u1(t)
R = 10000 ;
C = 4.7*10e-6 ;
Definition der Konstanten:
R = Widerstand
C = Kapazität des Kondensators
udot=0 ;
Vorinitialisierung
udot = -(1/(RC))u + (1/(RC))u1(t) ;
die Gleichung 1. Ordnung
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:
function [sprng] = u1(t)
Funktion u1
Aufruf über:  [sprng]= u1(t)
sprng = t>=0 ;
Implementierung der Einheitssprungfunktion
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.
 Zurück
Vor