¨ 170 005 Ubungen zu Numerische Methoden I ¨ Neunte Ubungseinheit

9., 10. und 11. Juni 2008

¨ Inhalt der neunten Ubungseinheit: • Gew¨ ohnliche Differentialgleichungen • KNW-Musterbeispiele

9 9.1

Gew¨ohnliche Differentialgleichungen Einschrittverfahren Bei der numerischen L¨ osung einer DG bestimmt man ausgehend von den Anfangsbedingungen eine Folge von Wertepaaren (x0 ,y0 ), (x1 ,y1 ), (x2 ,y2 ), . . ., die den Verlauf der gesuchten Funktion y = y(x) ann¨ ahern sollen. Schema: W¨ ahle Schrittweite h und maximale Schrittzahl N ; setze x0 und y0 laut Anfangsbedingung; f¨ ur i = 0,1, . . . ,N xi+1 = xi + h ; yi+1 = yi + hF (xi ,yi ,h) . Dabei bedeutet F (x,y,h) beim Euler-Verfahren F (x,y,h) = f (x,y), beim modifizierten Euler-Verfahren F (x,y,h) = f



 h h x + , y + f (x,y) , 2 2

beim Verfahren von Heun F (x,y,h) =

1 (k1 + k2 ) 2

mit k1 = f (x,y) k2 = f (x + h,y + hf (x,y)).

9.2

MATLAB-Funktionen zum L¨osen gew¨ohnlicher Differentialgleichungen Selbstverst¨ andlich ist es m¨ oglich, diverse Verfahren zum L¨ osen von Anfangswertproblemen selbst in MATLAB zu implementieren. MATLAB bietet aber auch einige eingebaute Solver f¨ ur Anfangswertprobleme an. Sie arbeiten im Prinzip wie die obigen Einschrittverfahren, aber • die Fortschreitrichtung wird aus mehreren, optimal gew¨ ahlten Zwischenpunkten bestimmt, 1

• der (globale) Fehler wird abgesch¨ atzt, • die Schrittweite wird automatisch angepasst. In der Version 6.5 sind sieben Solver verf¨ ugbar: ode23, ode23s, ode23t, ode23tb, ode45, ode113, ode15s. Eine genauere Beschreibung k¨ onnen Sie in der MALAB-Hilfe nachlesen. Die Funktion ode45 ist das Standardverfahren, welches man f¨ ur ein neues“ Anfangswert” problem zuerst probieren wird. Sie ist die Implementierung eines expliziten EinschrittRunge-Kutta-Verfahrens (Fehlerordnung 5). Es folgt ein kleines Beispiel, das die Verwendung von ode45 demonstriert. Angenommen, Sie wollen das Anfangswertproblem y ′ = g(x,y) = 2x(1 + y 2 ),

y(0) = 0,

im Intervall [0,1] mit diesem Solver l¨ osen. Sie speichern daf¨ ur die Funktion g in einer Datei g.m ab function ydot=g(x,y) ydot=2*x*(1+y^2); und fassen dann die L¨ osung des Anfangswertproblems in einer Scriptdatei (test.m) zusammen: y0=0.0; xspan=[0, 1]; ode45(@g,xspan,y0) So einfach ist das! Sie erhalten eine graphische Darstellung der L¨ osung. Wenn Sie die L¨ osung numerisch f¨ ur bestimmte x-Werte auswerten wollen, schreiben Sie beispielsweise sol=ode45(@g,xspan,y0); x = [0 0.25 0.5 0.75 1]; y = deval(sol,x)

9.3

Differentialgleichungen h¨oherer Ordnung, Systeme von DG erster Ordnung Differentialgleichungen h¨ oherer Ordnung lassen sich als Systeme von Differentialgleichungen erster Ordnung schreiben. Als Beispiel werden hier die Bewegungsgleichungen eines Systems aus drei durch Federn gekoppelten Massen diskutiert. Die Positionen x1 ,x2 ,x3 der drei Massen sind als Funktion der Zeit t durch ein System von Differentialgleichungen 2. Ordnung bestimmt. x ¨1 = −2x1 + x2 x ¨2 = x1 − 2x2 + x3 x ¨3 = x2 − 2x3 Anfangswerte sind x1 (0) = x2 (0) = x3 (0) = 1 , x˙ 1 = x˙ 2 = x˙ 3 = 0 . 2

