Frage zu 3D Plot und zugehöriger for-loop

[kowi]

Hallo zusammen,

ich bin neu hier im Forum und auch in Sachen Python. Aufgrund meiner Abschlussarbeit habe ich mich (so gut es in der kurzen Zeit ging) in die Sprache eingearbeitet, dennoch hapert es oft an Kleinigkeiten. Ich bin mir sicher, dass das nicht mein einziger Post bleiben wird Trotz Benutzen der Suchfunktion konnte ich bisher keine Lösung finden, daher hoffe ich, dass ihr mir auf die Sprünge helfen könnt. Ich habe bisher folgenden Code erstellt:

Code: Alles auswählen

 from mpl_toolkits.mplot3d import Axes3D 
 from matplotlib import cm
 from math import *
 from scipy.special import *
 import matplotlib.pyplot as plt
 import numpy as np

 ## Definition der Parameter für Druckgleichung nach Rudnicki (1986) ##
 q = 6.0/1000                                                                       
 lameu = 11.2*10**9                              
 lame = 8.4*10**9                               
 pi                                            
 alpha = 0.65                                    
 G = 8.4*10**9                                   
 k = 1.0e-15                                      
 eta = 0.001                                                 
 t = 1000*365*24600        

 kappa = k/eta                                                    
 print "kappa ist:",kappa                                        
 c = ((kappa*(lameu-lame)*(lame+2*G))/((alpha**2)*(lameu+2*G)))  
 print "c ist:",c                                                


 xmin = -10
 xmax = 10
 ymin = -10
 ymax = 10
 for x in range (xmin,xmax):
      for y in range (ymin,ymax):
          r=sqrt(x**2+y**2)
          P=(q/(rhof*4*pi*kappa))*(expn(1,r**2/(4*c*t)))
          z = P/1e6
          print x,  y,  z


 x, y = np.meshgrid(x, y)


 ## Plotting in 3D ##
 fig = plt.figure()
 ax = fig.gca(projection='3d')
 surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0,
                   antialiased=False, vmin=np.nanmin(z), vmax=np.nanmax(z))
 fig.colorbar(surf, shrink=0.5, aspect=5)

 ## Achsenskalen ##
 ax.set_xlim(xmin,xmax)      # x-Achsenskala vorgeben
 ax.set_ylim(ymin,ymax)      # y-Achsenskala vorgeben

 ## Beschriftung der Achsen ##
 ax.set_title('Druckverteilung')
 ax.set_xlabel('Distanz zu Well [m]')
 ax.set_ylabel('Distanz zu Well [m]')
 ax.set_zlabel('Druck in [MPa]')

 plt.show()  

Davor habe ich eine andere Variante mit np.arange statt der for-Schleife erstellt, da hat es mit dem 3D Plot geklappt. Da meine Werte für P bzw. z] nahe null gegen unendlich gehen, konnte ich das Problem mittels z[z==np.inf] = np.nan umgehen. Diesen Schritt kann ich hier ja nicht anwenden. Generell habe ich daher 2 Fragen:

- habe ich die for-Schleife richtig benutzt? Soweit ich weiß, muss man teilweise statt P z.B. P[x,y] coden, um eine Abhängigkeit dieser Werte zu definieren? Kann mir da jemand ein Beispiel geben?

- was muss ich tun, damit ich meinen 3D-Plot vernünftig hinbekomme? Mittels numpy und np.arange habe ich da bisher erfolg gehabt, aber hier stehe ich auf dem Schlauch. Muss ich vorher leere Felder à la np.zeros erzeugen? Oder gehören die gar nicht hierher? Auch hier wäre ein kurzes Beispiel echt klasse.

So, ein langer erster Post. Diese Fragen habe ich zudem bereits auf stackoverflow.com gestellt, die Antwort (nur Bahnhof verstanden) aber leider aufgrund meiner geringen Kenntnisse nicht nachvollziehen können. Hoffe auf Verständnis und Hilfe

Grüße
Alex

P.S: falls jemand nachschauen möchte, hier der Originalpost: http://stackoverflow.com/questions/2721 ... nge-python

[MagBen]

Die Numpy Funktion die Du brauchst ist nicht arange, sondern linspace:

