%matplotlib inline
import pylab as pl
import numpy as np
现有一个 m*n 的矩阵,矩阵外围一周的点的值是已知的,矩阵当中每个点的值都是上下左右四个点的平均值,请求出每个点的值。
M, N = 20, 40
np.random.seed(42)
m = np.random.rand(M, N)
m[1:-1, 1:-1] = 0
count = 0
while True:
count += 1
avg = 0.25 * (m[:-2, 1:-1] + m[2:, 1:-1] + m[1:-1, :-2] + m[1:-1, 2:])
err = avg - m[1:-1, 1:-1]
if np.max(np.abs(err)) < 1e-15:
break
m[1:-1, 1:-1] = avg
print count
3490
pl.imshow(m, interpolation="nearest");
from scipy import sparse
from scipy.sparse import linalg
idx = np.arange(M * N).reshape(M, N)
A = sparse.dok_matrix((M*N, M*N))
b = sparse.dok_matrix((M*N, 1))
i = 0
for row in range(M):
for col in range(N):
if row == 0 or col == 0 or row == M-1 or col == N-1:
A[i, idx[row, col]] = 1
b[i, 0] = m[row, col]
else:
A[i, idx[row, col]] = 4
A[i, idx[row-1, col]] = -1
A[i, idx[row+1, col]] = -1
A[i, idx[row, col-1]] = -1
A[i, idx[row, col+1]] = -1
i += 1
x = linalg.spsolve(A, b)
x.shape = M, N
np.allclose(x, m)
True