Lösung von Aufgabe 22

Vergleich der Zeiten und Schrittzahlen

Das Lösen der Differentialgleichung sollte keine Probleme bereiten. Am bequemsten definiert man sich eine Systemfunktion mit dem Skalenfaktor s als zusätzlichem Parameter:
- function dydt = f2d(t, y, s)
  % Schwingungs-DGL 2d, Skalenfaktor s für c2/b2 als Parameter
  m1 = 5; m2 = 0.5;
  c1 = 10; c2 = s*1;
  b1 = 0.1; b2 = s*0.1;
  F1 = 1;
  Om = 1.414;
  M = [m1, 0; 0, m2];
  C = [c1+c2,-c2;-c2,c2];
  B = [b1+b2,-b2;-b2,b2];
  Fhat = [F1; 0];
  dydt = [y(3:4); inv(M)*(Fhat*cos(Om*t) - B*y(3:4) - C*y(1:2))];
und verwendet bei gegebenem Wert von s die entsprechende Funktion "mit eingesetztem s"
- s = 1e2;
  f = @(t,y) f2d(t,y,s)
Die Rechenzeit bestimmt man mit
- tic; [t, y] = odeXY(...); toc
Die Zahl der Schritte ist einfach die Länge des Vektors t (oder y)
- nSteps = length(t)

Damit erhält man folgende Ergebnisse:

Skala	Solver	Zeit	nSteps
1.00e+00	ode45	0.40	57
1.00e+01	ode45	0.05	61
1.00e+02	ode45	0.01	69
1.00e+03	ode45	0.09	521
1.00e+04	ode45	0.80	5289
1.00e+05	ode45	7.96	53009
1.00e+06	ode45	111.61	530289
1.00e+00	ode113	0.20	27
1.00e+01	ode113	0.01	30
1.00e+02	ode113	0.02	46
1.00e+03	ode113	0.09	286
1.00e+04	ode113	0.92	2793
1.00e+05	ode113	9.14	27590
1.00e+06	ode113	108.33	275570
1.00e+00	ode15s	0.38	33
1.00e+01	ode15s	0.10	39
1.00e+02	ode15s	0.02	45
1.00e+03	ode15s	0.02	42
1.00e+04	ode15s	0.02	43
1.00e+05	ode15s	0.02	38
1.00e+06	ode15s	0.02	36

Die Solver ode45 und ode113 verhalten sich grundsätzlich gleich: Mit steiferer Feder nimmt die Zahl der Schritte und die Rechenzeit sehr stark zu. Dabei braucht ode113 im Schnitt nur etwa halb so viele Schritte, aber die gleiche Rechenzeit.
ode15s dagegen verhält sich wesentlich anders: Hier ist die Zahl der Schritte nahezu unabhängig von der Federsteifigkeit.

Vorwärts- und Rückwärtsrechnen

Um die tatsächliche Genauigkeit der verschiedenen Solver zu bestimmen, rechnet man zunächst bis zu t_e = 2 vorwärts:
- [t1, y1] = odeXY(@(t,y) f2d(t,y,s), [0 2], [0 0 0 0];
Der Wert des Zustandsvektors y1 am Endzeitpunkt t_e ist y1(end,:). Diesen wählt man als Anfangswert und gibt dem Solver t_e als Startzeit, 0 als Endzeit an:
- [t2, y2] = odeXY(@(t,y) f2d(t,y,s), [2 0], y1(end,:));
Bei idealer Genauigkeit sollte der letzte Wert von y2 wieder mit dem Anfangswert [0 0 0 0] übereinstimmen. Insbesondere ist der erhaltene Gesamtfehler für x₁
- d = abs(y2(end,1))

Für die verschiedenen Werte der Skala s und Solver erhält man

Für kleines s klappt alles, aber schon bei s = 100 ist der Fehler riesig. Für noch größere Werte explodiert die Lösung und überschreitet den Zahlenbereich. Auch der Solver ode15s verhält sich prinzipiell ähnlich, bekommt aber für sehr große s-Werte plötzlich wieder sinnvolle Resultate. Um das seltsame Verhalten zu verstehen, werden seine Ergebnisse (hin und zurück gerechnet) geplottet:
Für s = 10 ist alles in Ordnung, beide Lösungen liegen im Plot übereinander. Bei s = 100 wächst die Rücklösung (von rechts nach links gelesen!) plötzlich stark an, bei s = 1000 explodiert sie völlig. Versuche zeigen, dass auch eine erhöhte Genauigkeitsschranke (RelTol oder AbsTol) für die Solver nichts am prinzipiellen Verhalten ändern.
Die Ursache liegt in der Struktur der Differentialgleichung: Bei normaler Rechnung (t wächst) bewirkt der Reibungsterm B eine exponentielle Dämpfung, es geht Energie verloren und das System würde ohne äußere Anregung bald zur Ruhe kommen. In umgekehrter Richtung dagegen (t → -t) wird durch diesen Term die Reibungsenergie dem System wieder zugeführt. Bei idealer Genauigkeit sollte sich trotzdem der Anfangswert ergeben; leichte Abweichungen aber werden durch diese Zufuhr stark angefacht, umso mehr, je größer die Reibung ist. Dieses Problem lässt sich auch durch höhere Genauigkeit nur sehr begrenzt beherrschen, die ansteigende e-Funktion bringt selbst kleinste numerische Abweichungen schnell zur Explosion.

Alle Berechnungen können einfach, aber mühsam durch mehrfache Wiederholung ähnlicher Kommandos durchgeführt werden. Das Matlab-Skript ex22.m zeigt, wie man alle Aufgaben geschickt automatisieren kann.