Code: Alles auswählen

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from scipy.special import *
import matplotlib.pyplot as plt
import numpy as np

q = 6.0/1000                                                                      
lameu = 11.2*10**9                              
lame = 8.4*10**9                              
alpha = 0.65                                    
G = 8.4*10**9                                  
k = 1.0e-15                                      
eta = 0.001                                                
t = 1000*365*24600        
rhof = 8.314 # <==== geraten, ist so richtig?
kappa = k/eta                                                    
c = ((kappa*(lameu-lame)*(lame+2*G))/((alpha**2)*(lameu+2*G)))  

x = np.linspace(-10,10,21)
y = np.linspace(-10,10,21)
x, y = np.meshgrid(x, y)
r=np.sqrt(x**2+y**2)
P=(q/(rhof*4*np.pi*kappa))*(expn(1,r**2/(4*c*t)))
z = P/1e6

## Plotting in 3D ##
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0, vmin=700, vmax=1100)
fig.colorbar(surf, shrink=0.5, aspect=5)

## Beschriftung der Achsen ##
ax.set_title('Druckverteilung')
ax.set_xlabel('Distanz zu Well [m]')
ax.set_ylabel('Distanz zu Well [m]')
ax.set_zlabel('Druck in [MPa]')

plt.show()  

Mit np.linspace(-10,10,201) wird der Plot glatter, braucht aber auch deutlich länger.

- habe ich die for-Schleife richtig benutzt? Soweit ich weiß, muss man teilweise statt P z.B. P[x,y] coden, um eine Abhängigkeit dieser Werte zu definieren? Kann mir da jemand ein Beispiel geben?

Auf diese Art von Schleifen kann und soll man man weitgehend verzichten, wenn man Numpy benutzt. Die Schleifen gibt es jetzt natürlich immer noch, aber sie laufen ca. 100 mal Schneller im C-Code von Numpy ab.

Die Grenzen (vmin=700, vmax=1100) habe ich von Hand eingefügt, die Farbberechnung des Plots kommt mit der Polstelle sonst nicht gut zurecht, alles wird gleichmäßig Blau. Ansonsten brauchst Du Dir keine Sorgen bezüglich Nulldivision, np.inf und np.nan. Numpy und Matplotlib kommt damit sehr gut zurecht.

Benutze nicht das Paket math zusammen mit Numpy und schon gar nicht

Code: Alles auswählen

from math import *

Alles aus math findest Du auch bei Numpy, die Funktionen von math funktionieren aber nicht mit Numpy-Arrays.

Das ist in unnötig:

Code: Alles auswählen

#ax.set_xlim(xmin,xmax)      # x-Achsenskala vorgeben
#ax.set_ylim(ymin,ymax)      # y-Achsenskala vorgeben

Wenn Du es weglässt, macht Matplotlib genauso automatisch.

rhof hat gefehlt, soll das die gaskonstante sein?

Code: Alles auswählen

rhof = 8.314

[kowi]

Hi MagBen,

erstmal vielen Dank für dein Feedback! Ich habe die for-schleife nur benutzt, da meine Betreuerin meinte, dass man damit besser die Werte für P kontrollieren könne, wenn man sie sich ausgeben lässt. Dies ist in einem Array als Ausgabe ja nur bedingt möglich, oder gibt es Wege sich die Werte "komplett" anzeigen zu lassen? Wenn ich mir den P-array ausgeben lasse, bekomme ich eben immer z.B.
[[Wert, Wert,..., Wert, Wert]]. Das war der einzige Hintergrund, sozusagen als Kontrolle für meine Werte.

Ich habe zu deinen Erläuterungen noch Fragen, und zwar:

- mir ist der Unterschied zwischen np.arange und np.linspace bekannt, aber wieso ist hier die Benutzung von zweitem besser geeignet?
- du hast die Grenzen (vmin=700, vmax=1100) von Hand eingefügt. Wie kommst du auf die Werte?
- math habe ich importiert, um mir den exakten Wert für Pi zu holen. Sollte ich Pi lieber so definieren?

Danke übrigens für den Hinweis, rhof habe ich in der Tat vergessen. Die Variable stellt die Dichte meines injizierten Fluides dar, in dem Fall die von Wasser (1000 kg/m³).

