MATHEMATISCHES INSTITUT DER UNIVERSITÄT MÜNCHEN Dr. W. Spann, G. Kleen WS 2022/2023

Programmieren II für Studierende der Mathematik

Aufgabe 11

Zur numerischen Lösung eines Systems gewöhnlicher Differentialgleichungen 2. Ordnung der Gestalt

\begin{matrix}y′′(t) = & f(t, y(t)) (t ∈ [t_{0}, T]), \\ y(t_{0}) = & y_{0}, \\ y′(t_{0}) = & y′_{0} \end{matrix}

mit vorgegebenen Anfangswerten y_{0}, y_{0}′ ∈ ℝ^{n} und rechter Seite f : [t_{0}, T] × ℝ^{n} → ℝ^{n} und gesuchter Funktion y : [t_{0}, T] → ℝ^{n} kann ein Nyström–Verfahren 4. Ordnung eingesetzt werden. Ein Verfahrensschritt mit Schrittweite h_{i} ist gegeben durch

k_{i}^{(1)} = f(t_{i}, y_{i})
k_{i}^{(2)} = f(t_{i} + \frac{h_{i}}{2} , y_{i} + \frac{h_{i}}{2} y_{i}′ + \frac{h_{i}^{2}}{8} k_{i}^{(1)})
k_{i}^{(3)} = f(t_{i} + h_{i}, y_{i} + h_{i}y_{i}′ + \frac{h_{i}^{2}}{2} k_{i}^{(2)})
y_{i+1} = y_{i} + h_{i}y′_{i} + \frac{h_{i}^{2}}{6} (k_{i}^{(1)} + 2k_{i}^{(2)})
y′_{i+1} = y′_{i} + \frac{h_{i}}{6} (k_{i}^{(1)} + 4k_{i}^{(2)} + k_{i}^{(3)}).

Programmieren Sie dieses Verfahren (ohne Schrittweitensteuerung) und testen Sie es an

(a) \begin{aligned}y′′_{1} & = & −5y_{1} + 4y_{2}, & y_{1}(0) = 0, y′_{1}(0) = 2, T = 6π,  h = \frac{π}{500} \\ y′′_{2} & = & 4y_{1} − 5y_{2}, & y_{2}(0) = 0, y′_{2}(0) = 0 \end{aligned}

Es sollen für die ersten und letzten 10 Schritte y_{i} und y′_{i} ausgedruckt werden.

(b) ¨ r⃗_{i} = −k^{2}(1 + \frac{m_{i}}{m_{0}} )\frac{r⃗_{i}}{|r⃗_{i}|^{3}} + k^{2} ∑_{{j=1 \atop j≠i}}^{9} \frac{m_{j}}{m_{0}} (\frac{r⃗_{j} − r⃗_{i}}{|r⃗_{j} − r⃗_{i}|^{3}} − \frac{r⃗_{j}}{|r⃗_{j}|^{3}} )  (i = 1, …, 9)

Dieses Differentialgleichungssystem beschreibt die Bewegung der 9 Planeten des Sonnensystems in heliozentrischen kartesischen Koordinaten, wobei lediglich die Anziehungskraft der Sonne und die gegenseitige Anziehung berücksichtigt ist. m_{0} ist hierin die Sonnenmasse, m_{1}, …, m_{9} die Massen der Planeten (Erde: Masse von Erde + Masse des Mondes), r⃗_{1}, …, r⃗_{9} die Ortsvektoren der Planeten–Schwerpunkte (Erde: Schwerpunkt des Erde/Mondsystems) und k die Gaußsche Gravitationskonstante (k^{2} = γm_{0}). Die Anfangswerte zum Zeitpunkt t_{0} (1.1.1997, 0^{h} Weltzeit) in heliozentrischen kartesischen Äquatorialkoordinaten zur Epoche 2000 und die Massenverhältnisse \frac{m_{0}}{m_{i}} zum Zeitpunkt t_{0} (1.1.1997, 0^{h} Weltzeit) finden Sie auf der Vorlesungshomepage. Die Länge wird dabei in AE (= 1.4959787066·10^{11}m), die Masse in Sonnenmassen und die Zeit in Tagen (= 86400s) gemessen. Die Gaußsche Konstante k beträgt 0.01720209895. Geben Sie die mit Schrittweite 1 Tag berechneten Positionen und Geschwindigkeiten der Planeten 200 Tage und 520 Tage nach bzw. 360 Tage vor dem Zeitpunkt t_{0} aus.

Versuchen Sie bei Programmierung, das Lösungsverfahren mit der Ausgabe, den Nyströmschritt und die Daten voneinander getrennt zu halten.

Bearbeitungszeitraum: bis Donnerstag, 2.2.2023, 16^{00}