Lösungen zu Übungsblatt 8

1. In den Übungen 6 hatten Sie ein Modell zur Lösung eines Anfangswertproblems mit einem System linearer Differentialgleichungen konstruiert:

   

   und u.a. das Anfangswertproblem:

      y1'(t) +    y1(t) = 0;               y1(0) = 1
      y2'(t) - 24*y1(t) + 25*y2(t) = 0;    y2(0) = 1

   simuliert.

   Setzen Sie in diesem Modell die folgenden Simulations-Parameter-Werte:

   • stop time=50
   • relative tolerance=1.0e-8
   • absolute tolerance=1.0e-8
   • max step size=20
   und das Abbruchkriterium im Subsystem 'Abbruch' auf 1.0e-10.

   Simulieren Sie das Problem mit allen Ihnen verfügbaren Integratoren und halten Sie fest, wann die Simulation jeweils beendet wird.

   Hinweis:
   Die allgemeine Lösung des Anfangswertproblems ist:

      y1(t) = a*exp(-t),
      y2(t) = b*exp(-t) + c*exp(-25*t)

   Durch die Vorgabe der speziellen Anfangswerte [1;1] ergeben sich a=b=1 und c=0. Das vorgegebene spezielle Anfangswertproblem hat daher die Lösung:

      y1(t) = exp(-t),
      y2(t) = exp(-t)

   Die Komponente exp(-25*t) taucht in der exakten Lösung nicht auf, kann sich aber infolge von Rundungsfehlern in der numerischen Lösung bemerkbar machen und dort insbesondere die Wahl der Schrittweiten beeinflussen.

   Wir verwenden zur Simulation das folgende Modell:

   

   und erhalten für die Abbruchzeit die folgenden Werte:

   Dormand-Prince	50.0
   Bogacki-Shampine	50.0
   Adams	27.11
   ode15s	21.06
   Rosenbrock	23.74
   Trapez	22.01
   TR-BDF	22.17

   
   Es gilt: exp(-18.42)=1.0e-8 und exp(-23.02)=1.0e-10. Bei exakter Integration müsste diese also bei 23.02 abgebrochen werden. Warum ist nicht zu erwarten, dass die Simulation immer bei 23.02 abgebrochen wird?

   Die numerische Integration erfolgt in zwei Bereichen (s.a. Aufgabe 2). Im ersten Bereich muss die numerische Lösung sich an der exakten Lösung orientieren. Es werden nur Abweichungen toleriert, die nicht über 1.0e-8*y_i(t) hinausgehen.
   Sobald jedoch abs(y_i(t)) < absolute tolerance = 1.0e-8 ist, wird jeder Integrationsschritt akzeptiert, der innerhalb deieser Schranken verbleibt. Die Kopplung der numerischen Lösung an die exakte Lösung wird also aufgegeben. Es ist deshalb nicht zwingend erforderlich, dass eine numerische Lösung überhaupt kleiner als 1.0e-10 wird.

   Auffallend ist, dass die steifen Integratoren die exakte Abbruchschranke wesentlich besser treffen als die nicht-steifen Integratoren. Die Ergebnisse hängen allerdings von der Rechnerarithmetik, der Simulink-Version und vielen anderen Faktoren ab und lassen sich nicht in jeder Umgebung so wie hier angegeben reproduzieren.

2. Entfernen Sie im Modell von Aufgabe 1 das Subsystem zum Abbruch der Simulation und ändern Sie die folgenden Simulations-Parameter:
   • stop time=100
   • max step size=50
   An die Integratoren werden dann, ein wenig vereinfacht formuliert, die folgenden Bedingungen gestellt:
   • von 0 bis etwa 18.42:
     die numerische Lösung muss in einem Schlauch um die exakte Lösung y(t)=exp(-t) bleiben, die Breite des Schlauches ist: 1.0e-8*exp(-t).
   • von 18.42 bis zum Ende:
     die numerische Lösung muss in einem Schlauch der Breite 1.0e-8 um 0 herim bleiben. Beachten Sie, dass diese Bedingung unabhängig von der Lösung y(t) ist.
   Benutzen Sie alle Ihnen zugänglichen Integratoren zur Behandlung des Problems. Welche Integratoren benutzen im zweiten Bereich in etwa dieselben Schrittweiten und welche Integratoren benutzen grössere Schrittweiten?

   Zur Simulation verwenden wir das folgende Modell:

   

   und erhalten damit die folgenden Schrittweiten:

   Dormand-Prince	Bogacki-Shampine
   Adams	ode15s
   Rosenbrock	Trapezregel

   Wie man sieht, benutzen die nicht-steifen Integratoren in beiden Bereichen in etwa dieselben Schrittweiten, während die steifen Integratoren ihre Schrittweiten im zweiten Bereich deutlich erhöhen. Beachten Sie, dass das Wachstum der Schrittweiten am Ende des Simulationsintervalles dadurch eingeschränkt wird, dass der Zielpunkt 100 erreicht werden muss.

