A&A Title Image Startseite | Programm-Archiv | Lektion  Licht als Welle 2, py-Programm | Unserer Vorläufer | Kontakt, Datenschutz

Die Menschen drängen sich zum Lichte, nicht um besser zu sehen,
sondern um besser zu glänzen.
Friedrich Nietzsche

Das Licht als Welle, Python-Simulationsprogramm

Ich habe ein kleines Simulationsprogramm geschrieben, womit man experimentieren kann. Das vollständige Programm enthält ein paar Möglichkeiten mehr als die Kurzvariante hier. Es ist  im Programmarchiv abgelegt.

Hier die Gegenüberstellung der Programmzeilen mit den Formeln des vorigen Teils:

Parameter, Vektoren und Auslenkung (Störung)

Die Details bitte den Kommentaren entnehmen.


def sq(x):
    return x * x
 
# Main
##############################
N = 100 # Anzahl der Saitenstücke
L = 1 # Saitenlänge
h = 0.005 # Auslenkung beim Zupfen
F0 = 1 # Spannkraft
mu = 1 # Massendichte, ein knappes Gramm pro Meter
timesteps = 75
# Anzahl der Zeitschritte bis zum Programm-Ende
zp = int(N / 2) # der Zupfpunkt N/2 für mittig, N/4 für 2 Stück
initKoor(0,L,0.1,-h ,h,h / 10)
x = [0 for i in range(N + 1)] # Koordinate entlang der Saite
Y = [0 for i in range(N + 1)] # die Querauslenkung der Saite (vorige Runde)
y = [0 for i in range(N + 1)] # die Querauslenkung der Saite, berechnte Runde
v = [0 for i in range(N + 1)] # Geschwindigkeit
dx = L / N
c2 = F0 / mu # c²
dt = 0.8 * sqrt(dx * dx / c2) # Faktor 1,0 ist die Stabilitätsgrenze
print("dt=%5.3f  c=%5.3f" % (dt, sqrt(c2)))

for i in range(N + 1): # Koordinate der Saite setzen
    x[i] = i * L / N
Y[zp - 3] = h * 0.1
Y[zp - 2] = h * 0.4
Y[zp - 1] = h * 0.8
Y[zp + 0] = h * 0.9
Y[zp + 1] = h * 0.8
Y[zp + 2] = h * 0.4
Y[zp + 3] = h * 0.1

Berechnung der Kräfte und der Beschleunigung

Schleife über alle Zeitschritte, und dort für jedes Stückchen Seil:

  for j in range(timesteps): # alle Zeiten
    for i in range(1,N): # alle Masse-Elemente

Die Einzelkräfte der Seilstücke:

F1,y = F0 · (y1 - y0 ) / Δx ;   F2,y = F0 · (y2 - y1 ) / Δx
FQ = ( F2,y -  F1,y
        F1 = F0 * (Y[i] - Y[i - 1]) / dx
        F2 = F0 * (Y[i + 1] - Y[i]) / dx
        Fq = (F2 - F1) 

Die Beschleunigung benötigt die Masse des Seilstücks FQ = µ · Δx · a . Aus dieser Beschleunigung lassen sich die Orts- und die Geschwindigkeitsänderung bestimmen: 

        m = mu * dx # die infinitesimale Masse
        a = Fq / m 
        v[i] = v[i] + a * dt
        y[i] = Y[i] + v[i] * dt

Die so berechneten Koordinaten werden grafisch ausgegeben. Vorher müssen die Punkte der vorigen Runde gelöscht werden (unplot). Zum Schluss wird auf der Konsole die simulierte Zeit ausgegeben und der Grafikbildschirm aufgefrischt.

     for i in range(N): # berechnete Ergebnisse setzen
        unplot(x[i],Y[i])
        plot(x[i],y[i])
        Y[i] = y[i]
    print(j * dt)
    update() # Bildschirm anzeigen

Das war's.

--> 1 Einführung
--> 2 Funktionsweise des Programms
--> 3 Rechenergebnisse
--> 4 Sinuswelle
--> 5 Das Soliton
--> 6 Die Kugelwelle
--> 7 Interferenz zweier Wellenzüge (Formel)

Uwe Pilz, November 2020