RegistrierenRegistrieren   LoginLogin   FAQFAQ    SuchenSuchen   
Poisson Gl. durch Fouriertrafo mit Python lösen
 
Neue Frage »
Antworten »
    Foren-Übersicht -> Elektrik
Autor Nachricht
St3fan85



Anmeldungsdatum: 04.05.2015
Beiträge: 18

Beitrag St3fan85 Verfasst am: 30. Apr 2021 13:05    Titel: Poisson Gl. durch Fouriertrafo mit Python lösen Antworten mit Zitat

Hallo,
ich bin gerade dabei mir das Programmieren mit Python beizubringen und wollte dafür die Poissongleichung mit Hilfe der diskreten Fourier-Transformation(DFT) lösen.



dazu verwende ich die DFT:


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



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 Ergebnis habe ich als Plot hochgeladen. 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:

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



poisson_with_Fourier.png
 Beschreibung:
 Dateigröße:  22.9 KB
 Angeschaut:  647 mal

poisson_with_Fourier.png




Zuletzt bearbeitet von St3fan85 am 03. Mai 2021 23:51, insgesamt einmal bearbeitet
St3fan85



Anmeldungsdatum: 04.05.2015
Beiträge: 18

Beitrag St3fan85 Verfasst am: 02. Mai 2021 22:21    Titel: Antworten mit Zitat

Moin,
ich habe mir jetzt die Fourier-Transformierte Ladungsverteilung vom ersten Post plotten lassen und habe als Ergebnis folgendes bekommen:
https://i.ibb.co/wKDG12w/Alt-verteilung.png

