下面复制了演示。G中的编号是不同的，但数字只是标签(标签网格让我困惑)。在import numpy as np

from scipy import sparse

from scipy.sparse import linalg

import matplotlib.pyplot as plt

def numgrid(n):

"""

NUMGRID Number the grid points in a two dimensional region.

G = NUMGRID('R',n) numbers the points on an n-by-n grid in

an L-shaped domain made from 3/4 of the entire square.

adapted from C. Moler, 7-16-91, 12-22-93.

Copyright (c) 1984-94 by The MathWorks, Inc.

"""

x = np.ones((n,1))*np.linspace(-1,1,n)

y = np.flipud(x.T)

G = (x > -1) & (x < 1) & (y > -1) & (y < 1) & ( (x > 0) | (y > 0))

G = np.where(G,1,0) # boolean to integer

k = np.where(G)

G[k] = 1+np.arange(len(k[0]))

return G

def delsq(G):

"""

DELSQ Construct five-point finite difference Laplacian.

delsq(G) is the sparse form of the two-dimensional,

5-point discrete negative Laplacian on the grid G.

adapted from C. Moler, 7-16-91.

Copyright (c) 1984-94 by The MathWorks, Inc.

"""

[m,n] = G.shape

# Indices of interior points

G1 = G.flatten()

p = np.where(G1)[0]

N = len(p)

# Connect interior points to themselves with 4's.

i = G1[p]-1

j = G1[p]-1

s = 4*np.ones(p.shape)

# for k = north, east, south, west

for k in [-1, m, 1, -m]:

# Possible neighbors in k-th direction

Q = G1[p+k]

# Index of points with interior neighbors

q = np.where(Q)[0]

# Connect interior points to neighbors with -1's.

i = np.concatenate([i, G1[p[q]]-1])

j = np.concatenate([j,Q[q]-1])

s = np.concatenate([s,-np.ones(q.shape)])

# sparse matrix with 5 diagonals

return sparse.csr_matrix((s, (i,j)),(N,N))

if __name__ == '__main__':

print numgrid(6)

print delsq(numgrid(6)).todense()

G = numgrid(32)

D = delsq(G)

N = D.shape[0]

rhs = np.ones((N,1))

u = linalg.spsolve(D, rhs) # vector solution

U = np.zeros(G.shape) # map u back onto G space

U[G>0] = u[G[G>0]-1]

plt.contour(U)

plt.show()

结果：