Simulation anziehender und abstoßender Felder in Python 3 (Software)

Martin Vogel , Dortmund / Bochum, Sun, 01.01.2017, 17:54 (vor 2840 Tagen)

Eine Simulation ohne allzu konkreten physikalischen Bezug zeigt die Wirkung zweier entgegengesetzt wirkender Felder auf eine Anzahl Partikel. Ähnlich der Gravitation werden die Partikel mit einer Kraft angezogen, die sich umgekehrt quadratisch zum Abstand der einzelnen Teilchen verhält. Gleichzeitig werden sie von einer Kraft mit stärkerer Wirkung aber geringerer effektiver Reichweite abgestoßen. Ganz entfernt erinnert das Verhalten an die Wechselwirkung der Elementarkräfte im Atomkern.

Sinn der Übung ist weniger die Darstellung eines realen physikalischen Vorgangs (insofern muss ich alle enttäuschen, die das als Last-Minute-Idee für die aktuelle Programmieraufgabe aufgreifen wollen) als vielmehr, ein Programmierbeispiel für den Einfluss eines Reglerpaneels auf eine laufende Simulation zu sein. Diese wird zu Beginn des Programms einmal angestoßen und läuft dann kontinuierlich weiter, während alle möglichen Parameter bis hin zur Partikelanzahl lustig geändert werden können. Gelegentlich läuft die Simulation total aus dem Ruder, was sich in der Regel durch entschlossenes Hochziehen des Dämpfungsreglers schnell beheben lässt.

Wer kann mir das kuriose Verhalten der Partikel erklären, wenn man den Regler für die Abstoßung ganz oder fast ganz auf null herunterzieht?

Download: ZIP-Archiv_2ABK27ET1.zip

#!/usr/bin/python3
 
# Ein Haufen Partikel, der sich auf große Entfernung anzieht
# und auf geringe Entfernung abstößt
 
#
# Version vom 1. Januar 2017
# Autor: Martin Vogel, martinvogel.de
# Lizenz: cc-by-3.0 https://creativecommons.org/licenses/by/3.0/de/
#
 
from tkinter import *
from random import randint, choice, random
from time import perf_counter
from math import hypot
 
 
class Einstellungen:
    """Einige programmweit gültige Einstellungen, die meisten sind über
    # Schieberegler beeinflussbar."""
 
    def __init__(self):
        # Die Anfangsgröße der Canvas
        self.Höhe = 600
        self.Breite = 800
 
        # Die Anfangszahl der Partikel
        self.Anzahl = 200
 
        # Die Elastizität der Wand:
        # 1: Impuls wird voll zurückgegeben, 0: Partikel "klebt" an der Wand.
        self.Wandelastizität = 0.5
 
        # Anteil, um den die Geschwindigkeit pro Rechenschritt vermindert wird.
        self.Dämpfung = 0.01
        # Anmerkung: da die Zahl der Rechenschritte pro Sekunde von der Anzahl
        # der Partikel abhängt (der Rechenaufwand wächst quadratisch), wirkt
        # sich die Dämpfung bei geringer Partikelzahl viel stärker aus als bei
        # großer Partikelzahl. Wer mag, kann das gerne mal verbessern.
 
        # Die Stärke von Anziehung und Abstoßung.
        self.Anziehungsfaktor = 9
        self.Abstoßungsfaktor = 18
 
        # Bezugswert für Zeitmessungen:
        self.Zeit = perf_counter()
 
    def setze_dämpfung(self, event):
        # Lies den Schieberwert!
        self.Dämpfung = SD.get()
 
    def setze_anzahl(self, event=None):
        # Wo steht der Schieber?
        self.Anzahl = SA.get()
        while len(Partikel) < self.Anzahl:
            # Brauchen wir noch Partikel? Neuen erzeugen!
            C_Partikel()
        while len(Partikel) > self.Anzahl:
            # Haben wir zu viele Partikel? Irgendeins löschen!
            choice(Partikel).lösche()
 
    def setze_anziehungsfaktor(self, event):
        # Hole den Schieberwert!
        self.Anziehungsfaktor = SAn.get()
        # Die Grafik im Kontrollpaneel aktualisieren
        C.after_idle(ZeigeVerlauf)
 
    def setze_abstoßungsfaktor(self, event):
        # Welcher Schieberwert wurde eingestellt?
        self.Abstoßungsfaktor = SAb.get()
        # Die Grafik im Kontrollpaneel aktualisieren
        C.after_idle(ZeigeVerlauf)
 
    def setze_wandelastizität(self, event):
        # Lies den Schieber, setze das Attribut:
        self.Wandelastizität = SW.get()
 
    def setze_canvasgröße(self, event):
        # Falls die Größe der Canvas verändert wird, neue Abmessungen auslesen
        # und merken:
        self.Breite = C.winfo_width()
        self.Höhe = C.winfo_height()
 
 
