3D-Plot

[sliderchrisss]

Hallo,

habe ein Gitternetz(x-y-Ebene) Auf diesem Gitternetz habe ich nun zu jedem Gitterpunkt den entsprechenden Wert in einer Matrix abgespeichert ( Es geht um eine Temperaturverteilung).
Nun möchte ich daraus einen 3-D Plot machen, also die Temperaturverteilung über diese Fläche grafisch darstellen. D.h. jedem Gitterpunkt (i*x,j*x) den entrechenden Wert T[i,j] zuweisen.
Habe es schon mit mplot3d versucht, aber da lassen sich nur Skalare als Funktionswerte eingeben...und keine Matrixeinträge, so wie ich es gerne hätte.
...hoffentlich ist mein Problem verständlich.

Danke im Voraus.

[problembär]

Habe es schon mit mplot3d versucht, aber da lassen sich nur Skalare als Funktionswerte eingeben...und keine Matrixeinträge, so wie ich es gerne hätte.

Kann man das nicht umrechnen?

[gkuhl]

Hallo sliderchrisss und Willkommen im Forum,

sowas sollte eigentlich mit mplot3d gehen. Hier gibt es einige Beispiele wie. Zeig am besten mal etwas Code, damit wir dein Problem nachvollziehen können.

Grüße
Gerrit



PS: Für Temperaturfelder finde ich 2d-Darstellungen (Beispiele) meist schöner.

[xpilz]

Hier gibt es einige Beispiele wie.

Ich glaube der Link ist kaputt . Edit: Nein nur ein Tippfehler am Ende das l weg dann gehts.

@sliderchrisss

Auf diesem Gitternetz habe ich nun zu jedem Gitterpunkt den entsprechenden Wert in einer Matrix abgespeichert

Kannst du das nicht umgehen? Ich meine einzeln abspeichern. Wäre zwar viel Aufwand und würde glaub ich auch sehr unpythonisch sein, aber dann funktionierts doch mit mplot3d. Das würde ich aber nur tun wenn es gar keine andere Möglichkeit gibt, was ich bezweifle.

[sliderchrisss]

Hier ist der Code.
Vorweg: Ich bin absoluter Einsteiger, was das Programmieren angeht. Ihr werdet sicher das Grauen bekommen. Wichtig ist mir momentan nur, dass die grafische Ausgabe wie oben beschrieben in 3D funktioiert. Alles andere ist erstmal egal.

Code: Alles auswählen

from sympy import*
from numpy import*
N=11.0 #Number of space grid points ex)N=11 -> i = 0,1,2,3,4,5,6,7,8,9,10
r=0.25 # (alpha*deltaT)/(deltaX)²
alpha=0.2 
print "------Parameter------"
A=0.5*r
B=1+r
C=A
D=B
dx=1.0/(N-1.0) # space increment x direction
print "Spacestep dx:",dx
dy=dx        # space increment y direction
print "Spacestep dy:",dy
dt=(r*(dx**2))/alpha
print "Timestep dt: ",dt
t=0.5/dt
print "required # of timesteps t: ",t
print "------Spacesteps in x & y direction-----"
i=0
while i<=(N-1):
    step=(i)*dx
    print step
    i=i+1
def f(i,j): # function to derive the tridiagonal matrix form
    if i==j:
        return -B
    elif i==(j+1):
        return A
    elif i==(j-1):
        return A
    else:
        return 0
M=Matrix(N-2,N-2,f)  # construction of a (N-2)x(N-2) tridiagonalmatrix with A(C) on the diagonal and B(D) on both secondary diagonals
print "--------Matrix------"
print M
def q(i,j):
    if i==0:
        return 0
    else:
        return 0
T=Matrix(N-2,N-2,q) 
result=Matrix(N-2,N-2,q)
def u(i,j): # defining the boundary & initial conditions
    u=(100*sin(pi*j*dx))+(100*sin(pi*i*dy))
    return u
def K(i,j):
    K=-u(i,j)-(A*(u(i,j+1)-(2*u(i,j))+u(i,j-1)))
    return K
def g(i,j):
    if i==0:
        return (K(1,j+1)-(A*u(0,j+1)))
    if i==(N-3):
        return (K((N-2),j+1)-(A*(u((N-1),j+1))))
    else:
        return K(i+1,j+1)
