Licht als Welle

--> Lösungen der Wellengleichung

Seilwellen

An Hand einer Seilwelle werden wichtige Eigenschaften des Lichts mit einem Simulationsprogramm zugänglich gemacht. Dieses Programm gehort zur Lektion Licht als Welle . Im selben Verzeichnis muss die Programmsammlung plot liegen.

# Wellengleichung - Seilwelle

from math import *
from turtle import *
from random import *
from plot import *
import sys

def plot(x,y):
    Plot(x,y,3, 'black')
    
def unplot(x,y):
    Plot(x,y,3, 'white')

def sq(x):
    return x * x
 
# Main
##############################
# Die Modi:
# 1 - schmaler Störung (Dirac-Impuls) in der Mitte
# 2 - wellenförmige Störung in der Mitte
# 3 - zwei wellenförmige Störungen
# negativ: wie oben, aber zwei lineare Dichten
modus = 3
N = 100 # Anzahl der Saitenstücke
L = 1 # Saitenlänge
h = 0.05 # Auslenkung beim Zupfen
F0 = 1 # Spannkraft
mu = 1 # Massendichte, ein knappes Gramm pro Meter
timesteps = 100 # Anzahl der Zeitschritte bis zum Programm-Ende
if abs(modus) < 3:
    zp = int(N / 2) # der Zupfpunkt N/2 für mittig, N/4 für 2 Stück
    zp2 = zp
else:
    zp = int(N / 4) # der Zupfpunkt N/2 für mittig, N/4 für 2 Stück
    zp2 = int(3 * N / 4) # der Zupfpunkt
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
if modus < 0:
    dt=dt/2 # dt ist für den schnellen Abschnitt sonst zu hoch
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
# 1: eine Dirac-Auslenkung führt zu Rauschen in der Mitte
if (1 == abs(modus)):
    Y[zp] = h
else: # wellenförmige Störung
    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] = Y[zp - 1]
    Y[zp + 2] = Y[zp - 2]
    Y[zp + 3] = Y[zp - 3]
    # zweite Störung, wirkt nur bei Modus 3
    Y[zp2 - 3] = Y[zp - 3]
    Y[zp2 - 2] = Y[zp - 2]
    Y[zp2 - 1] = Y[zp - 1]
    Y[zp2 + 0] = Y[zp + 0]
    Y[zp2 + 1] = Y[zp - 1]
    Y[zp2 + 2] = Y[zp - 2]
    Y[zp2 + 3] = Y[zp - 3]


# Ausgangssituation plotten
for i in range(N): plot(x[i],Y[i])
update() # Bildschirm anzeigen
print("Zeit:")
# rechnen
for j in range(timesteps): #
    if (j % 20 == 10): # die vielen Objekte wieder loswerden
        initKoor(0,L,0.1,-2 * h, 2 * h, 2 * h / 10)
    for i in range(1,N):
        # die beiden Winkel links und rechts des Massestückes
        F1 = F0 * (Y[i] - Y[i - 1]) / dx
        F2 = F0 * (Y[i + 1] - Y[i]) / dx
        Fq = (F2 - F1) 
        m = mu * dx # die infinitesimale Masse
        # Massenkorrektur für negative Modi
        if modus < 0:
            if i < N / 3:
                m = 2 * m 
            if i > 2 * N / 3:
                m = 0.5 * m
        a = Fq / m 
        v[i] = v[i] + a * dt
        y[i] = Y[i] + v[i] * dt
    for i in range(N): # berechnete Ergebnisse setzen
        unplot(x[i],Y[i])
        plot(x[i],y[i])
        Y[i] = y[i]
    print("%5.3f " % (j * dt),end="")
    sys.stdout.flush()
    update() # Bildschirm anzeigen
hideturtle()
name = input("ENDE!")

Lichtwellen

# Lichtwelle
# Elliptische Differentialgleichung - Lichtwelle durch aufgeprägem Sinus u.ä.
from math import *
from turtle import *
from random import *
from plot import *

def plot(x,y):
    Plot(x,y,3, 'black')
    
def unplot(x,y):
    Plot(x,y,3, 'white')

def sq(x):
    return x * x
 
def initK(L, h):
    initKoor(0,L,0.1,-2 * h ,2 * h,h / 2)

# Main
##############################
modus = 5
N = 200 # Anzahl der Saitenstücke
L = 1 # Saitenlänge
EMax = 0.01 # max.  Feldstärke
timesteps = 250
# 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
initK(L,EMax)
x = [0 for i in range(N + 1)] # Koordinate entlang der Saite
E = [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
c = 1 
c2 = c * 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


# Ausgangssituation plotten
for i in range(N): plot(x[i],E[i])
update() # Bildschirm anzeigen

# rechnen
for t in range(timesteps): #
    # primitiver Dirac-Impuls
    if modus == 0:
        if (t==0):
            E[0] = EMax 
    # Dirac-Impuls
    if modus == 1:
        tLokal = tLokal = (t - 25) / 5
        sigma = 0.2
        E[0] = EMax * sigma / (sigma + tLokal * tLokal)
    # Sinuswelle
    if modus == 2:
        T = 100
        E[0] = EMax * sin(2 * pi * (t / T)) 
    # ansteigende Sinuswelle
    if modus == 3:
        E[0] = t / 100 * EMax * sin(t / (4 * pi)) 
    # abschwellender Sinuswelle
    if modus == 4:
        E[0] = 2 * EMax / (1 + t / 100) * sin(t / (4 * pi)) 
    # Soliton, die stabilste aller Wellen, keine Dispersion
    if modus == 5:
        if (t < 50):
            tLokal = (t - 25) / 5
            E[0] = 3 * EMax / (exp(-tLokal) + exp(tLokal))
    # Wellenpaket
    if modus == 6:
        if (t < 100):
            tLokal = (t - 50) / 10
            E[0] = 3 * cos(3 * tLokal) * EMax / (exp(-tLokal) + exp(tLokal)) # Kosinus, damit es symmetrisch zu t=0 ist
    if (t % 20 == 10): # die vielen Objekte wieder loswerden
        initK(L,EMax)
    for i in range(1,N):
        # die beiden Winkel links und rechts des Massestückes
        a1 = c / (dx * dx) * (E[i] - E[i - 1]) 
        a2 = c / (dx * dx) * (E[i + 1] - E[i]) 
        a = a2 - a1
        v[i] = v[i] + a * dt
        y[i] = E[i] + v[i] * dt
    for i in range(N): # berechnete Ergebnisse setzen
        unplot(x[i],E[i])
        plot(x[i],y[i])
        E[i] = y[i]
    E[N] = y[N - 1]
    update() # Bildschirm anzeigen

Uwe Pilz, Dezember 2020