K. Taubert, W. Wiedl: Matlab Differentialgleichungen

1. Integratoren für den allgemeinen Fall.

die im einfachsten Falle alle durch:

   [t,y] = Integrator(M-Function, t-Intervall, Anfangswerte);

• t, ein Spaltenvektor der Länge n mit den t-Koordinaten der Ausgabepunkte, t und die seine Länge werden vom Integrator besetzt,
• y, ein n*m Feld mit den y-Koordinaten der Lösungsvektoren y(:,i), m ist die Anzahl der Differentialgleichungen im System und wird mit Hilfe der M-Function im 1. Argument festgelegt. Die Matrix y wird vom Integrator besetzt.
• Integrator, der Name eines der vorhandenen Integratoren. Zur Auswahl für den allgemeinen Fall, für den keine speziellen Integratoren erforderlich sind, stehen:
  • ode23, Bogacki-Shampine Integrator, ein expliziter Runge-Kutta Integrator niedriger (2,3) Ordnung. Einschrittverfahren, benötigt 3 Auswertungen des Differentialgleichungssystems dy = F(t,y) pro Integrationsschritt.
    Dieser Integrator wird in der Matlab-Dokumentation als erste Wahl empfohlen, wenn nur geringe Genauigkeit erforderlich ist.
    Literatur: Bogacki,P. and Shampine,L.F.: "A 3(2) pair of Runge-Kutta formulas.", Appl.Math.Letters,2,1989,pp 1-9)
  • ode45, Dormand-Prince Integrator, ein expliziter Runge-Kutta Integrator mittlerer (4,5) Ordnung. Einschrittverfahren,
    Dieser Integrator wird in der Matlab-Dokumentation als erste Wahl empfohlen, wenn keine Informationen über die Eigenschaften des Systems vorliegt.
  • ode113, Adams-Bashforth-Moulton PECE-Integrator variabler Ordnung. Mehrschrittverfahren, benötigt 2 Auswertungen des Differentialgleichungssystems dy = F(t,y) pro Integrationsschritt.
    Dieser Integrator wird in der Matlab-Dokumentation als erste Wahl bei Problemen empfohlen, die keine speziellen Integratoren und hohe Genauigkeit erfordern oder wenn die Berechnung des M-Files, das die Ableitungen liefert, sehr aufwendig ist.
  ode23 und ode45 sind explizite Verfahren, d.h. der Ergebnisvektor y_n+1 eines Integrationsschrittes von t_n nach t_n+1 wird durch arithmetische Ausdrücke bestimmt, die von y_n+1 unabhängig sind.
  Als Alternative zu den expliziten Verfahren, gibt es implizite Verfahren, bei denen der Ergebnisvektor y_n+1 durch arithmetische Ausdrücke bestimmt wird, die von der Ergebnisvektor y_n+1 abhängen. Im Normalfall bedeutet dies, daß ein Gleichungssystem gelöst werden muß. Ein Integrationsschritt mit einem impliziten Verfahren erfordert daher meist wesentlich mehr Aufwand als ein Integrationsschritt mit einem expliziten Verfahren.
  ode23 und ode45 sind Einschrittverfahren, im Integrationsschritt von t_n nach t_n+1 wird von den bisherigen Ergebnisvektoren y_i nur der Vektor y_n verwendet, aber keiner der weiter zurückliegenden Vektoren y_n-1, y_n-1.
• M-Function, der Name einer M-Function, die im einfachsten Fall nur das System der gewöhnlichen Differentialgleichungen 1. Ordnung beschreibt. Diese Funktion wird vom Integrator und i.d.R. nicht vom Benutzerprogramm aufgerufen. Der Funktionskopf hat im einfachsten Fall die Form:

  function dy = Funktionsname(t,y)

  Die Argumente t und y müssen auch dann angegeben werden, wenn sie zur Berechnung der Ableitungen nicht benötigt werden. Diese beiden Argumente sind Spaltenvektoren. Der Funktionswert, hier dy genannt, muß die Ableitungen der m Differentialgleichungen liefern.
• t-Intervall, ein Zeilen- oder Spalten-Vektor der Länge 2, der das Integrationsintervall festlegt.
• Anfangswerte als Zeilen- oder Spalten-Vektor der Länge m, der die Anfangswerte des Systems festlegt.

1. die Formulierung des Problems in der angegebenen Form,
2. die Wahl eines geeigneten Integrators.

