Mathesis Wiki

Bis jetzt sind wir: Anton und Philipp.

Idee: Position von Körpern des Sonnensystems vorhersagen.

Ziel: a) die nächste Mondfinsternis berechnen

alternativ, falls das zu schwer ist b) berechnen, wann andere Planeten von der Erde aus sichtbar sind c) berechnen, wann man von wo die ISS sieht

spontane Inspirationen: https://de.wikipedia.org/wiki/Mondfinsternis#Schwierigkeiten_der_Vorausberechnung https://de.wikipedia.org/wiki/Mondfinsternis#Schwierigkeiten_der_Vorausberechnung http://www.diesonnenseite.de/startseite2/finsternisse_allg/berechnungen_zur_mondfinsternis/pdf/mond.pdf

Protokoll

Grundlagen des Projekts

In diesem Labor soll versucht werden eine Mondfinsternis so akkurat wie möglich vorherzusagen. Um dies zu Erreichen soll das System Erde-Sonne-Mond mittels einer Python Simulation berechnet werden um so die relative Positionen dieser Himmelskörper zu jedem Zeitpunkt zu ermitteln. Dabei müssen die zugrunde liegenden Differentialgleichungen gelöst werden und daraufhin geprüft werden ob die neu errechneten Positionen die Kriterien für eine Mondfinsternis erfüllen.

Differentialgleichungen Differentialgleichung sind Gleichungen, die nicht einen Wert als Unbekannte haben sondern eine Funktion. So ist in unserem Fall die Veränderung der Geschwindigkeit über die Zeit von der Beschleunigung abhängig. Von der Geschwindigkeit ist nun wieder die Veränderung der Position abhängig. Also müssen wir in jeden Zeitschritt zwei Differentialgleichungen lösen. Da Differentialgleichung in den meisten und auch in unserem Fall nicht analytisch exakt lösbar sind muss hier ein numerisches Verfahren angewendet werden.

Das heißt, dass ein Algorithmus benötigt wird der in jedem Zeitschritt immer den direkt folgenden Wert aus dem Ursprungswert berechnet. Dies kann man durch verschiedene Verfahren bewerkstelligen, wir haben in der eigentlichen Simulation das Runge-Kutta Verfahren 4.Ordnung genutzt. Zur Übung haben wir aber auch des einfachere, aber ungenauere Euler Verfahren implementiert

Im Euler Verfahren wird grundlegend die Ausgangsfunktion genommen und mit dem Zeitschritt multipliziert und das Ergebnis auf den Ausgangswert aufaddiert.

Das Runge-Kutta Verfahren 4.Ordung hingegen splittet die Gleichung zunächst in 4 Differenzenquotienten aufgeteilt und diese anschließen zusammenaddiert. Sie erfordert einen höheren Rechenaufwand bietet bei gleichem Zeitschritt aber ein deutlich genaueres Ergebnis.

Hier wird die Differentialgleichung x' = x mit den beiden Verfahren im Intervall 0-3,5 und dt = 0,175 mit den beiden Verfahren gelöst. Dies hat als analytische Lösung die bekannt e-Funktion und so lässt sich leicht die Genauigkeit überprüfen.

import numpy as np
import matplotlib.pyplot as plt
 
 
def euler (f, startwert, intervallende, anzahlschritte, x): # Euler Verfahren
	h = (intervallende - startwert)/anzahlschritte # berechnet den Zeitschritt abhängig von der den Anzahl der Schritte und dem Intervall
	tpoints2 = np.linspace(startwert, intervallende, anzahlschritte) # legt die Werte der x-Achse fest. teilt das Intervall in gleich große Schritte auf
	xpoints2 = [] # Liste für die Funktionswerte
	for t in tpoints2:
		xpoints2.append(x) # erweitert xliste in jedem Schritt um den jeweils nächsten x-Wert
		x += h * f(t,x) # berechnet den neuen Funktionswert nach Euler 
	return tpoints2, xpoints2
 
 
def rk4(f, startwert, intervallende, anzahlschritte, x): # Runge-Kutta-Verfahren 4.Ordnung
	h = (intervallende - startwert)/anzahlschritte # genau wie bei Euler
	tpoints = np.linspace(startwert, intervallende, anzahlschritte)
	xpoints = []
	for t in tpoints:
		xpoints.append(x)
		k1 = h*f(t,x) # berechnet die Teilfunktionen nach Runge-Kutta
		k2 = h*f(t+0.5*k1,x+0.5*h)
		k3 = h*f(t+0.5*k2,x+0.5*h)
		k4 = h*f(t+k3,x+h)
		x += (k1+2*k2+2*k3+k4)/6 # fasst die 4 Teilfunktionen zusammen und berechnet neuen Funktionswert
	return tpoints, xpoints
 