K=Matrix(N-2,N-2,g) # Initial right hand side Vektor from Equation (*) - look at my notes
print "Initialvektor K: "
print K
timestep=0
print " Starting the iteration "
while timestep<=1:
    print "------------------------------------------------------ Timestep:",timestep,"--------------------------------------"
    #-------------------------ADI 1.Step----------------------"
    for j in range (0,N-2):
        U=M.LUsolve(K[:,j])
        T[:,j]=U #"intermediate Solution @ timestep n+0.5"
    #-------------------------ADI 2.Step-----------------------"
    def L(i,j):
        L=-T[i,j]-(C*(T[i+1,j]-(2*T[i,j])+T[i-1,j]))
        return L
    def h(i,j):
        if j==0 and i!=0 and i!=N-3:
            return (-T[i,0]-(A*(T[i+1,0]-(2*T[i,0])+T[i-1,0])))-(C*(u(i+1,0)))
        if j==0 and i==0:
            return (-T[0,0]-(A*(T[1,0]-(2*T[0,0])+u(0,1))))-(C*(u(i+1,0)))
        if j==0 and i==N-3:
            return (-T[N-3,0]-(A*(u(0,1)-(2*T[N-3,0])+T[N-4,0])))-(C*(u(i+1,0)))
        if j==(N-3) and i!=0 and i!=N-3:
            return (-T[i,N-3]-(A*(T[i+1,N-3]-(2*T[i,N-3])+T[i-1,N-3])))-(C*(u(i+1,N-1)))
        if j==(N-3) and i==0:
            return (-T[0,N-3]-(A*(T[1,N-3]-(2*T[0,N-3])+u(0,N-2))))-(C*(u(i+1,N-1)))
        if j==(N-3) and i==(N-3):
            return (-T[N-3,N-3]-(A*(u(0,N-2)-(2*T[N-3,N-3])+T[N-4,N-3])))-(C*(u(i+1,N-1)))
        if i==0 and j!=0 and j!=N-3:
            return (-T[0,j]-(A*(T[1,j]-(2*T[0,j])+u(0,j+1))))
        if i==N-3 and j!=0 and j!=N-3:
            return (-T[N-3,j]-(A*(u(N-1,j+1)-(2*T[N-3,j])+T[N-4,j])))
        else:
            return L(i,j)
    L=Matrix(N-2,N-2,h) # the right hand vector of the second step of the ADI scheme - referred to my notations
    for i in range (0,N-2):
        P=M.LUsolve(L[:,i])
        B=P.transpose()
        T[i,:]=B
    print"Solution @ timestep n+1"
    print T
    def L(i,j):
        L=-T[i,j]-(C*(T[i+1,j]-(2*T[i,j])+T[i-1,j]))
        return L
    def h(i,j):
        if j==0 and i!=0 and i!=N-3:
            return (-T[i,0]-(A*(T[i+1,0]-(2*T[i,0])+T[i-1,0])))-(C*(u(i+1,0)))
        if j==0 and i==0:
            return (-T[0,0]-(A*(T[1,0]-(2*T[0,0])+u(0,1))))-(C*(u(i+1,0)))
        if j==0 and i==N-3:
            return (-T[N-3,0]-(A*(u(0,1)-(2*T[N-3,0])+T[N-4,0])))-(C*(u(i+1,0)))
        if j==(N-3) and i!=0 and i!=N-3:
            return (-T[i,N-3]-(A*(T[i+1,N-3]-(2*T[i,N-3])+T[i-1,N-3])))-(C*(u(i+1,N-1)))
        if j==(N-3) and i==0:
            return (-T[0,N-3]-(A*(T[1,N-3]-(2*T[0,N-3])+u(0,N-2))))-(C*(u(i+1,N-1)))
        if j==(N-3) and i==(N-3):
            return (-T[N-3,N-3]-(A*(u(0,N-2)-(2*T[N-3,N-3])+T[N-4,N-3])))-(C*(u(i+1,N-1)))
        if i==0 and j!=0 and j!=N-3:
            return (-T[0,j]-(A*(T[1,j]-(2*T[0,j])+u(0,j+1))))
        if i==N-3 and j!=0 and j!=N-3:
            return (-T[N-3,j]-(A*(u(N-1,j+1)-(2*T[N-3,j])+T[N-4,j])))
        else:
            return L(i,j)
    K=Matrix(N-2,N-2,h)
    #print " The refreshed right hand side vector for the next timestep -> K(new):"
    #print K
    timestep=timestep+1
def d(i,j):
    if i==0 or i==N-1 or j==0 or j==N-1:
        return (100*sin(pi*j*dx))+(100*sin(pi*i*dy))
    else:
        return 0
R=Matrix(N,N,d) # final solution with the temperautre distribution: along the ROWS of Matrix R -> find the temperatures along the X-AXIS.
R[1:N-1,1:N-1]=T # -------------------------------------------------along the COLUMNS of Matrix R -> find the temperature along Y-AXIS.
print R
def f(X,Y):
    f=R[Y,X]
    return f
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import matplotlib.pyplot as plt

fig = plt.figure()
ax = Axes3D(fig)
X = arange(-5, 5, 0.25)
Y = arange(-5, 5, 0.25)
X, Y = meshgrid(X, Y)
Z = f
ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet)

plt.show()[code=]

[/code]

[sliderchrisss]

Zeile 144 ist dann mein Problem!
Er will mir " f " einfach nicht als Variable geben. Wobei "f" ja jeweils ein Eintrag aus meiner Matrix R darstellt.

[xpilz]

Für den ersten Blick ist f für mich eine Funktion. Du hast nur in der Funktion
f, noch ein Array definiert das auch f heißt. Versuche doch mal das f als Funktion aufzurufen und nicht als Variable?.

[mkesper]

Hallo sliderchriss, für Code bitte den "Python"-Button benutzen, Code in dieser Länge kann das Forum allerdings nicht verwenden, deshalb ist es besser, ihn z.B. nach http://paste.pocoo.org auszulagern und hier nur als Link einzubinden. Das Problem mit dem doppelten "f" umgehst du, wenn du versuchst, den Dingen "sprechende" Namen zu geben. Was "macht" die Funktion "f"? So sollte sie auch heißen.