1.1 Formulierung des Problems.

Beispiel 1.1

Ein Standardbeispiel, das wir später auch für spezielle Integratoren verwenden werden, ist die 'van der Pol-Differentialgleichung'. Dieses System 2. Ordnung muß in einer M-Function als System 1. Ordnung beschrieben werden. Wir setzen:

   y1'(t) = y2(t)

und erhalten für µ=2 damit:

   y1' = y2
   y2' = 2*(1-y12) * y2 - y1;

Bedingt durch die Forderung, daß der Funktionswert als Spaltenmatrix geliefert werden muß, formulieren wir die M-Function also:

   function dy = VDP(t,y)
   dy(1,1) = y(2);
   dy(2,1) = 2*(1 - y(1)^2)*y(2) - y(1);

Der Aufruf des Integrators und die

   [t,y] = ode23('VDP',[0 20],[2; 0]);
   plot(t,y)

Beispiel 1.2

   function dy=RIGID(t,y)
   dy = [ y(2)*y(3);
         -y(1)*y(3);
         -0.51*y(1)*y(2)];

   [t,y] = ode45('RIGID',[0;12],[0 1 1]);

In den neuesten Versionen von Matlab (seit Release 11) ist die möglich, die obigen Integratoren auch für Probleme zu nutzen, die in der Form:

M(t) y'(t) = F(t,y);  y(0) = y0

Wenn ein Problem mit Massenmatrix formuliert ist, muß dies den Integratoren mitgeteilt werden. Man benutzt dazu beim Aufruf der Integratoren einen 4. Stellungsparameter, dessen Werte man mit Hilfe einer zusätzlichen odeset-Anweisung besetzt:

   opt = odeset(...Eigenschaften festlegen...);
   [t,y] = integrator('M-File',Integrationsintervall,Anfangswerte,opt); 
   plot(t,y)

   'none' | 'M' | 'M(t)' | 'M(t,y)' 

Die Massenmatrix selbst fragt der Integrator vom M-File ab, das zu diesem Zwecke mit einem dritten Stellungsparameter formuliert werden m&uszlig;:

   function dy=name(t,y,flag)

Beispiel 1.3

   function wert=DRIGID(t,y,flag)
   switch flag
     case ''
       wert = [ y(2)*y(3);
             -y(1)*y(3);
             -0.51*y(1)*y(2)];
     case 'mass'
       wert = [ 1 0; 0 1];
   end

Der Aufruf erfolgt jetzt durch:

   opt = odeset('Mass','M');
   [t,y] = ode113('DRIGID',[0;12],[0 1 1]);
   plot(t,y);

Bemerkung: In früheren Matlab-Versionen wurden die Eigenschaften der Massenmatrix anders beschrieben. Die hier gegebene Darstellung ist ab Version 5.2 gültig.

1.2 Fehlerkontrolle.

abs(LokalerFehleri) <= max(RelTol*abs(yi),AbsTol(i))

Beispiel 1.4:

y'(t) = µ ( y(t) - sin(t) ) + cos(t),   y(0) = 1

y(t) = sin(t) + exp( µt )

y'(t) = µ ( y(t) - sin(t) ) + cos(t),   y(0) = 1+d

y(t) = sin(t) + (1+d) exp( µt )

1.3 Schrittweitenkontrolle.

Die Vorgabe eines Wertes für die maximale Schrittweite erfolgt über den Wert von MaxStep und ist immer eine positive skalare Zahl. Sinnvoll ist eine Vorgabe z.B. dann, wenn man die Periode einer Lösung kennt und sicher sein will, daß die Schwingung

1.5 Integrator Ausgabe

   opt = odeset('Refine',1,'AbsTol',1.0e-9,'RelTol',1.0e-6);
   [t1,y1] = ode45('VDP',[0 20],[2 0],opt);
   [t2,y2] = ode23('VDP',[0 20],[2 0],opt);
   subplot(2,1,1)
   plot(t2,y2)
   title('Lösung der VDP-Gleichung mit µ=2')
   subplot(2,1,2)
   n1 = length(t1);
   h1 = t1(2:n1) - t1(1:n1-1);
   n2 = length(t2)
   h2 = t2(2:n2) -t2(1:n2-1);
   plot(t1(1:n1-1),h1,t2(1:n2-1),h2)
   xlabel('ode45/ode23-Schrittweite bei Fehlertoleranz 1.0e-9/1.0e-6')