def f(t, x): # Funktion x' = x
    return x
 
vt, vx = rk4(f, 0., 3.5, 20, 1.) # legt Intervall, Anzahl der Schritte und Ausgangswert fest
vx1, vt1 = euler (f, 0., 3.5, 20, 1.)
 
 
 
plt.plot(vx1, vt1)  
plt.show() # zeigt das geplottete Ergebnis für Euler
plt.plot(vt, vx)
plt.show() # zeigt das geplottete Ergebnis für Runge-Kutta 4.Ordnung

Euler:

Hier erreicht die Funktion bei x = 1 den Funktionswert 2,4 was noch deutlich vom korrekten Wert e = 2,71… abweicht.

Runge-Kutta 4.Ordnung:

Hier erreicht die Funktion bei x = 1 einen Funktionswert von 2,52 was deutlich näher an dem gesuchten Wert liegt.

Vpython

Um eine graphische Ausgabe zu erhalten haben wir mit der Python Erweiterung VPython gearbeitet. In dieser Erweiterung ist eine Vektorklasse implementiert die in einem 3D-Modul Objekten Positionen, Form, Farbe etc übergibt und dieses dann darstellt. Als Übung haben wir zunächst ein kleines Programm mit sich statisch bewegenden Kugeln geschrieben.

from visual import *
import numpy as np
 
floor = box (pos=(0,0,0), length=50, height=0.5, width=4, color=color.blue)
ball = sphere (pos=(0,4,0), radius=1, color=color.red) # Ball wird erzeugt
ball.velocity = vector(5,-1,0) # Geschwindigkeit des Balls wird festgelegt
dt = 0.01
wall_left= box (pos= (25,5,0), length=1, height=10, widht=4, color=color.blue) # Begrenzungen werden festgelegt
wall_right= box (pos= (-25,5,0), length=1, height=10, widht=4, color=color.blue)
 
schnee=[] # List mit diversen kleinen "Schneekugeln"
for i in range (5):
	schnee.append(sphere(pos= ((-25-25)*np.random.rand()+25,10,0) , radius=0.1 , color=color.white)) # Position wird partiell zufaellig erzeugt
 
for flocke in schnee:
	flocke.velocity= vector (0.,-1.,0.) # Geschwindigkeit wird für jede Schneeflocke festgelegt
	flocke.timer=1
 
while 1: # Simulation läuft endlos
	rate (100) # 100 Schleifendurchläufe werden pro Sekunde durchgeführt
	ball.pos = ball.pos + ball.velocity*dt # Position des Balls wird nach Euler geupdated
 
	if ball.y < ball.radius:
		ball.velocity.y = abs(ball.velocity.y)
	else:
		ball.velocity.y = ball.velocity.y - 9.8*dt # Verhalten des Balls im an den Rändern der erzeugten Box wird festgelegt
 
	if ball.x > 25-ball.radius:
		ball.velocity.x = -(ball.velocity.x)
	if ball.x < -25+ball.radius:
		ball.velocity.x = -(ball.velocity.x)
 
	if len(schnee)>0 and schnee[-1].y < 9: # Berechnung der Schneeflocken
		for i in range (5):
			schnee.append(sphere(pos= ((-25-25)*np.random.rand()+25,10,0) , radius=0.1 , color=color.white))
			schnee[-1].velocity=vector (0,-1,0)
			schnee[-1].timer=1
	zu_loeschend=[]
	for i,flocke in enumerate(schnee):
		#flocke.velocity= vector (0,-5,0)
		flocke.pos= flocke.pos+flocke.velocity*dt
		flocke.timer=flocke.timer+1
		if flocke.y < flocke.radius+0.18: # legt fest wenn die Scheeflocke den "Boden" erreicht hat
			flocke.velocity=vector(0.,0.,0.)
		else:
		    flocke.velocity= vector (0.,-1.,0.)
		if flocke.timer > 500: # loescht Schneeflocken wenn ihre Anzahl 500 übersteigt
			zu_loeschend.append(i)
 
	schnee_neu=[ f for (i,f) in enumerate(schnee) if i not in zu_loeschend]

Simulation des Sonnensystems

Um das S