3. Lassen Sie bei dem Modell von Aufgabe 2 jeweils anzeigen, ob die L2-norm der Lösung kleinergleich der absoluten Toleranz ist. Welche Integratoren erfüllen die Bedingung L2-Norm(y(t)) bis zum Simulationsende, nachdem sie sie zum ersten mal erfüllt haben? Welche Integratoren haben Probleme diese Bedingung einzuhalten?

   Hinweise:

   • die L2-Norm eines zweidimensionalen Wertes [f1,f2] wird gegeben durch:

        sqrt(f1^2+f2^2)

   • in Matlab wird der logische Wert 'wahr' durch 1, der logische Wert 'nein' durch 0 dargestellt. Sie können also das Ergebnis eines 'Relational Operator-Blockes' direkt in einen Scope-Block schicken.

   Wir benutzen das Modell:

   

   und erhalten damit die folgenden Ergebnisse:

   Dormand-Prince	Bogacki-Shampine
   Adams	ode15s

   Wie man sieht, haben es die nicht-steifen Integratoren schwer, die numerische Lösung in dem durch den Wert zu 'absolute tolerance' vorgegebenen Schlauch zu halten, alle drei Integratoren liefern Lösungen, die diesen Schlauch gelegentlich oder auch häufig wieder verlassen. Die steifen Integratoren haben dieses Problem nicht. Wir haben hier nur das Ergebnis zu ode15d wiedergegeben, die Ergebnisse der anderen steifen Integratoren sind von diesem nicht zu unterscheiden.

4. Die numerisch ermittelten Lösungskurven zum Anfangswertproblem aus Aufgabe 1 weisen bei der üblichen Darstellung mittels eines Scope-Blockes keine erkennbaren Unterschiede in Abhängigkeit vom benutzten Integrator auf. Um solche Unterschiede sichtbar zu machen, schalten Sie die Scope-Ausgabe erst dann ein, wenn die L2-Norm der Lösungskurven zum ersten mal kleiner als 'absolute tolerance', also kleiner als 1.0e-8 wird.

   Lösen Sie das Problem mit steifen und nicht-steifen Integratoren.

   Hinweise:

   • das Signal zum Scope-Block können Sie durch ein enabled Subsystem laufen lassen, das durch einen Block gesteuert wird, welcher hochzählt, wie oft die beobachtete L2-Norm kleiner als 'absolute tolerance' war.

   Im nachfolgenden Modell wird das Modell aus Aufgabe 3 durch einen Zähler ergänzt, der festhält, wie oft bisher schon die Bedingung 'L2-Norm < abs. tol' eingetreten ist. Dieser Zähler wird mit 0 initialisiert und steuert die Ausgabe der Ergebnisse in die 'Lupe'. Da die Steuerung über ein 'enabled subsystem' erfolgt, das lediglich die Ergebnisausgabe durchlässt oder nicht, wird die Ausgabe eingeschaltet, sobald die L2-Norm zum ersten mal kleiner als die absolute Toleranz wird.

   

   Das Bild im 'normal'-Scope sieht bei allen Integratoren gleich aus:

   

   In der 'Lupe' erhalten wir dagegen:

   Dormand-Prince	Bogacki-Shampine
   Adams	ode15s
   Rosenbrock	Trapezregel

   Man sieht, dass die nicht-steifen Integratoren die Komponente y₂ in Schritten berechnen, deren Richtung mit der der exakten Lösung nicht übereinstimmt. Beachten Sie allerdings die unterschiedliche Skalierung der Achsen bei der Bewertung der Richtungen. Die Bedingung, dass die Lösungskomponenten im Schlauch um 0 mit der Breite 'abs. tol.' bleiben sollen, erfordert daher sehr kleine Schritte.

   Die Lösungskomponente y₁ ist teilweise am Anfang des Ausschnittes zu erkennen, bei ihrer Berechnung treten obige Probleme offensichtlich nicht auf.

   Bemerkung: Je nach Rechnerumgebung kann das Dormand-Prince Verfahren auch Lösungen liefern, die in etwa spiegelbildlich zur x-Achse in negativen Teil des 'abs.tol'-Schlauches liegen.

   Bei der Trapezregel fällt auf, dass sie im Schlussbereich eine schwingende Lösungskurve mit alternierendem Vorzeichen liefert und die Schwingungen nur langsam herausgedämpft werden, eine Folge der fehlenden L-Stabilität der Trapezregel. Zur Verdeutlichung zeigen wir diesen Bereich noch stärker vergrößert und mit längerer Laufzeit.

ode15s	Trapezregel