Diese Testprogramme helfen uns dabei gewisse Funktionen auf ihre Richtigkeit zu prüfen, wie auch zum Vergleichen welche Funktion bessere Ergebnisse liefert.

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.

from functions import *
import pygame
 
class Sun(object):
 
    def __init__(self):
        self.x = 0
        self.y = 0
        self.radius = 20
        self.mass = 1.989e30
        self.colour = (255, 255, 0)
 
class Earth(object):
 
    def __init__(self, color):
        self.x = 1.49e11
        self.y = 0.
        self.radius = 10
        self.mass = 5.97e24
        self.velocity = np.array([0., 2.978e4])
        self.force = np.array([0, 0])
        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):
        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
 
    # euler method for DE-solving
    def move2(self):
        self.x += self.velocity[0] * dt
        self.y += self.velocity[1] * dt
        self.velocity += self.force/self.mass * dt
 
def draw(a, points):
    """Draws an object on the screen
    """
    pygame.draw.circle(screen, a.colour, (int(a.x * 10e-10) + width/2,
    int(a.y * 10e-10) + height/2), a.radius, 0)
    pygame.draw.aalines(screen, a.colour, False, points, 2)
 
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)
 
dt = 36000
(width, height) = (1000,600)
screen = pygame.display.set_mode((width, height))
pygame.display.set_caption(" ")
background_colour = (10,10,10)
 
earth = Earth((0, 0, 255))
kerbin = Earth((50, 150, 255))
sonne = Sun()
 
count = 0
running = True
while running:
    for event in pygame.event.get():
        if event.type == pygame.QUIT:
            running = False
 
    # display everything
    screen.fill(background_colour)
    draw(erde, linepoints1)
    draw(kerbin, linepoints2)
    draw2(sonne)
    pygame.display.update()
 
    # calculate everything
    erde.force = get_force(erde, sonne)
    kerbin.force = get_force(kerbin, sonne)
    erde.move()
    kerbin.move2()
 
    count += 1   
 

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()

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()

# -*- 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)