Poisson Gl. durch Fouriertrafo mit Python lösen

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
st3fan85
User
Beiträge: 8
Registriert: Mittwoch 12. August 2020, 10:15

Hallo,
ich bin gerade dabei mir das lösen von Differentialgleichungen mit Python etwas näher zu bringen und wollte dafür die Poissongleichung mit Hilfe der diskreten Fourier-Transformation(DFT) lösen.

Bild

dazu verwende ich die DFT:

Bild

Wenn ich die DFT auf die PoissonGl anwende, dann komm ich auf das Ergebnis:

Bild

Jetzt habe ich also eine Ladungsverteilung erstellt. Zum Anfang eine einzelne Ladung in der Mitte meines Feldes. Danach wurde durch Numpy eine FFT auf das Array durchgeführt und ich habe alle Werte des Transformierten Arrays mit dem Faktor von Oben geteilt. Danach wurde die inverse FFT angewendet und sollte damit die Lösung für erhalten.
Das bekomme ich als Ergebnis raus:
Bild
Man erkennt eine nicht symmetrische Verteilung für , wobei die Lösung doch Gaußglockenförmig um der Position der Ladung sein sollte.
Hat Jemand eine Idee was ich falsch gemacht habe?

Das ist der Code den ich verwende:

Code: Alles auswählen

import numpy as np
import matplotlib.pyplot as plt

length = 10
steps = 100

rho = np.zeros((steps,steps))

rho[steps//2-1:steps//2,steps//2-1:steps//2] = -1

plt.imshow(rho, origin='lower', interpolation='none')
plt.show()

rho_fft = np.fft.fft2(rho)

for x in range(len(rho_fft)):
    for y in range(len(rho_fft)):
        if x==0 and y==0: rho_fft[x,y]=0
       
        else:
            k_sqr = (2 * np.pi * x / length)**2 + (2 * np.pi * y / length)**2
            rho_fft[x,y] = -rho_fft[x,y] / k_sqr

phi = np.fft.ifft2(rho_fft)

phi_abs = np.real(phi)

plt.imshow(phi_abs,extent=(0,length,0,length ))
plt.show() 


Vielen Dank für eure Tipps und
viele Grüße

Stefan
Benutzeravatar
__blackjack__
User
Beiträge: 12984
Registriert: Samstag 2. Juni 2018, 10:21
Wohnort: 127.0.0.1
Kontaktdaten:

@st3fan85: Keine Lösung des Problems, aber eine Anmerkung zum bisherigen Code: ``for``-Schleifen in Python und Numpy-Arrays sind in der Regel ein „code smell“. Man verwendet Numpy-Arrays ja gerade wegen der Operationen die Arrays und Numpy bieten um keine Schleifen schrieben zu müssen. Also gleiches Ergebnis mit Numpy-Mitteln:

Code: Alles auswählen

#!/usr/bin/env python3
import numpy as np
from matplotlib import pyplot as plt


def main():
    length = 10
    steps = 100

    rho = np.zeros((steps, steps))

    rho[steps // 2 - 1 : steps // 2, steps // 2 - 1 : steps // 2] = -1

    plt.imshow(rho, origin="lower", interpolation="none")
    plt.show()

    rho_fft = np.fft.fft2(rho)

    xx, yy = np.meshgrid(range(steps), range(steps))
    k_sqr = (2 * np.pi * xx / length) ** 2 + (2 * np.pi * yy / length) ** 2
    rho_fft = np.divide(
        -rho_fft, k_sqr, out=np.zeros_like(rho_fft), where=k_sqr != 0
    )

    phi = np.fft.ifft2(rho_fft)
    phi_abs = np.real(phi)
    plt.imshow(phi_abs, extent=(0, length, 0, length))
    plt.show()


if __name__ == "__main__":
    main()
“Most people find the concept of programming obvious, but the doing impossible.” — Alan J. Perlis
tonikae
User
Beiträge: 90
Registriert: Sonntag 23. Februar 2020, 10:27

Ich weiss nicht ob ich dich richtig verstehe, aber die Poisson-Verteilung ist ja eigentlich nie identisch mit der
Gauss'schen Normalverteilung. Also wären die Asymetrien eigentlich normal.

https://drive.google.com/file/d/1AkyWb1 ... k4ILv/view
https://drive.google.com/file/d/1lhj6DV ... hLSOS/view
st3fan85
User
Beiträge: 8
Registriert: Mittwoch 12. August 2020, 10:15

Hallo,
erstmal vielen Dank für eure Antworten.
@tonikae
Ich gebe Dir da vollkommen recht, aber du sprichst von der Poisson-Verteilung aus dem Gebiet Wahrscheinlichkeits-Theorie. Das ist tatsächlich derselbe Herr Poisson. Sein Name tritt allerdings auch in Verbindung mit der oben genannten Poisson-Gleichung auf. Diese wird z.B. in der Elektrodynamik in der Physik verwendet. Damit lässt sich ein Elektrisches Potential einer gegebenen Ladungsverteilung berechnen. In meinem Beispiel befindet sich in einem Raumbereich ein geladenen Teilchen z.B ein Elektron. Die Lösung sollte eigentlich in dieser Art aussehen:
Bild

@__blackjack__
Danke für den Tipp! Den werde ich jetzt immer beherzigen.
Zeitunterschied für die Division:
Mit for-Schleife: 0.016994 s
Mit np.divide: 0.000982 s


Ich habe mir die Fourier-Transformation der Punktladung plotten lassen. In der Version von oben hat die Punktladung eine Größe von 2x2 Feldpunkten und es kam ein merkwürdiges Bild raus:
Bild

Nachdem ich die Größe der Punktladung vergrößert habe, sieht es gut aus:
rho[steps//2-2:steps//2+2,steps//2-2:steps//2+2] = -1
Bild

Aber leider sieht die Lösung nach der Rücktransformation genauso aus wie vorher:
Bild

Schöne Grüße
Stefan
__deets__
User
Beiträge: 14480
Registriert: Mittwoch 14. Oktober 2015, 14:29

Zur Physik und Herleitung kann ich mangels Grundlagen nichts sagen. Aber da ich versuche, DSP zu lernen, und da die FFT ja zentral ist, fällt mir was auf: du behandelst deine fft2 so, als ob das in jeder Richtung steps Frequenzen sind. Das stimmt aber nicht! Nyquist lehrt uns, das wir nur steps/2 Frequenzen haben können, danach aliasen wir. Trotzdem hat die FFT steps viele komplexe Einträge. Und das liegt AFAIK daran, das eine reellwertige Frequenz aus zwei Phasoren gebildet wird. Der cos f zb aus der Summe zweier komplex konjugierter. Und die liegen üblicherweise an den Stellen i und i + steps.

Daher wendest du deine Koeffizienten asymmetrisch an. Ich denke du musst die auf steps//2 beschränken, und dann 4 mal anwenden, auf (0, 0), (0, steps), (steps, 0) und (steps, steps). Ob das schon reicht, oder du ggf noch die Wurzel ziehen musst, weil sonst dein K quadratisch eingeht, musst du mal überlegen.
st3fan85
User
Beiträge: 8
Registriert: Mittwoch 12. August 2020, 10:15

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.
Antworten