Divergenz in Python

mit matplotlib, NumPy, pandas, SciPy, SymPy und weiteren mathematischen Programmbibliotheken.
Antworten
divergence
User
Beiträge: 6
Registriert: Mittwoch 26. August 2015, 16:22

Guten Tag,

ich habe die Aufgabe, ein größeres Matlabprojekt in Python umzusetzen. Nun habe ich ein größeres Problem mit dem Berechnen der Divergenz in Numpy. Die einzelne Berechnung der Gradienten und dann das anschließende Aufsummieren klappt so weit ohne Probleme, jedoch ist die Divergenz in Matlab irgendwie anders definiert.

In Matlab gilt: div = divergence(X,Y,U,V)

wobei X und Y die beiden Koordinaten der Vektorkomponenten definieren (ich vermute im Zweifel einfach eine simple Auflistung (0,1,...,x) bzw. (0,1,...,y), bin mir jedoch nicht ganz sicher. U und V sind die jeweiligen Richtungsanteile. Hat jemand irgendeinen Ansatz, wie man das in Python umsetzen könnte? Ich sitze schon einige Stunden daran, aber habe keinerlei Idee.

Vielen Dank bereits im Voraus!
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

Hier gibt's eine gute Diskussion darüber
http://stackoverflow.com/questions/1143 ... ing-python
Beachte insbesondere den letzten Beitrag, der sagt, dass alles falsch ist, dass man die Divergenz nicht aus aufsummierten Gradienten berechnen kann.

Ich denke man baut sich die Divergenz mit numpy.diff selbst:

Code: Alles auswählen

from numpy import linspace, meshgrid, sin, cos, pi, diff
import matplotlib.pyplot as plt

# Argumentwerte als 1D Arrays erzeugen
x_1d = linspace(-2,2,801)
y_1d = linspace(-2,1,601)

# Argumentwerte als 2D Arrays erzeugen
x, y = meshgrid(x_1d, y_1d)

# Interessante Daten erzeugen
U = (x**2+y**2)*cos(pi*x)*cos(pi*y)
V = (x**2+y**2)*sin(pi*x)*sin(pi*y)

div = (diff(U,axis=1)/diff(x,axis=1))[:-1,:] + (diff(V,axis=0)/diff(y,axis=0))[:,:-1]

div_analytisch = (
    2*x*cos(pi*x)*cos(pi*y) - (x**2+y**2)*pi*sin(pi*x)*cos(pi*y) + # dU/dx
    2*y*sin(pi*x)*sin(pi*y) + (x**2+y**2)*pi*sin(pi*x)*cos(pi*y)   # dV/dy
                  )

# Plotten
plt.figure()
plt.title("U")
plt.pcolormesh(x, y, U)
plt.colorbar()
plt.gca().set_aspect("equal") 

plt.figure()
plt.title("V")
plt.pcolormesh(x, y, V)
plt.colorbar()
plt.gca().set_aspect("equal") 

plt.figure(figsize=(7,4))
plt.title("div")
plt.pcolormesh(x, y, div, vmin=-4, vmax=4)
plt.colorbar()
plt.gca().set_aspect("equal") 
plt.savefig("div.jpg", bbox_inches="tight") 

plt.figure(figsize=(7,4))
plt.title("div analytisch")
plt.pcolormesh(x, y, div_analytisch, vmin=-4, vmax=4)
plt.colorbar()
plt.gca().set_aspect("equal") 
plt.savefig("div_analytisch.jpg", bbox_inches="tight")

plt.show()
BildBild
a fool with a tool is still a fool, www.magben.de, YouTube
divergence
User
Beiträge: 6
Registriert: Mittwoch 26. August 2015, 16:22

Erstmal möchte ich mich bei dir bedanken! Total toll, der Code. Ich war richtig begeistert. Den Artikel mit dem Aufsummieren hatte ich auch gelesen. Ist ja so das einzige, was man überhaupt zu der Thematik findet.

Eine Frage, bzw. ein Problem habe ich noch. Für unsere Problemstellung ist ja die numerische Divergenz vollkommen ausreichend. Wir haben jedoch als Eingabeparameter Matrizen (48,48) jeweils für X,Y,U und V. Die Ausgabematrix der Divergenz ist jedoch in der Form (47,47). Das ist in Matlab nicht so. Habe leider den Fehler nicht gefunden. Hast du da noch eine Lösung dafür?

Vielen Dank und Liebe Grüße
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

divergence hat geschrieben:Für unsere Problemstellung ist ja die numerische Divergenz vollkommen ausreichend.
Die analytische Divergenz ist als Test gedacht. Wenn man etwas numerisch berechnet, dann macht man das i.A. deshalb, weil es nicht analytisch geht. Will man die Numerik aber testen, dann rechnet man ein Problem, wofür es eine analytische Lösung gibt und vergleicht die Ergebnisse.
divergence hat geschrieben:Wir haben jedoch als Eingabeparameter Matrizen (48,48) jeweils für X,Y,U und V. Die Ausgabematrix der Divergenz ist jedoch in der Form (47,47).
Das ist in der Tat ein Problem. Die Numpy Funktion diff macht die Arrays kürzer:
http://docs.scipy.org/doc/numpy/referen ... .diff.html
"The n order differences. The shape of the output is the same as a except along axis where the dimension is smaller by n."
Es wird bei n=1 immer die Differenz zwischen zwei Werten berechnet. Ein Array der Länge N, hat aber nur N-1 Zwischenräume.

In der Zeile 15 habe ich auch einen Hack drin: [:-1,:] und [:,:-1]
diff(U,axis=1)/diff(x,axis=1) erzeugt Werte zwischen 2 x-Werten aber an den ursprünglichen y-Werten und
diff(V,axis=0)/diff(y,axis=0) erzeugt Werte zwischen 2 y-Werten aber an den ursprünglichen x-Werten.
Um die Dimension der beiden Arrays passend zu machen habe ich den Hack mit [:-1,:] und [:,:-1] eingebaut.
Eine saubere Lösung wäre diff(U,axis=1)/diff(x,axis=1) in y-Richtung zu interpolieren und
diff(V,axis=0)/diff(y,axis=0) in x-Richtung zu interpolieren.

Wenn Du die ursprüngliche Dimensionierung haben möchtest (wie in Matlab), dann müsttest Du
diff(U,axis=1)/diff(x,axis=1) in x-Richtung zu extrapolieren und
diff(V,axis=0)/diff(y,axis=0) in y-Richtung zu extrapolieren.
Extrapolieren ist aber problematischer als Interpolieren.
a fool with a tool is still a fool, www.magben.de, YouTube
divergence
User
Beiträge: 6
Registriert: Mittwoch 26. August 2015, 16:22

Ohje, vielen Dank. Weißt du zufällig, ob Matlab da intern auch extrapoliert?

Wie kann ich denn extrapolieren? Habe jetzt spontan bei der Suche nur Wege gefunden,
eindimensionale Arrays zu interpolieren.

Vielen Dank auch weiterhin.
divergence
User
Beiträge: 6
Registriert: Mittwoch 26. August 2015, 16:22

Ganz kurz noch ein Zwischenergebnis: Ich experimentiere momentan mit scipy.interpolate.RectBivariateSpline - absolut keine Ahnung, ob ich da richtig liege, aber ich würde versuchen, damit die letzten fehlenden Punkte zu bekommen. Performancetechnisch ist das aber wahrscheinlich ne Katastrophe, oder? :-/

Mit einem simplen resize würde ich wohl auch nichts erreichen, richtig?
Benutzeravatar
MagBen
User
Beiträge: 799
Registriert: Freitag 6. Juni 2014, 05:56
Wohnort: Bremen
Kontaktdaten:

Ist es wirklich ein Problem, dass das Array in jeder Dimension 1 Element kleiner ist?
Numerisch gibt es nun einmal nur die Steigung zwischen Punkten, eine Punktsteigung bekommst Du nur, wenn Du es analytisch rechnest.
divergence hat geschrieben:Mit einem simplen resize würde ich wohl auch nichts erreichen, richtig?
Richtig, resize ist nicht mit KI ausgestattet.
a fool with a tool is still a fool, www.magben.de, YouTube
divergence
User
Beiträge: 6
Registriert: Mittwoch 26. August 2015, 16:22

Ansich würde ich sicherlich behaupten, dass es wohl auch ohne gegen würde.. Da ich hier jedoch dabei bin,
viel alten, schlecht geschriebenen Programmcode, an dem Mitarbeiter gearbeitet haben, die schon lange nicht
mehr hier arbeiten in Python neu zu schreiben, ist es leider etwas schwierig, da in der Hinsicht, den Code noch
nachzuvollziehen. Es handelt sich im Grunde um eine physikalische Simulation, bei der in jedem
Schritt ein gewisses Delta auf die Grundmatrix dazugerechnet wird (mit teilweise positiven, teilweise negativen
Werten). Hierbei kommen an vielen Stellen Gradienten vor, aber eben auch die Divergenz. Ich bin mir nicht ganz
sicher, wie ich damit arbeiten kann, wobei ich deine Begründung natürlich nachvollziehen kann. Ich frage mich,
wie es Matlab macht.

Hast du vielleicht eine Idee?

Weiterhin: Vielen Dank für deine auch bis jetzt schon sehr große Hilfe!
Antworten