Hallo,
Danke für Deine Antwort.
Auf diese Lösung bin ich gestern Abend gestoßen. Ich habe mich nochmal hingesetzt und habe die Fourier-Transformation angeguckt. Bei der diskreten Transformation erhält man neben den positiven Frequenzen eben auch die negativen. Es war ein Fehler in der Verwendung der FFT Funktion.
Die FFT-Funktion gibt ein Array mit einer bestimmten Sortierung zurück:
Wenn
a ein 1D Array mit
n Werten ist und
A=FFT(a) angewendet wird, dann gilt:
A[0] : Summe der Werte oder Null-Frequenz Term
A[1] bis A[n/2]: Positive Frequenz-Werte (aufsteigend sortiert)
A[n/2+1] bis A[n-1] Negative Frequenz-Werte (absteigend sortiert)
Dementsprechend habe ich immer bei der Division durch k^2 falsche Werte berechnet.
Ich habe den richtigen Code angefügt. Die Berechnung der richtigen k-Werte folgt wie immer im k-Raum:
Der Abstand zwischen zwei Gitterpunkten im k-Raum: 2 pi/Lx
Lx: Länge des gesamten Raums. Bei einer Transformation von t -> f wäre Lx die Dauer des Signals, das Transformiert werden soll. Im k-Raum kommt der Faktor 2 pi noch dazu aus der Definition der Wellenzahl bzw. Wellenvektors.
Bei
N Gitterpunkten habe ich eine x-Achse erstellt im Bereich: [-(N+1) * 2 pi/Lx : N * 2 pi/Lx]
Dasselbe für die y-Achse und dann die passenden kx^2, ky^2-Werte für das 2D-Gitter Addiert und mit der Transformierten Ladungsverteilung dividiert.
Praktischer Weise gibt es einen passenden Befehl dafür in Numpy der die passenden k-Werte erstellt. Im Code habe ich beider Varianten angegeben. Von Hand und per Numpy (Zeitlich kein Unterschied)
Code: Alles auswählen
import numpy as np
import matplotlib.pyplot as plt
length = 1 # Seitenlänge des Bereichs
steps = 100 # Anzahl der Gitterpunkte
h = length/steps # Abstand zw. Gitterpunkten
rho = np.zeros((steps,steps)) # Leeres Array für Ladungsverteilung
rho[steps//2,steps//2] = -1 # Eine Ladung in die Mitte des Feldes setzen
rho_fft = np.fft.fft2(rho) # Fourier-Transformation der Ladungsverteilung
# Variante 1: Von Hand Berechnung der k_x und k_y und Sortierung
# passend für die Transformierte Ladungsverteilung
k_positiv = np.arange(0, 2*np.pi / length * steps//2, 2*np.pi / length) #Positive k Werte im Abstand von 2pi/length
k_negativ = np.arange(2*np.pi / length, 2*np.pi / length * (steps+1)//2, 2*np.pi/length)
k = np.concatenate((k_positiv, np.flip(-k_negativ)), axis = None) # Zusammenfügen pos und neg k Werte in der passenden Reihenfolge
'''
# Variante 2: Automatisiert durch Numpy
k = 2 * np.pi * np.fft.fftfreq(steps, d=h) #Befehl eigentlich für Frequenzen, daher mit 2pi multiplizieren
'''
k_x, k_y = np.meshgrid(k, k) # kx und ky Werte für das 2D Gitter
k_sqr = k_x**2 + k_y**2 # k^2 für jeden Gitterplatz
rho_fft = np.divide(-rho_fft, k_sqr, out=np.zeros_like(rho_fft), where=k_sqr != 0) # Division jedes Gitterplatzes mit k^2 außer bei k^2=0
phi = np.fft.ifft2(rho_fft) # Rücktransformation
phi_real = np.real(phi) # Realteil der Rücktransformation
# Darstellung des Ergebnisses
plt.imshow(phi_real,extent=(0,length,0,length ))
plt.colorbar()
plt.show()
Ich meine gerade zu sehen, dass man den Algorithmus wohl noch verschnellern kann, mit der Methode die Du vorgeschlagen hast. Es kommen Doppelte Einträge für k^2 vor.