Nachdem ich die Größe des Punktladung (2x2 Feldpunkte) vergrößert habe sieht das Ergebnis jetzt gut aus:
rho[steps//2-2:steps//2+2,steps//2-2:steps//2+2] = -1

https://i.ibb.co/Jx1dnGv/FFT-Punktladung.png

Leider ist aber das Ergebnis nach der Rücktransformation immer noch dasselbe...

https://i.ibb.co/MDLChz1/poisson-unsym.png

Ich versteh nicht so richtig, woran das liegen kann.
TomS
Moderator


Anmeldungsdatum: 20.03.2009
Beiträge: 14391

Beitrag TomS Verfasst am: 03. Mai 2021 08:35    Titel: Antworten mit Zitat

In deiner Poisson-Gleichung fehlt ein Quadrat beim Nabla-Operator.

Zu deinem Code kann ich erst mal nichts sagen.

So wie du die Gleichung ansetzt, löst du das Problem auf einem 2-Torus bzw. der 2-dim. Ebene mit der Zusatzbedingungen periodischer Funktionen in x und y. Diese Bedingung muss auch deine Ladungsverteilung erfüllen.

Hast du das Problem mal in einer Dimension gelöst?

Ich fasse nochmal zusammen:





Daraus folgt



und somit



Dabei ist zu beachten, dass konstante und lineare Funktionen in x,y durch diesen Ansatz nicht festgelegt werden, d.h. auch, dass m=0, n=0 gesondert betrachtet werden müssen.

Soweit ist jedenfalls alles gut.

_________________
Niels Bohr brainwashed a whole generation of theorists into thinking that the job (interpreting quantum theory) was done 50 years ago.
TomS
Moderator


Anmeldungsdatum: 20.03.2009
Beiträge: 14391

Beitrag TomS Verfasst am: 03. Mai 2021 09:19    Titel: Re: Poisson Gl. durch Fouriertrafo mit Python lösen Antworten mit Zitat

St3fan85 hat Folgendes geschrieben:



Das ist der Code den ich verwende:
Code:

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

in deiner Formel stehen richtigerweise x, y sowie k, l; in deinem Code sehe ich dagegen

Code:

k_sqr = (2 * np.pi * x / length)**2 + (2 * np.pi * y / length)**2


kein k, l.

Ich kenne die Methode der FFT in Python nicht, aber irgendwie sieht das seltsam aus. Wieso packst du in der Schleife die Werte in ein Array in [x,y]? Du rechnest doch im k-Raum, also erwarte ich ein Array in [k,l], anschließend die Rücktrf. nach [x,y].

Aus eigener Erfahrung: verwende konsistente Bezeichner, Variablen etc. zwischen deinen Rechnungen und dem Code.

_________________
Niels Bohr brainwashed a whole generation of theorists into thinking that the job (interpreting quantum theory) was done 50 years ago.
St3fan85



Anmeldungsdatum: 04.05.2015
Beiträge: 18

Beitrag St3fan85 Verfasst am: 03. Mai 2021 18:10    Titel: Antworten mit Zitat

Hallo,
Das Quadrat am Nabla habe ich vergessen.
Mit den Bezeichnungen hast Du natürlich vollkommen Recht, das führt nur zu Verwirrungen.
Der Code wurde allerdings auch schon verändert, damit nur numpy-Funktionen zu verwenden, um Schleifen zu vermeiden. Das bringt eine deutlich bessere Effizienz:
Code:

import numpy as np
import matplotlib.pyplot as plt

length = 10
steps = 100

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

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

plt.imshow(rho, origin='lower', interpolation='none')
plt.show()
rho_fft = np.fft.fft2(rho)
rho_fft_shift = np.fft.fftshift(rho_fft)
plt.imshow(np.abs(rho_fft_shift))
plt.show()
   
kk, ll = np.meshgrid(range(steps), range(steps))

k_sqr = (2 * np.pi * kk / length) ** 2 + (2 * np.pi * ll / 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()


Ich habe mir das irgendwie zu einfach vorgesellt. Ich versteh auch nicht genau wie ich die Randbedingungen mit einbringen kann.

Ich habe gelesen den Laplace-Op zu diskretisieren und dann die Fourier-Transformation anzuwenden:



Das ist dabei der Abstand der diskretisierungs-Punkten.
Wenn die diskrete Fourier trafo in die `s eingesetzt wird, heben sich fast alle Terme aus den Summen auf und man kommt auf die Gleichung:


mit


Diese Methode liefert tatsächlich die richtige Lösung der Poisson-Gl. für eine Punktladung:
https://i.ibb.co/VLDXnGt/Poisson-Lsg.png

Ich verstehe nicht was der Unterschied zwischen den beiden Varianten ist.
TomS
Moderator


Anmeldungsdatum: 20.03.2009
Beiträge: 14391

Beitrag TomS Verfasst am: 03. Mai 2021 22:22    Titel: Antworten mit Zitat

Jetzt musst du erst mal erklären, worauf du hinauswillst.

Willst du die kontinuierliche Poisson-Gleichung in einem endlichen Intervall bzw. kompaktifizierten Bereich lösen? Das funktioniert mittels der Fouriertransformation so wie oben von dir eingeführt und von mir nochmals bestätigt. Die Fouriertransformation erzeugt aus dem kontinuierlichen Problem im Ortsraum (x,y) ein diskretes Problem im Impulsraum (k,l). Dazu musst du keine Ableitung diskretisieren, du berechnest sie ja exakt.

Oder willst du die kontinuierliche Poisson-Gleichung im Ortsraum diskretisieren? Wozu benötigst dann noch die Fouriertransformation?

Irgendwie vermischt du jetzt beides. Warum?

_________________
Niels Bohr brainwashed a whole generation of theorists into thinking that the job (interpreting quantum theory) was done 50 years ago.
St3fan85



Anmeldungsdatum: 04.05.2015
Beiträge: 18

Beitrag St3fan85 Verfasst am: 03. Mai 2021 22:47    Titel: Antworten mit Zitat

Mein Ziel ist es die Poisson-Gleichung numerisch im Ortsraum zu lösen. Dies ist, wenn ich mich nicht irre, nur möglich auf einem begrenzten Raum der diskretisiert ist. Ich habe schon einmal die Poisson-Gl. mit Dirichlet-RB mit Hilfe der Finiten-Differenzen-Methode numerisch gelöst. Da aber diese Methode rechnerisch sehr aufwendig ist und daher viel Rechenzeit braucht, wollte ich es über den Weg mit der Fourier-Transformation versuchen, da dieser als sehr effizient gilt.
TomS
Moderator


Anmeldungsdatum: 20.03.2009
Beiträge: 14391

Beitrag TomS Verfasst am: 03. Mai 2021 23:18    Titel: Antworten mit Zitat

St3fan85 hat Folgendes geschrieben:
Mein Ziel ist es die Poisson-Gleichung numerisch im Ortsraum zu lösen. Dies ist, wenn ich mich nicht irre, nur möglich auf einem begrenzten Raum der diskretisiert ist.

Ja.

Die Bedingung des endlichen Raumes kann man durch eine geeignete Koordinatentransformation umgehen, allerdings erkauft läuft man damit in andere Schwierigkeiten.

Also bleiben wir zunächst beim endlichen Raum.

St3fan85 hat Folgendes geschrieben:
Da aber diese Methode rechnerisch sehr aufwendig ist und daher viel Rechenzeit braucht, wollte ich es über den Weg mit der Fourier-Transformation versuchen, da dieser als sehr effizient gilt.

Auch gut.

Allerdings, wenn du die Fourier-Transformation nutzt, rechnest du nicht mehr im kontinuierlichen Ortsraum, stattdessen verwendest du diskrete Vektoren und Matrizen in k,l. Letztlich steht alles, was mathematisch notwendig ist, bereits da, du musst es lediglich korrekt implementieren.

_________________
Niels Bohr brainwashed a whole generation of theorists into thinking that the job (interpreting quantum theory) was done 50 years ago.
TomS
Moderator


Anmeldungsdatum: 20.03.2009
Beiträge: 14391

Beitrag TomS Verfasst am: 03. Mai 2021 23:58    Titel: Antworten mit Zitat

Ich habe mir jetzt nochmal deinen Code angesehen. In meiner ersten Interpretation steckte ein Denkfehler, das ist jetzt klar. Abgesehen von den Zeilen, die ich mangels Python-Kenntnissen nicht ganz verstehe, sieht das doch sinnvoll aus.

Trotzdem: der Fehler steckt im Code; die Mathematik im ersten Ansatz ist ok und sinnvoll.

_________________
Niels Bohr brainwashed a whole generation of theorists into thinking that the job (interpreting quantum theory) was done 50 years ago.
St3fan85



Anmeldungsdatum: 04.05.2015
Beiträge: 18

Beitrag St3fan85 Verfasst am: 04. Mai 2021 13:02    Titel: Antworten mit Zitat

Hallo,
ich habe herausgefunden wo das Problem liegt. 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 falsche Werte berechnet.

Ich habe den richtigen Code angefügt. Die Berechnung der richtigen -Werte folgt wie immer im -Raum:

Der Abstand zwischen zwei Gitterpunkten im -Raum:
Bei Gitterpunkten habe ich eine x-Achse erstellt im Bereich
Das Selbe für die y-Achse und dann die passenden -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


Code:

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)   
k_negativ = np.arange(2*np.pi / length, 2*np.pi / length * (steps+1)//2, 2*np.pi/length)
# Zusammenfügen pos und neg k Werte in der passenden Reihenfolge
k = np.concatenate((k_positiv, np.flip(-k_negativ)), axis = None) 
k_x, k_y = np.meshgrid(k, k)  # kx und ky Werte für das 2D Gitter

'''
# Variante 2: Automatisiert durch Numpy   
#Befehl eigentlich für Frequenzen, daher mit 2pi multiplizieren
k = 2 * np.pi * np.fft.fftfreq(steps, d=h)     
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

# Division jedes Gitterplatzes mit k^2 außer bei k^2=0
rho_fft = np.divide(-rho_fft, k_sqr, out=np.zeros_like(rho_fft), where=k_sqr != 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()


Lösung:
https://i.ibb.co/1s2Hqd6/pois-lsg.png

Vielen Dank für die tolle Hilfe
TomS
Moderator


Anmeldungsdatum: 20.03.2009
Beiträge: 14391

Beitrag TomS Verfasst am: 04. Mai 2021 16:41    Titel: Antworten mit Zitat

Gerne - und danke für deine Rückmeldung und die Ergebnisse.
_________________
Niels Bohr brainwashed a whole generation of theorists into thinking that the job (interpreting quantum theory) was done 50 years ago.
Neue Frage »
Antworten »
    Foren-Übersicht -> Elektrik