Ich poste abschließend hier nochmals den Code, in dem ich np.arange verwendet habe:

Code: Alles auswählen

from mpl_toolkits.mplot3d import Axes3D 
from matplotlib import cm
from math import *
from scipy.special import *
import matplotlib.pyplot as plt
import numpy as np

#Definition der Parameter für Druckgleichung nach Rudnicki (1986)
q = 6.0/1000                   #Fluidmenge pro Fläche und Zeit [m³/s]
rhof = 1000                     #Dichte Flüssigkeit [kg/m³]
lameu = 11.2*10**9         #Lamé-Parameter, undrained [GPa]
lame = 8.4*10**9            #Lamé-Parameter, drained [GPa]
pi                                  #durch Pythonmodul "math" gegeben
alpha = 0.65                    #Biot-Willis-Koeffizient
G = 8.4*10**9                 #Schermodul [GPa]
k = 1.0*10**(-15)             #Permeabilität [m²] bzw. [Darcy] 
eta = 0.001                      #Viskosität des Fluids [Pa*s]            
t = 1000*24*3600               #Zeit in [s]

#Beziehungen der Parameter untereinander
kappa = k/eta                                                   #Berechnung der Permeabilität nach Rudnicki (1986), [m³*s/kg]
print "kappa ist:",kappa                                        #Ausgabe Permabilität
c = ((kappa*(lameu-lame)*(lame+2*G))/((alpha**2)*(lameu+2*G)))  #Berechnung der Diffusivität
print "c ist:",c                                                #Ausgabe der Diffusivität
    
#Bereiche der Achsen in [m]
xmin = -10
xmax = 10
ymin = -10
ymax = 10

#Wertebereich
x = np.arange(xmin,xmax,0.1)
x[(x >= -0.2) & (x <= 0.1)] = 0
y = np.arange(ymin,ymax,0.1)
y[(y > -0.2) & (y < 0.1)] = 0
x, y = np.meshgrid(x, y)

#Formeln für den Druck
r = np.sqrt(x**2+y**2)
P = (q/(rhof*4*pi*kappa)*(expn(1,(r**2)/(4*c*t))))      
z = P/1e6                                                                     
print z
z[z==np.inf] = np.nan

#Plotting in 3D
fig = plt.figure()
ax = fig.gca(projection='3d')
surf = ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap=cm.jet, linewidth=0,
                       antialiased=False, vmin=np.nanmin(z), vmax=np.nanmax(z))
fig.colorbar(surf, shrink=0.5, aspect=5)

#Achsenskalen
ax.set_xlim(xmin,xmax)      # x-Achsenskala vorgeben
ax.set_ylim(ymin,ymax)      # y-Achsenskala vorgeben

#Beschriftung der Achsen
ax.set_title('Druckverteilung')
ax.set_xlabel('Distanz zu Well [m]')
ax.set_ylabel('Distanz zu Well [m]')
ax.set_zlabel('Druck in [MPa]')

plt.show()

Hier habe ich beispielsweise auf Hinweis von einem User von stackoverflow folgende Zusätze benutzt:

Code: Alles auswählen

z[z==np.inf] = np.nan

und

Code: Alles auswählen

vmin=np.nanmin(z), vmax=np.nanmax(z)

.

Davor sah mein Plot auch nur einfarbig blau aus.

Viele Grüße,
Alex

[BlackJack]

Bezüglich `pi`: Nein, nicht selber definieren, aber:

Code: Alles auswählen

In [10]: numpy.pi, math.pi, numpy.pi - math.pi == 0
Out[10]: (3.141592653589793, 3.141592653589793, True)

Und Du wirst wohl mit keinem Modul den *exakten* Dezimalwert von `pi` bekommen.

[MagBen]

mir ist der Unterschied zwischen np.arange und np.linspace bekannt, aber wieso ist hier die Benutzung von zweitem besser geeignet?

Beim Plotten muss man bei arange mehr Kopfrechnen machen, z.B.

Code: Alles auswählen

x = np.arange(xmin,xmax,0.1)

geht nämlich nicht bis xmax, sondern nur bis xmax-0.1. Außerdem siehst Du nicht sofort wieviele Werte erzeugt werden und die Anzahl ist ein wichtiger Parameter beim Plotten für die Optik und die Rechenzeit.