0
我在計算python中的矩陣更新時發現了一個非常有趣的問題。我必須計算錯誤(這是前n個更新矩陣之間的差異)。python中的矩陣相關計算
import numpy as np
import matplotlib.pyplot as plt
#from matplotlib import animation
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
def update(A):
C=A
D=A
D[1:-1,1:-1]=(C[0:-2,1:-1]+C[2:,1:-1]+C[1:-1,0:-2]+C[1:-1,2:])/4
return(np.abs(D-C),D)
def error(A,B):
C=np.zeros(np.shape(A),np.float64)
#e=np.max(np.max(np.abs(C)))
e=(np.abs(C))
return (e.sum(dtype='float64'))
def initial(C):
C[0,:]=0 ## Top Boundary
C[-1,:]=0 ## Bottom Boundary
C[:,0]=0 ## left Boundary
C[:,-1]=100 ## Right Boundary
return(C)
def SolveLaplace(nx, ny,epsilon,imax):
## Initialize the mesh with some values
U = np.zeros((nx, ny),np.float64)
## Set boundary conditions for the problem
U=initial(U)
## Store previous grid values to check against error tolerance
UN=np.zeros((nx, ny),np.float64)
UN=initial(UN)
## Constants
k = 1 ## Iteration counter
## Iterative procedure
while k<imax:
err,U=update(U)
print(err.sum())
k+=1
return (U)
nx = 50.0
ny = 50.0
dx = 0.001
epsilon = 1e-6 ## Absolute Error tolerance
imax = 5000 ## Maximum number of iterations allowed
Z = SolveLaplace(nx, ny,epsilon,imax)
#x = np.linspace(0, nx * dx, nx)
#y = np.linspace(0, ny * dx, ny)
#X, Y = np.meshgrid(x,y)
##===================================================================
def PlotSolution(nx,ny,dx,T):
## Set up x and y vectors for meshgrid
x = np.linspace(0, nx * dx, nx)
y = np.linspace(0, ny * dx, ny)
fig = plt.figure()
ax = fig.gca(projection='3d')
X, Y = np.meshgrid(x,y)
ax.plot_surface(X, Y, T.transpose(), rstride=1, cstride=1, cmap=cm.cool, linewidth=0, antialiased=False)
plt.xlabel("X")
plt.ylabel("Y")
#plt.zlabel("T(X,Y)")
plt.figure()
plt.contourf(X, Y, T.transpose(), 32, rstride=1, cstride=1, cmap=cm.cool)
plt.colorbar()
plt.xlabel("X")
plt.ylabel("Y")
plt.show()
##===================================================================
PlotSolution(nx, ny, dx, Z)
我想解決拉普拉斯方程2-d片(溫度分佈),並且當誤差小於某一最小值,平衡將得以實現。但是在計算錯誤時,我總是得到0,但是當我打印我的矩陣時,我發現它不應該是零。夥計們,我認爲我在這裏有一些概念上的問題,所以請幫助。
我覺得'A.copy()'更清晰一點 – RodericDay
是的。我依靠http://stackoverflow.com/a/6435446/3088138。但是避免算術運算應該快一些。 – LutzL
謝謝@LutzL ..它現在真的很好用.. – Gaurav16