Zur Schreibweise: In der Physik heißt die durch Differentialgleichungen bestimmte Funktion h¨ aufig nicht y = y(x), sondern x = x(t). Dabei ist x eine Ortskoordinate und t die Zeit. Der Ausdruck x˙ bedeutet die erste, x ¨ die zweite Ableitung von x nach t. Wir transformieren durch Einf¨ uhren folgender Hilfsfunktionen auf ein System von sechs Differentialgleichungen 1. Ordnung. y1 (t) = x1 (t), y2 (t) = x2 (t), y3 (t) = x3 (t), y4 (t) = x˙ 1 (t), y5 (t) = x˙ 2 (t), y6 (t) = x˙ 3 (t) . Das ¨ aquivalente System 1. Ordnung lautet dann y˙ 1 = y4 y˙ 2 = y5 y˙ 3 = y6 y˙ 4 = −2y1 + y2 y˙ 5 = y1 − 2y2 + y3 y˙ 6 = y2 − 2y3 Der Standard-L¨ oser ode45 braucht eine Funktionsdatei dreimassen.m, welche f¨ ur einen Eingabevektor y die sechs Ableitungen als Vektor ydot berechnet. function ydot=dreimassen(t,y) ydot=[y(4); y(5); y(6); -2*y(1)+y(2); y(1)-2*y(2)+y(3); y(2)-2*y(3)] Damit ben¨ otigen wir nur noch eine Scriptdatei loesung1.m in der die Anfangsbedingungen f¨ ur die gesuchte vektorwertige Funktion y sowie deren Ableitung spezifiziert werden und das betrachtete Intervall definiert wird. Danach k¨ onnen wir schon den Solver ode45 aufrufen, der hier mit der Option, nur die ersten drei Komponenten des Systems auszugeben, gestartet wird. Unser Plot zeigt dann nur die Kurven f¨ ur die Auslenkungen osungen von y1 und y3 zusammeny1 (t) = x1 (t), y2 (t) = x2 (t), y3 (t) = x3 (t), wobei die L¨ fallen. y0=[1;1;1;0;0;0]; tspan=[0 2]; options=odeset(’OutputSel’,[1,2,3]); ode45(@dreimassen,tspan,y0,options) Versuchen Sie auch den Befehl ode45(@dreimassen,tspan,y0,options) durch [y,t]=ode45(@dreimassen,tspan,y0,options) zu ersetzen. Was ¨ andert sich an der Ausgabe?

3

Allgemeines Schema DG h¨oherer Ordnung −→ System von DG 1. Ordnung Umformen von y (d) = f (x,y,y ′ , . . . ,y (d−1) ) Man setzt y1 = y, y2 = y ′ , . . . , yd = y (d−1) und schreibt

y1′ = y2 y2′ = y3 .. . ′ yd−1 = yd

yd′ = f (x,y1 , . . . ,yd )

9.4

¨ Ubungsaufgaben, KNW-Musterbeispiele Aufgabe 65

Kepler-Ellipsen

Planeten durchlaufen elliptische Bahnen; bei jeder Ellipse steht die Sonne in einem Brennpunkt (1. Keplersches Gesetz). Gesucht sind die Koordinaten eines Planeten in der x − yEbene als Funktion der Zeit: x = x(t), y = y(t). Die Bewegungsgleichungen lauten: p x2 + y 2 r = x x ¨ = − 3 r y y¨ = − 3 r (Gravitationskonstante und Sonnenmasse hier der Einfachheit halber = 1 gesetzt). Anfangsbedingungen: x(0) = 10, x(0) ˙ = 0, y(0) = 0, y(0) ˙ = 0,2 Mit diesen Anfangswerten dauert ein Umlauf T = 31,25π ≈ 98,1747704. Bei exakter Rechnung w¨ urde sich die Umlaufbahn schließen, das heißt, die Werte zu dieser Zeit T w¨ aren gleich den Startwerten. Wie groß ist der Fehler des Standardl¨ osers? Stellen Sie die elliptische Umlaufsbahn in der (x,y)-Ebene dar. Aufgabe 66

Populationsdynamik

Ein mathematisches Modell, das Entwicklung des Bestandes zweier Arten ( R¨ auber“ und ” Beutetiere“) in vereinfachter Weise darstellt, wird durch das System nichtlinearer Diffe” rentialgleichungen y1′ (t) = k1 y1 (t) − k2 y1 (t)y2 (t) y2′ (t) = k3 y1 (t)y2 (t) − k4 y2 (t) ausgedr¨ uckt. Dabei stellen y1 (t) und y2 (t) den zeitlichen Bestand der Beutetiere beziehungsweise der R¨ auber dar. 4

