Hello. I got one question, when i want to streamplot this, i got eror: "raise ValueError("The columns of 'y' must be equal")", i can't figure it out. :/
Thank you for your time.
this is entire code:
import matplotlib.pyplot as plt
import numpy as np
# discretization of domain
x_start, x_end = 0. , 10.
y_start , y_end = 0. , 8.
i = np.linspace(x_start, x_end, 80)
j = np.linspace(y_start, y_end, 50)
X,Y = np.meshgrid(i,j)
Y[:,:40] += 1
Dx = x_end/(np.size(i)-1)
Dy = y_end/(np.size(j)-1)
beta = Dx/Dy
nx = np.size(i) #number of points on x axis
ny = np.size(j) #number of points on y axis
#matrix
nk = nx*ny
A = np.zeros([nk, nk])
B = np.zeros([nk,1])
# solving
for i in range (nx):
for j in range (ny):
k = j*nx + i
if i == 0: # left edge
k1 = j*nx + 1
k0 = j*nx + 0
A[k,k1] = 1
A[k,k0] = -1
B[k] = Dx
elif i == nx-1: # right edge
k_1 = j*nx + (nx-1)
k_2 = j*nx + (nx-2)
A[k,k_2] = -1
A[k,k_1] = 1
B[k] = Dx
elif j == 0: # bottom
K1 = nx + i
K0 = i
A[k,K0] = -1
A[k,K1] = 1
B[k] = 0
elif j == ny-1: # topper
K_1 = (ny-1)*nx + i
K0 = (ny-2)*nx + i
A[k,K_1] = 1
A[k,K0] = -1
B[k] = 0
else: # middle
k_top = (j+1)*nx + i # fi i,j+1
k_down = (j-1)*nx + i # fi i,j-1
k_left = j*nx + (i-1) # fi i-1,j
k_right = j*nx + (i+1) # fi i+1,j
A[k,k_top] = beta**2
A[k,k_down] = beta**2
A[k,k_left] = 1
A[k,k_right] = 1
A[k,k] = -2*(1+beta**2)
B[k] = 0
sol_vec = np.linalg.solve(A,B)
fi = np.reshape(sol_vec,(ny,nx))
# u and v
u = np.zeros([ny,nx])
v = np.zeros([ny,nx])
# velocity u i v
for i in range(nx):
for j in range (ny):
if i == 0: # left edge
v[j,i] = 0
u[j,i] = (fi[j,i+1] - fi[j,i])/Dx
elif i == nx-1: # right edge
v[j,i] = 0
u[j,i] = (fi[j,i] - fi[j,i-1])/Dx
elif j == 0: # bottom, 'v'
v[j,i] = 0
elif j == ny-1: # topper, 'v'
v[j,i] = 0
else: # middle
u[j,i] = (fi[j,i+1] - fi[j,i-1])/(2*Dx)
v[j,i] = (fi[j+1,i] - fi[j-1,i])/(2*Dy)
for i in range (1,nx-1):
u[0,i] = (fi[j,i+1] - fi[j,i-1])/(2*Dx) # 'u'
u[-1,i] = (fi[j,i+1] - fi[j,i-1])/(2*Dx)
u[15,40] = 0
v[15,40] = 0
plt.figure()
plt.quiver(X,Y,u,v)
plt.streamplot(X,Y,u,v,linewidth=1.5,arrowsize=1)
plt.show()