K. Taubert, W. Wiedl: Matlab Differentialgleichungen

2.6.1 Ruby Laser Oszillator.

   y1' = -y1*(alpha*y2 + beta) + gamma
   y2' = y2*(rho*y1 - sigma) + tau*(1+y1)

   y1(0) = -1,  y2(0) = 0

   alpha=1.5e-18, beta=2.5e-6, gamma=2.1e-6, rho=0.6, sigma=0.18, tau=0.016

   function w = DLOM(t,y,flag)
   % Laser Oscillator Model
   alpha = 1.5e-18;
   beta  = 2.5e-6;
   gamma = 2.1e-6;
   rho   = 0.6;
   sigma = 0.18;
   tau   = 0.016;
   switch flag
     case ''
       w = [ -y(1)*(alpha*y(2)+beta)+gamma;
              y(2)*(rho*y(1)-sigma)+tau*(1+y(1))];
     case 'jacobian'
       w  = [ -(alpha*y(2)+beta),    -y(1)*alpha ;
              y(2)*rho + tau,        (rho*y(1)-sigma) ];
   end

Szenario 1:
Physikalische Gründe lassen über dem Integrationsintervall den folgenden Verlauf erwarten. Dieser Verlauf weist deutlich zwei Phasen auf.
Im ersten Teil des Intervalls findet man einen recht einfachen Verlauf der Lösungskurven. Die Jacobi-Matrix hat am Startpunkt t=0, y=[-1 0] die Eigenwerte [-2.5e-6 -7.8e-1] und damit ein Steifigkeitsverhältnis von 3.0e+5. Viele Voraussetzungen für ein steifes Problem sind damit gegeben. Man wird diesen Teil z.B. mit ode15s integrieren.
Im zweiten Teil des Intervalls erfordert die Integration des Systems vermutlich von jedem Integrator erheblichen Aufwand. Steife Integratoren werden hier nicht von Vorteil sein. ode45 dürfte hier die erste Wahl sein.
Wir gehen hier davon aus, daß bekannt ist, daß der Wechsel vom steifen zum nicht-steifen Integrator bei 4.95e+5 zu erfolgen hat, dies ist ungefähr die Stelle, an der das System zu schwingen beginnt.

Der Matlab-Code mit dem Integrator-Aufruf lautet dann:

   opt = odeset('AbsTol',1.0e-8,'RelTol',1.0e-5,'Jacobian','on');
   tic;
   [ta,ya] = ode15s('DLOM',[0 4.95e+5],[-1 0],opt);
   n = length(ta);
   opt = odeset(opt,'Jacobian','off');
   [tb,yb] = ode45('DLOM',[ta(n) 1.0e6],ya(n,:),opt);
   toc
   t  = [ta;tb];
   ya(1,2) = ya(2,2);
   y1 = [ya(:,1);yb(:,1)];
   y2 = [ya(:,2);yb(:,2)];
   plot(t,y1,t,log(y2));

Das unterschiedliche Verhalten in den beiden Teilen des Integrationsintervalles macht sich auch bei den Eigenwerten der Jacobi-Matrix bemerkbar.
Der Realteil des einen Eigenwertes bleibt über dem ganzen Intervall bei etwa -2.5e-6, der Betrag des anderen Realteils sackt im ersten Teil des Integrationsintervalls von -7.8e-1 auf den Wert des ersten ab und schwingt dann um diesen Wert.
Im ersten Teil des Intervalls sind die Eigenwerte reell, im zweiten treten periodisch komplexe Eigenwerte auf, wobei Periode und Amplitude mit zunehmenden t immer kleiner werden. Dasselbe gilt auch für die Schwingungen von phi.

Insgesamt kann man für die Integration einigen Aufwand erwarten und es dürfte deshalb sinnvoll sein, für die beiden unterschiedlichen Phasen verschiedene Integratoren zu benutzen. Dies setzt allerdings voraus, daß man weiß, an welcher Stelle man von einem steifen zu einem nicht-steifen Integrator wechseln sollte.
Andererseits scheint das System nicht derart aufwendig zu sein, daß man das Problem nicht durchgehend mit ode45 integrieren könnte.

Durchlaufzeiten

	beide Intervalle	0 bis 4.95e5	4.95e5 bis 1.0e6
ode45	406.70	379.29	27.41
ode113	1138.70	1059.90	78.80
ode15s	831.23	1.85	829.38
ode15s/Jacobian	865.79	1.56	864.23
ode15s/BDF		1.93
ode15s/Jacobian/BDF	929.12	1.67
ode23s	(*)	21.07	(*)
ode23s/Jacobian		11.82