Listen mit Ausgabezeitpunkten
Das Integrationsintervall haben wir bisher immer in der Form [Anfangswert Endwert] angegeben. Die Ausgabewerte werden dann gesteuert durch:

• die Integrationsschritte
• den Refine-Wert, der mehrere Ausgabepunkte pro Integrationsschritt erzwingen kann.

   [t,y] = ode45('VDP',pi*(0:-1:-4),[2; 0]);

Statistiken zur Arbeit der Integrators
Mit Hilfe der Stats Eigenschaft kann jeder Matlab Integrator dazu veranlaßt; werden, genauere Angaben zu seiner Arbeit auszugeben, dies sind im Einzelnen:

1. die Anzahl der erfolgreichen Schritte,
2. die Anzahl der Fehlversuche,
3. wie oft dy = F(t,y) ausgewertet wurde,
4. wie oft die Jacobi-Matrix dF/dy bzw. eine Näherung dazu aufgestellt wurde,
5. die Anzahl der LU-Zerlegungen von dF/dy,
6. wie oft ein lineares Gleichungssystem gelöst wurde.

Beispiel 1.6:

   opt   = odeset('Mass','M','Stats','on');
   [t,y,S] = ode45('DYRIGID2',Intv,Ynull,opt);

    19 successful steps
     2 failed attempts
   127 function evaluations

Eigene Ausgabefunktionen
Matlab bietet die Möglichkeit, nach jedem Integrationssschritt eine eigene Ausgabefunktion aufzurufen, der Name dieser Funktion, genauer, der Name des M-Files, muß als Wert der Eigenschaft OutputFcn angegeben werden.

   opt = odeset('OutputFcn','Eigenw','Refine',1);   
   [t,y] = ode45('VDP',[0 20],[2 0],opt);

Der Rückgabewert einer Ausgabefunktion muß ein Skalar mit dem Wert 0 oder 1 sein. Er bestimmt, ob die Integration abgebrochen wird (1) oder fortgesetzt wird (0).

Der Integrator ruft die Ausgabefunktion auf:

• zu Beginn der Integration mit den drei Argumenten t, y und 'init',
• nach jedem Integrationsschritt mit den zwei Argumenten t und y,
• am Ende der Integration mit den drei Argumenten [], [] und 'done'.

Beispiel 1.7:

   function dy = DOUTFUN(t,y,flag)
   global my;
   switch flag
     case ''
       dy = [ y(2);
              my*(1 - y(1)^2)*y(2) - y(1) ];
     case 'Jacobian'
       dy = [ 0.0                     1.0;
              -my*2*y(1)*y(2)-1       my*(1-y(1).^2) ];
   end

   function out = OUTOUT2(t,y,flag)
   if nargin == 2 | flag == 'init'
     eiw = eig(DOUTFUN(t,y,'Jacobian'))';
     if imag(eiw(1)) ~= 0
       out = 1;
     else
       out = 0;
     end
   end

   global my eiw;
   my = 5.0;
   opt = odeset('Refine',1,'OutputFcn','Outout2', ...
                'AbsTol',1.0e-9,'RelTol',1.0e-6);
   [t,y] =ode45('VDP',[0 20],[2 0],opt);
   plot(t,y)
   title('Lösung der van der Pol-Gleichung für µ=5')
   xlabel('t-Intervall mit realen Eigenwerten')

Beispiel 1.8:

   function out = OUTOUTF(t,y,flag)
   global my eiw;
   if nargin == 2 | flag == 'init'
     [n1,n2] = size(eiw);
     eiw(n1+1,:) = eig(DOUTFUN(t,y,'Jacobian'))';
     out = 0;
   else
     out = 1;
   end

   clear all
   global my eiw;
   my = 5.0;
   opt = odeset('Refine',1,'OutputFcn','Outoutf', ...
                'AbsTol',1.0e-9,'RelTol',1.0e-6);
   [t,y] =ode45('VDP',[0 20],[2 0],opt);
   subplot(3,1,1)
   plot(t,y)
   title('Lösung der van der Pol-Gleichung mit µ=5');
   subplot(3,1,2)
   plot(t,real(eiw))
   subplot(3,1,3)
   plot(t,imag(eiw))
   xlabel('Eigenwerte der Jacobi-Matrix')