Hier werden die Unterschiede zwischen zwei Versionen gezeigt.
Beide Seiten der vorigen Revision Vorhergehende Überarbeitung Nächste Überarbeitung | Vorhergehende Überarbeitung | ||
ss17:viele_dinge_fliegen_im_weltall_durcheinander:testprogramme [2017/08/23 12:01] mtoppermann [Vergleich] |
ss17:viele_dinge_fliegen_im_weltall_durcheinander:testprogramme [2017/08/28 15:41] (aktuell) mtoppermann [Vergleich von Planetenbewegungen mit verschiedenen Methoden] |
||
---|---|---|---|
Zeile 1: | Zeile 1: | ||
====== Testprogramme ====== | ====== Testprogramme ====== | ||
- | ===== Vergleich ===== | + | Diese Testprogramme helfen uns dabei gewisse Funktionen auf ihre Richtigkeit zu prüfen, wie auch zum Vergleichen welche Funktion bessere Ergebnisse liefert. |
+ | |||
+ | ===== Vergleich von Planetenbewegungen mit verschiedenen Methoden===== | ||
+ | |||
+ | In diesem Skript wird die Genauigkeit des Euler-Verfahrens mit dem Leapfrog-Verfahren am Bespiel der Erdbewegung | ||
+ | um die Sonne verglichen. Die Bewegung der dunkelblauen Erde wir durch das Leapfrog-Verfahren berechnet, das der hellblauen durch das Euler-Verfahren. Nach einiger Zeit sieht man bereits deutliche Unterschiede in der Bahnbewegung: Die dunkelblaue Erde verlässt die richtige Bahn langsamer, da ihre die Berechnungsfehler kleiner sind. | ||
<code python> | <code python> | ||
from functions import * | from functions import * | ||
import pygame | import pygame | ||
+ | |||
class Sun(object): | class Sun(object): | ||
- | """Our Sun in the point of origin | + | |
- | """ | + | |
def __init__(self): | def __init__(self): | ||
+ | """initialize all necessary properties of the sun | ||
+ | """ | ||
self.x = 0 | self.x = 0 | ||
self.y = 0 | self.y = 0 | ||
self.radius = 20 | self.radius = 20 | ||
self.mass = 1.989e30 | self.mass = 1.989e30 | ||
- | self.colour = (255,255,0) | + | self.colour = (255, 255, 0) |
- | + | ||
- | class Planet(object): | + | class Earth(object): |
- | """An ordinary planet, in this case our Earth | + | |
- | """ | + | def __init__(self, color): |
- | def __init__(self,color): | + | """initialize all necessary properties of the earth |
+ | """ | ||
self.x = 1.49e11 | self.x = 1.49e11 | ||
self.y = 0. | self.y = 0. | ||
self.radius = 10 | self.radius = 10 | ||
self.mass = 5.97e24 | self.mass = 5.97e24 | ||
- | self.velocity = np.array([0.,2.978e4]) | + | self.velocity = np.array([0., 2.978e4]) |
- | self.force = np.array([0,0]) | + | self.force = np.array([0, 0]) |
self.colour = (color) | self.colour = (color) | ||
- | self.positions = [(int(self.x * 1e-9 + width/2), | + | |
- | int(self.y * 1e-9 + height/2))] | + | |
- | + | ||
- | # leapfrog method for DE-solving | + | |
def move(self): | def move(self): | ||
+ | """leapfrog method for differential-equation-solving | ||
+ | """ | ||
self.x += 0.5 * self.velocity[0] * dt | self.x += 0.5 * self.velocity[0] * dt | ||
self.y += 0.5 * self.velocity[1] * dt | self.y += 0.5 * self.velocity[1] * dt | ||
Zeile 38: | Zeile 44: | ||
self.x += 0.5 * self.velocity[0] * dt | self.x += 0.5 * self.velocity[0] * dt | ||
self.y += 0.5 * self.velocity[1] * dt | self.y += 0.5 * self.velocity[1] * dt | ||
- | + | ||
- | # euler method for DE-solving | + | |
def move2(self): | def move2(self): | ||
+ | """euler-method for differential-equation-solving | ||
+ | """ | ||
self.x += self.velocity[0] * dt | self.x += self.velocity[0] * dt | ||
self.y += self.velocity[1] * dt | self.y += self.velocity[1] * dt | ||
self.velocity += self.force/self.mass * dt | self.velocity += self.force/self.mass * dt | ||
+ | |||
def draw(a, points): | def draw(a, points): | ||
- | """Draws an object and its trail on the screen | + | """Draws an object on the screen |
""" | """ | ||
pygame.draw.circle(screen, a.colour, (int(a.x * 10e-10) + width/2, | pygame.draw.circle(screen, a.colour, (int(a.x * 10e-10) + width/2, | ||
int(a.y * 10e-10) + height/2), a.radius, 0) | int(a.y * 10e-10) + height/2), a.radius, 0) | ||
pygame.draw.aalines(screen, a.colour, False, points, 2) | pygame.draw.aalines(screen, a.colour, False, points, 2) | ||
+ | |||
def draw2(a): | def draw2(a): | ||
"""Draws an object in the middle of the screen | """Draws an object in the middle of the screen | ||
""" | """ | ||
pygame.draw.circle(screen, a.colour, (width/2, height/2), a.radius, 0) | pygame.draw.circle(screen, a.colour, (width/2, height/2), a.radius, 0) | ||
- | + | ||
- | def last_elems(elems, k): | + | |
- | """Returns the last k elements of a list | + | |
- | """ | + | |
- | tmp = [] | + | |
- | for i in xrange(k): | + | |
- | if i > len(elems): | + | |
- | break | + | |
- | elif i < len(elems): | + | |
- | tmp.append(elems[-(i+1)]) | + | |
- | return tmp | + | |
dt = 36000 | dt = 36000 | ||
- | (width, height) = (1000,600) | + | (width, height) = (1000, 600) |
screen = pygame.display.set_mode((width, height)) | screen = pygame.display.set_mode((width, height)) | ||
pygame.display.set_caption(" ") | pygame.display.set_caption(" ") | ||
- | background_colour = (10,10,10) | + | bgcolour = (10, 10, 10) |
- | + | ||
- | erde = Planet((0,0,255)) | + | earth = Earth((0, 0, 255)) |
- | kerbin = Planet((50,150,255)) | + | kerbin = Earth((50, 150, 255)) |
- | sonne = Sun() | + | sun = Sun() |
- | + | ||
- | count = 0 | + | |
running = True | running = True | ||
while running: | while running: | ||
Zeile 84: | Zeile 79: | ||
if event.type == pygame.QUIT: | if event.type == pygame.QUIT: | ||
running = False | running = False | ||
- | + | ||
- | #every ten steps, the trail is updated | + | |
- | if count % 10 == 0: | + | |
- | erde.positions.append((int(1e-9 * erde.x) + width/2, | + | |
- | int(1e-9 * erde.y) + height/2)) | + | |
- | kerbin.positions.append((int(1e-9 * kerbin.x) + width/2, | + | |
- | int(1e-9 * kerbin.y) + height/2)) | + | |
- | linepoints1 = last_elems(erde.positions, 200) | + | |
- | linepoints2 = last_elems(kerbin.positions, 200) | + | |
# display everything | # display everything | ||
- | screen.fill(background_colour) | + | screen.fill(bgcolour) |
draw(erde, linepoints1) | draw(erde, linepoints1) | ||
draw(kerbin, linepoints2) | draw(kerbin, linepoints2) | ||
draw2(sonne) | draw2(sonne) | ||
pygame.display.update() | pygame.display.update() | ||
+ | |||
# calculate everything | # calculate everything | ||
- | erde.force = get_force(erde, sonne) | + | erde.force = get_force(earth, sun) |
- | kerbin.force = get_force(kerbin, sonne) | + | kerbin.force = get_force(kerbin, sun) |
erde.move() | erde.move() | ||
kerbin.move2() | kerbin.move2() | ||
- | + | ||
- | count += 1 | + | |
</code> | </code> | ||
+ | =====Euler-Verfahren vs Verlet(Leapfrog)-Verfahren===== | ||
+ | |||
+ | <code python> | ||
+ | import matplotlib.pyplot as plt | ||
+ | import numpy as np | ||
+ | |||
+ | # x .. = f(x) | ||
+ | def f(x): | ||
+ | return -x | ||
+ | |||
+ | def leapfrog(startwertx, startwertv, intervallende, anzahlschritte): | ||
+ | h = (intervallende - startwertx)/anzahlschritte | ||
+ | tpoints = np.linspace(0, intervallende, anzahlschritte) | ||
+ | xpoints = [] | ||
+ | vpoints = [] | ||
+ | x = startwertx | ||
+ | v = startwertv | ||
+ | for t in tpoints: | ||
+ | xpoints.append(x) | ||
+ | vpoints.append(v) | ||
+ | x += 0.5 * v * h | ||
+ | v += f(x) * h | ||
+ | x += 0.5 * v * h | ||
+ | #x += h * f(x) + f(x) * h**2/2 | ||
+ | return tpoints, xpoints | ||
+ | |||
+ | def euler(startwertx, startwertv, intervallende, anzahlschritte): | ||
+ | h = (intervallende - startwertx)/anzahlschritte | ||
+ | tpoints = np.linspace(0, intervallende, anzahlschritte) | ||
+ | xpoints = [] | ||
+ | vpoints = [] | ||
+ | x = startwertx | ||
+ | v = startwertv | ||
+ | for t in tpoints: | ||
+ | xpoints.append(x) | ||
+ | vpoints.append(v) | ||
+ | v_neu = v + h * f(x) | ||
+ | x += h * v | ||
+ | v = v_neu | ||
+ | #x += h * f(x) + f(x) * h**2/2 | ||
+ | return tpoints, xpoints | ||
+ | |||
+ | |||
+ | t1, x1 = leapfrog(0.,1., 100, 500) | ||
+ | t2, x2 = euler(0.,1., 100, 500) | ||
+ | pointsin = np.sin(t1) | ||
+ | abweichung_1 = [] | ||
+ | abweichung_2 = [] | ||
+ | for i in range(len(t1)): | ||
+ | abweichung_1.append(pointsin[i]-x1[i]) | ||
+ | abweichung_2.append(pointsin[i]-x2[i]) | ||
+ | |||
+ | fig = plt.figure(figsize = (18, 8.5)) | ||
+ | ax = fig.add_subplot(111) | ||
+ | |||
+ | ax.plot(t1, x1, label="Leapfrog") | ||
+ | ax.plot(t2, x2, label="Euler") | ||
+ | ax.plot(t1, np.sin(t1), label="np.Sinus") | ||
+ | ax.plot(t1, abweichung_1, label="Abweichung Leapfrog") | ||
+ | ax.plot(t1, abweichung_2, label="Abweichung Euler") | ||
+ | ax.plot(t1, 0*t1, color="black") | ||
+ | plt.legend(loc=best) | ||
+ | plt.show() | ||
+ | </code> | ||
+ | |||
+ | =====Asteoridengürtel nur mit Sonnengravitation===== | ||
+ | |||
+ | <code python> | ||
+ | import numpy as np | ||
+ | import pygame | ||
+ | from functions import get_force | ||
+ | |||
+ | class Asteorid(object): | ||
+ | def __init__(self,phi,radius): | ||
+ | self.x = np.cos(phi) * radius | ||
+ | self.y = np.sin(phi) * radius | ||
+ | self.radius = np.random.randint(1,3) | ||
+ | self.mass = 6.687e15 | ||
+ | self.velocity = np.sqrt(G * 1.989e30/radius) | ||
+ | * np.array([-np.sin(phi), np.cos(phi)]) | ||
+ | self.force = np.array([0.,0.]) | ||
+ | self.colour = (100,100,100) | ||
+ | |||
+ | def move(self): | ||
+ | self.x += 0.5 * self.velocity[0] * dt | ||
+ | self.y += 0.5 * self.velocity[1] * dt | ||
+ | self.velocity += self.force/self.mass * dt | ||
+ | self.x += 0.5 * self.velocity[0] * dt | ||
+ | self.y += 0.5 * self.velocity[1] * dt | ||
+ | |||
+ | class Sun(object): | ||
+ | """Our Sun in the point of origin | ||
+ | """ | ||
+ | def __init__(self): | ||
+ | self.x = 0 | ||
+ | self.y = 0 | ||
+ | self.radius = 20 | ||
+ | self.mass = 1.989e30 | ||
+ | self.colour = (255,255,0) | ||
+ | |||
+ | def draw(a): | ||
+ | """Draws an object and its trail on the screen | ||
+ | """ | ||
+ | pygame.draw.circle(screen, a.colour, (int(a.x * 0.5e-9) + width/2, | ||
+ | int(a.y * 0.5e-9) + height/2), a.radius, 0) | ||
+ | |||
+ | def draw2(a): | ||
+ | """Draws an object in the middle of the screen | ||
+ | """ | ||
+ | pygame.draw.circle(screen, a.colour, (width/2, height/2), a.radius, 0) | ||
+ | |||
+ | G = 6.673e-11 | ||
+ | AE = 1.496e11 | ||
+ | n = 2000 | ||
+ | dt = 36000 | ||
+ | liste = [] | ||
+ | |||
+ | (width, height) = (1000,600) | ||
+ | screen = pygame.display.set_mode((width, height)) | ||
+ | pygame.display.set_caption(" ") | ||
+ | background_colour = (10,10,10) | ||
+ | |||
+ | for i in xrange(n): | ||
+ | phi = 2 * np.random.random() * np.pi | ||
+ | radius = np.random.uniform(2.,3.4) * AE | ||
+ | obj = Asteorid(phi, radius) | ||
+ | liste.append(obj) | ||
+ | |||
+ | sonne = Sun() | ||
+ | |||
+ | running = True | ||
+ | while running: | ||
+ | for event in pygame.event.get(): | ||
+ | if event.type == pygame.QUIT: | ||
+ | running = False | ||
+ | |||
+ | # display everything | ||
+ | screen.fill(background_colour) | ||
+ | draw2(sonne) | ||
+ | for i in liste: | ||
+ | draw(i) | ||
+ | i.force = get_force(i, sonne) | ||
+ | i.move() | ||
+ | pygame.display.update() | ||
+ | </code> | ||
+ | |||
+ | =====Neue Berechnung des Drehimpulses nach einem Zusammenstoß===== | ||
+ | |||
+ | <code python> | ||
+ | # -*- coding: utf8 -*- | ||
+ | from functions import * | ||
+ | import numpy as np | ||
+ | n = 2 | ||
+ | liste_teilchen = init_particle_list(n) | ||
+ | part_1,part_2=liste_teilchen | ||
+ | part_1.velocity = np.array([1.,0.]) | ||
+ | part_2.velocity = np.array([0.,1.]) | ||
+ | |||
+ | |||
+ | def get_torque(part_1, part_2): | ||
+ | """ | ||
+ | Drehimpulsberechnung für neu entstandenes Teilchen nach Kollision | ||
+ | """ | ||
+ | newpoint_x = (part_1.mass * part_1.x + part_2.mass * part_2.x) | ||
+ | / (part_1.mass + part_2.mass) | ||
+ | newpoint_y = (part_1.mass * part_1.y + part_2.mass * part_2.y) | ||
+ | / (part_1.mass + part_2.mass) | ||
+ | newpoint = ([newpoint_x, newpoint_y]) #Schwerpunkt zwischen beiden Teilchen | ||
+ | print part_1.x,part_1.y | ||
+ | print part_2.x,part_2.y | ||
+ | print newpoint | ||
+ | |||
+ | r_1 = ([part_1.x - newpoint_x, part_1.y - newpoint_y]) #Richtungsvektoren | ||
+ | r_2 = ([part_2.x - newpoint_x, part_2.y - newpoint_y]) | ||
+ | |||
+ | print r_1 | ||
+ | print r_2 | ||
+ | |||
+ | torq_1 = part_1.mass * np.cross(r_1, part_1.velocity) #Drehimpulse der einzelnen Teilchen zum Schwerpunkt | ||
+ | torq_2 = part_2.mass * np.cross(r_2, part_2.velocity) | ||
+ | torq = torq_1 + torq_2 #Gesammtdrehimpuls (+ Eigendrehimpulse falls vorhanden) | ||
+ | |||
+ | print torq | ||
+ | return newpoint,torq | ||
+ | |||
+ | for i,part_1 in enumerate(liste_teilchen): | ||
+ | for j,part_2 in enumerate(liste_teilchen): | ||
+ | if i < j: | ||
+ | get_torque(part_1, part_2) | ||
+ | </code> | ||
+ |