L¨ osen Sie dieses System f¨ ur 0 ≤ t ≤ 4 unter der Annahme, daß der Anfangsbestand an Beutetieren 1000, der der R¨ auber 200 und die Konstanten k1 = 3, k2 = 0,002, k3 = 0,0006 und k4 = 0,5 betragen. Methode: ode45. Um die Ergebnisse zu veranschaulichen, stellen Sie die L¨ osung, also die Anzahl der beiden Populationen versus Zeit, graphisch dar. Wiederholen Sie die Rechnung f¨ ur einen Anfangsbestand von 2000 Beutetieren und 1500 R¨ aubern. Aufgabe 67

Schmetterlingseffekt

Der Meteorologe Edward N. Lorenz formulierte 1963 ein System von drei gekoppelten, nichtlinearen gew¨ ohnlichen Differentialgleichungen zur Modellierung der Erdatmosph¨ are zum Zweck einer Langzeitvorhersage. dx/dt = a(y − x) dy/dt = x(b − z) − y dz/dt = xy − cz Die numerische L¨ osung des Systems zeigt bei bestimmten Parameterwerten chaotisches Verhalten. Die typische Parametereinstellung mit chaotischer L¨ osung lautet: a = 10,

b = 28,

c=

8 3

L¨ osen Sie dieses System f¨ ur 0 ≤ t ≤ 50, Anfangsbedingungen x = 20, y = 20, z = 40. Stellen Sie die L¨ osung als Kurve in Parameterform in einem xyz-Koordinatensystem (Befehl plot3) dar. 50 45 40 35 30 25 20 15 10 5 0 40 20 0 −20 −40

−20

−15

−10

−5

0

5

10

Die L¨ osung ist f¨ ur ihre typische Schmetterlings-Form bekannt.

5

15

20

Aufgabe 68

Datenanpassung

Zus¨ atzlich zu Differentialgleichungen k¨ onnen auch Beispiele zur Datenanpassung zum KNW kommen. Typische Aufgaben w¨ aren vom Typ Aufgabe 51 oder 52, oder die folgende: Die Zentralanstalt f¨ ur Meteorologie und Geodynamik gibt f¨ ur die Landeshauptst¨ adte folgende Werte der magnetischen Deklination an: Deklinationswerte der einzelnen Landeshauptst¨ adte (¨ ostliche Deklination bezogen auf Jahresmitte 2006) Stadt Wien (WIK) Eisenstadt St.P¨ olten Graz Linz Klagenfurt Salzburg Innsbruck Bregenz

Breite x 48.20 47.85 48.20 47.07 48.30 46.62 47.80 47.27 47.50

L¨ ange y 16.37 16.52 15.63 15.45 14.30 14.31 13.03 11.40 9.77

Deklination D 2.83 2.68 2.73 2.60 2.35 2.32 2.05 1.67 1.20

Finden Sie f¨ ur diese Daten eine bilineare Anpassung der Form z(x,y) = a1 + a2 x + a3 y + a4 xy Geben Sie die Differenz D − z(x,y) zwischen den tats¨ achlichen Deklinationswerten und den angepassten Werten an. Welcher Wert wird am schlechtesten approximiert? Aufgabe 69 F¨ ur die Blies (einen Nebenfluss der Saar) sollen die Hochwasserst¨ ande am Pegel Neunkirchen aus den Wasserst¨ anden des Pegels Ottweiler und des Pegels Hangard vorhergesagt werden. Es liegen die Daten der Scheitelwasserst¨ ande von 12 Winterhochw¨ assern aus den Jahren 1963–1971 vor: Wasserstand in cm Neunkirchen y Ottweiler x1 Hangard x2

172 93 120

309 193 258

302 187 255

283 174 238

443 291 317

298 184 246

319 205 265

419 260 304

361 212 292

267 169 242

337 216 272

230 144 191

Quelle: U. Maniak, Hydrologie und Wasserwirtschaft, Springer, 1988

Wir unterstellen den Daten das lineare Modell a0 + a1 x1 + a2 x2 = y. In diesem Ansatz sind a0 , a1 und a2 unbekannte Koeffizienten, die aus den zw¨ olf gegebenen Werte-Tripeln m¨ oglichst gut bestimmt werden sollen. Berechnen Sie die Koeffizienten des linearen Modells, und geben Sie auch den Differenzenvektor zwischen Modellvorhersage und Messwerten an. (Dieses Beispiel ist in den Erg¨ anzungen zum Skriptum Kapitel 5.4 ausf¨ uhrlich durchgerechnet auf der VorlesungsfolienInternetseite verf¨ ugbar.)

6