# Die einzige Instanz dieser Klasse:
E = Einstellungen()
 
# Alle Partikel werden in dieser Liste verwaltet.
Partikel = []
 
 
def zufallsfarben():
    "Erzeugt einen nicht zu dunklen RGB-Farbwert und eine dunklere Variante"
    # Ein Farbpaar für die Partikel: Heller Kern, dunkler Rand
    r, g, b = randint(64, 255), randint(64, 255), randint(64, 255)
    return "#%02x%02x%02x" % (r, g, b), "#%02x%02x%02x" % (r//2, g//2, b//2)
 
 
class C_Partikel:
    "Massepunkt mit Ort und Geschwindigkeit."
    def __init__(self, x=None, y=None):
        # Ort
        if x is None:
            x = randint(0, E.Breite)
        if y is None:
            y = randint(0, E.Höhe)
        self.x = x
        self.y = y
        # Geschwindigkeit
        self.vx = 0
        self.vy = 0
        # Farbe
        hell, dunkel = zufallsfarben()
        # Canvasobjekt
        self.ID = C.create_oval(self.x-3, self.y-3,
                                self.x+3, self.y+3,
                                fill=hell,
                                outline=dunkel)
        # Eintragen in die Liste aller Partikel
        Partikel.append(self)
        # Der Punkt kann mit der Maus verschoben werden.
        C.tag_bind(self.ID, "<B1-Motion>", self.bewegezu)
 
    def lösche(self):
        "Punkt wird entfernt"
        # Lösche Punkt aus der Partikelliste
        Partikel.remove(self)
        # Lösche Canvas-Objekt
        C.delete(self.ID)
 
    def bewegezu(self, event):
        "Punkt wird verschoben"
        C.move(self.ID, event.x-self.x, event.y-self.y)
        self.x, self.y = event.x, event.y
 
 
def Animationsschleife(event=None):
    """Einsammeln aller Kräfte, die auf jedes Partikel wirken und anschließende
    Berechnung der Kraftwirkung"""
 
    # Berechnung der seit dem letzten Aufruf verstrichenen Zeit delta_t
    jetzt = perf_counter()
    delta_t = jetzt-E.Zeit
    E.Zeit = jetzt
 
    Anz, Abs = E.Anziehungsfaktor, E.Abstoßungsfaktor
    Dmp = E.Dämpfung
    Wez = E.Wandelastizität
 
    # Kopie der Partikelliste erzeugen, damit es keine Verfälschungen durch
    # Neuberechnungen gibt.
    P = Partikel[:]
    # Für alle Partikel:
    for i, Pi in enumerate(P):
        # Aktuelle Koordinaten
        x, y = Pi.x, Pi.y
        # Beschleunigungssumme in x- und y-Richtung
        ax = ay = 0
        # Einwirkungen aller anderen Partikel einsammeln:
        for Pj in P:
            # Abstände in x- und y-Richung
            dx = Pj.x-Pi.x
            dy = Pj.y-Pi.y
            # Abstand diagonal
            d = hypot(dx, dy)
            if d:
                # Falls Pj nicht dasselbe Partikel wie Pi ist:
                # Die Beschleunigungsgleichung:
                # Die Anziehung lässt mit dem Quadrat der Entfernung nach
                # (ähnlich der Gravitation).
                # Die Abstoßung soll hier mit der dritten Potenz der
                # Entfernung schwächer werden (ohne physikalische Grundlage).
                a = (Anz/d)**2 - (Abs/d)**3
                # Beschleunigungskomponenten aufaddieren
                ax += dx*a
                ay += dy*a
        # Nachdem die resultierende Beschleunigung aufaddiert wurde, kann
        # die Geschwindigkeitsänderung berechnet werden:
        Partikel[i].vx = (ax*delta_t+Pi.vx)*(1-Dmp)
        Partikel[i].vy = (ay*delta_t+Pi.vy)*(1-Dmp)
        # Daraus ergibt sich der neue Ort des Partikels.
        x += Partikel[i].vx*delta_t
        y += Partikel[i].vy*delta_t
        # Bei den folgenden Kollisionen wird das Partikel um eine zufällige
        # Entfernung von der Wand abgesetzt, um kuriose Effekte zu vermeiden,
        # die entstehen, wenn ein Haufen benachbarter Partikel dieselben
        # x- oder y- Werte aufweisen.
 
        if x < 0:
            # Kollision mit linker Wand?
            x = random()/10
            Partikel[i].vx *= -Wez
            Partikel[i].vy *= Wez
 
        elif x > E.Breite:
            # Kollision mit rechter Wand?
            x = E.Breite-random()/10
            Partikel[i].vx *= -Wez
            Partikel[i].vy *= Wez
 
        if y < 0:
            # Kollision mit oberer Wand?
            y = random()/10
            Partikel[i].vx *= Wez
            Partikel[i].vy *= -Wez
 
        elif y > E.Höhe:
            # Kollision mit unterer Wand?
            y = E.Höhe-random()/10
            Partikel[i].vx *= Wez
            Partikel[i].vy *= -Wez
 
        # Canvasobjekt des Partikels auf neue Koordinaten setzen
        C.move(Partikel[i].ID, x-Partikel[i].x, y-Partikel[i].y)
        Partikel[i].x, Partikel[i].y = x, y
 
    # Sobald die CPU wieder Luft geholt hat, geht’s von vorne los.
    C.after_idle(Animationsschleife)
 
# Das Hauptfenster mit der Canvas
T = Tk()
T.title("Wechselwirkung nah und fern")
 
# Die Zeichenfläche. Sie passt sich Größenänderungen des Fensters an.
C = Canvas(width=E.Breite, height=E.Höhe, bg="black", highlightthickness=0)
C.pack(expand=True, fill="both")
C.bind("<Configure>", E.setze_canvasgröße)
 
# Das Kontrollfenster
KF = Toplevel()
KF.title("Reglerpult")
# Das Reglerpult
SLF = LabelFrame(KF, text="Parameter", padx=5, pady=5)
SLF.grid(row=0, column=0, columnspan=2, sticky="wens", padx=5, pady=5)
 
# Zeile 1 im Labelframe (die Regler) soll dessen Größenanpassung übernehmen,
# außerdem soll sie mindestens 200 Pixel hoch sein, damit die
# Skalenzahlen gut zu erkennen sind.
KF.rowconfigure(0, weight=1)
SLF.rowconfigure(1, weight=1, minsize=200)
# Die 5 Spalten sollen sich gleichmäßig verteilen können.
SLF.columnconfigure(0, weight=1)
SLF.columnconfigure(1, weight=1)
SLF.columnconfigure(2, weight=1)
SLF.columnconfigure(3, weight=1)
SLF.columnconfigure(4, weight=1)
 
# Schieberegler für die Anzahl der Partikel
Label(SLF, text="Anzahl").grid(row=0, column=0)
SA = Scale(SLF, from_=500, to=0, width=20, orient="vertical",
           tickinterval=100, resolution=1, command=E.setze_anzahl)
SA.set(E.Anzahl)  # Anfangswert des Reglers
SA.grid(row=1, column=0, sticky="nsew")
 
# Schieberegler für den Anziehungsfaktor
Label(SLF, text="Anziehung").grid(row=0, column=1)
SAn = Scale(SLF, from_=50, to=0, width=20, orient="vertical",
            tickinterval=10, resolution=0.1, command=E.setze_anziehungsfaktor)
SAn.set(E.Anziehungsfaktor)  # Anfangswert des Reglers
SAn.grid(row=1, column=1, sticky="nsew")
 
# Schieberegler für den Abstoßungsfaktor
Label(SLF, text="Abstoßung").grid(row=0, column=2)
SAb = Scale(SLF, from_=50, to=0, width=20, orient="vertical",
            tickinterval=10, resolution=0.1, command=E.setze_abstoßungsfaktor)
SAb.set(E.Abstoßungsfaktor)  # Anfangswert des Reglers
SAb.grid(row=1, column=2, sticky="nsew")
 
# Schieberegler für die Wand-Elastizität
Label(SLF, text="Wand-nelastizität").grid(row=0, column=3)
SW = Scale(SLF, from_=1, to=0, width=20, orient="vertical",
           tickinterval=0.2, resolution=0.01, command=E.setze_wandelastizität)
SW.set(E.Wandelastizität)  # Anfangswert des Reglers
SW.grid(row=1, column=3, sticky="nsew")
 
# Schieberegler für die Dämpfung
Label(SLF, text="Dämpfung").grid(row=0, column=4)
SD = Scale(SLF, from_=1, to=0, width=20, orient="vertical",
           tickinterval=0.2, resolution=0.001, command=E.setze_dämpfung)
SD.set(E.Dämpfung)  # Anfangswert des Reglers
SD.grid(row=1, column=4, sticky="nsew")
 
 
def ZeigeVerlauf(event=None):
    """In einer Canvasgrafik unterhalb der Regler wird der Wirkungsradius
    und die Stärke von Anziehung und Abstoßung (innerhalb gewisser Grenzen)
    als Diagramm angezeigt."""
    # Einige Werte zur Darstellung ausrechnen und zwischenspeichern
    Werte = [0]
    for d in range(1, 600):
        Anziehung = (E.Anziehungsfaktor/d)**2
        Abstoßung = (E.Abstoßungsfaktor/d)**3
        Werte.append(Anziehung-Abstoßung)
 
    # Extrema?
    Wmin, Wmax = min(Werte)-0.001, max(Werte)+0.001
    # Der betragsmäßig kleinere der beiden Werte bestimmt den Maßstab
    Wmin, Wmax = max(Wmin, -2*Wmax), min(Wmax, -2*Wmin)
    # … noch etwas Randabstand
    Wmax += 0.05*(Wmax-Wmin)
    Wmin -= 0.05*(Wmax-Wmin)
 
    # Lage der Nulllinie in y-Richtung
    Null = 100*(Wmax)/(Wmax-Wmin)
 
    # Punkte des Funktionsgraphen
    Punkte = []
 
    # Skalierung der Werte auf den Darstellungsbereich
    for x in range(1, 600):
        Punkte.append(x)
        y = min(101, max(-1, 100*(Wmax-Werte[x])/(Wmax-Wmin)))
        Punkte.append(y)
 
    # Radius der Grenze zwischen Anziehung und Abstoßung
    try:
        r0 = E.Abstoßungsfaktor**3/E.Anziehungsfaktor**2
    except:
        r0 = float("inf")
 
    # Radius wird für die Darstellung auf Anzeigebereich begrenzt
    x0 = min(600, r0)
 
    # Grafik neu aufbauen
    SC.delete("all")
    # Nullinie
    SC.create_line(0, Null, 600, Null, fill="#008000")
    # An-Ab-Grenze
    SC.create_line(x0, 0, x0, 100, fill="#008000")
    # untenliegende Beschriftung
    SC.create_text(x0+2, 2, text="Anziehung", anchor="nw", fill="#40C000")
    SC.create_text(x0-2, 2, text="Abstoßung", anchor="ne", fill="#40C000")
    # Funktionsgraph
    SC.create_line(Punkte, fill="#E0E000")
    # Beschleunigungsgleichung
    SC.create_text(598, 98, text="a = (%.1f/x)² - (%.1f/x)³" % (
                   E.Anziehungsfaktor, E.Abstoßungsfaktor),
                   anchor="se", fill="#E0E000")
    # Radius als gestrichelter Kreis mit symbolischem Partikel im Zentrum
    SC.create_oval(-3, Null-3, 3, Null+3, fill="red", outline="red")
    SC.create_oval(-x0, Null-x0, x0, Null+x0, outline="red", dash=(3, ))
    SC.create_text(x0/2, Null-1, text="r0 = %.1f" % r0, fill="red", anchor="s")
 
# Die Canvas im Kontrollfenster
SC = Canvas(KF, width=600, height=100, bg="black", highlightthickness=0)
SC.grid(row=1, column=0, padx=5, pady=5)
 
# Anzeige des Diagramms
ZeigeVerlauf()
 
# Canvas anzeigen
C.update()
 
# Erstellung der Anfangsbestückung an Partikeln
E.setze_anzahl()
 
# Animation starten
Animationsschleife()
 
# Kontrollfenster in den Vordergrund holen
KF.lift()
 
# Lift off!
T.mainloop()

--
Dipl.-Ing. Martin Vogel
Leiter des Bauforums

Bücher:
CAD mit BricsCAD
Bauinformatik mit Python