python中的矩阵相关计算

Question

我在计算python中的矩阵更新时发现了一个非常有趣的问题。我必须计算错误（这是前n个更新矩阵之间的差异）。

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，但是当我打印我的矩阵时，我发现它不应该是零。伙计们，我认为我在这里有一些概念上的问题，所以请帮助。

Answer 1

您的问题是，在更新功能中分配C=A; D=A时，您使用浅拷贝，即仅拷贝参考。实质上，在D的构建之后，所有三个变量A,C,D指向相同的对象。使用

def update(A): 
    C=1.0*A 
    D=1.0*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 update(A): 
    D=A.copy() 
    D[1:-1,1:-1]=(A[0:-2,1:-1]+A[2:,1:-1]+A[1:-1,0:-2]+A[1:-1,2:])/4 

    return(np.abs(D-A),D) 

在深副本自动传递参数和进行算术运算的结果。

---

你知道（几何，一阶）收敛速度是一样的东西max(1-C/(nx^2), 1-C/(ny^2))，即速度很慢甚至中等大电网？对于实际应用，最好使用共轭梯度，其他与Krylov相关的算法或多网格方法（或稀疏求解器库，UMFpack ...）。

---

在（未使用）error过程中，应该有不一样的东西

e = abs(A-B) 

目前，我们返回新生成的零矩阵C的常态。