Die Durchlaufzeiten wurden auf einem 200 MHz Pentium mit 32 MB RAM und den Funktionen tic bzw. toc gemessen.
Man sieht, daß die Integration mit ode15s/Jacobian bis 4.95e+5 und danach mit ode45 bis 1.0e6 etwa 1.56 + 27.41 = 28.97 Sekunden dauert, wobei auch schnellere Durchlaufzeiten, etwa 27.8 Sekunden gemessen wurden.
Den besten Wert ohne Wechsel des Integrators liefert ode45 mit 406.70 Sekunden, das ist fast 15 mal mehr als die obigen 28.97 Sekunden.
Ganz ungeschickt wäre es, die erste, steife Phase mit ode113 und die zweite, nicht-steife Phase mit einem steifen Integrator zu rechnen, man käme so leicht an eine Laufzeit von 2000 Sekunden heran. Die Markierung (*) besagt, daß der betreffende Integrator im betreffenden Intervall Warnungen bezüglich der schlechten Kondition der Jacoby-Matrix ausgegeben hat.

Szenario 2:
Wenn man keine Vorstellungen über den Verlauf der Lösung und die Eigenschaften des Systems hat, wird man in einem ersten Versuch vermutlich ode45 als Integrator wählen und nach einiger Zeit die Integration abbrechen, da sie zu lange dauert. Man wird in einem solchen Falle auch eine Output Function verwenden, um etwas über das Verhalten der Integratoren bei diesem System zu erfahren.

   opt = odeset('AbsTol',1.0e-8,'RelTol',1.0e-5,'OutputFcn','odeprint');
   [t,y] = ode45('DLOM',[0 1.0e6],[-1 0],opt);

Zum Vergleich probieren wir ode15s/Jacobian mit diesem Problem. Als Integrationsintervall wählen wir 0 bis 100, als Output Funktion 'odeplot'. Die Graphikausgabe zeigt, daß ode15s in diesem Startintervall schon Schrittweiten der Länge 10 erreicht. Das vorliegende Problem könnte steif sein.
Eine Berechnung der Eigenwerte der Jacoby-Matrix zur Zeit t=0 ergibt -7.8e-1 und -2.5e-6, also ein Verhältnis von etwa 3e+5, so daß vieles für ein steifes Problem spricht.
Allerdings führt auch eine Integration mit ode15s/Jacobian nicht rasch zum Ziel, so daß man sich genauer mit dem Problem beschäftigen wird. Hat man ode15s/Jacobian mit einer Output Function, etwa 'odeprint' laufen lassen, so hat man gesehen,daß die Integration bis etwa 5.0e+5 problemlos ist.
Eine Möglichkeit wäre, daß sich die Eigenschaften des Problems im Laufe des Integrationsintervalls ändern, so daß man über einigen Teilintervallen vorteilhaft mit steifen Integratoren, über anderen dagegen mit nicht steifen Integratoren rechnen sollte.
Man kann nun bis 5.0e+5 oder 4.95e+5 einen steifen Integrator verwenden und danach die Integration mit ode45 fortsetzen, was wir in Szenario 1 getan haben.
Falls die Übergangsstelle von steifer zu nicht-steifer Integration aus den vorangegangenen Versuchen nicht klar ist, kann man versuchen, sie mit Hilfe der Eigenwerte zu ermitteln, was bei diesem System nicht allzu aufwendig ist. Wir verwenden dazu die Ereignissteuerung in den Matlab Integratoren.
Wir haben gesehen, daß die Jacobi-Matrix zum Zeitpunkt t=0 eine Steifigkeitskoeffizienten von etwa 3e+5 aufweist und vermuten, daß das Absinken dieses Koeffizienten dazu führt, daß das System nach einiger Zeit nicht mehr steif ist. Wir integrieren daher in einem ersten Schritt solange bis dieser Koeffizient kleiner 5.0 wir. Der betreffende Matlab-Code lautet:

   opt = odeset('AbsTol',1.0e-8,'RelTol',1.0e-5);
   Ival   = [ -1; 0 ];  
   opt = odeset('AbsTol',1.0e-8,'RelTol',1.0e-5, ...
                             'Jacobian','on','Events','on');
   [ta,ya,te] = ode15s('RubyODE',[0 1.0e+6],Ival,opt);
   n = length(ta);
   opt = odeset(opt,'Events','off');
   [tb,yb] = ode45('RubyODE',[ta(n) 1.0e+6],ya(n,:),opt);
   toc
   te

   function [w,ist,dir] = DLOM(t,y,flag)
% Laser Oscillator Model
   alpha = 1.5e-18;
   beta  = 2.5e-6;
   gamma = 2.1e-6;
   rho   = 0.6;
   sigma = 0.18;
   tau   = 0.016;
   switch flag
     case ''
       w  = [ -y(1)*(alpha*y(2)+beta)+gamma;
               y(2)*(rho*y(1)-sigma)+tau*(1+y(1))];
     case 'jacobian'
       w  = [ -(alpha*y(2)+beta),    -y(1)*alpha ;
              y(2)*rho + tau,        (rho*y(1)-sigma) ];
     case 'events'
       eiw  = eig([ -(alpha*y(2)+beta),    -y(1)*alpha ;
                    y(2)*rho + tau,        (rho*y(1)-sigma) ]);
      w  = min(eiw)/max(eiw)-5.0;
      ist= [1];        
      dir= [0]; 
   end