python - 尝试将 Matlab 转换为 python：（https://en.wikipedia.org/wiki/Laplacian_matrix#Example_of_the_operator_on_a_grid）

Question

我试图了解拉普拉斯矩阵的维基百科页面上的示例代码。

它是用 MatLab 编写的，我只能使用开源工具。

我认为其中大部分内容都相当简单，但我在这一行遇到了一些麻烦。

C0V = V'*C0; % Transform the initial condition into the coordinate system
% of the eigenvectors

现在我的数学还没有达到真正理解评论的含义。从 matlab 网站判断，这似乎是一个转置矩阵乘以另一个矩阵（内积） 。左侧矩阵（转置后）是（m x p），右侧是（p x n）。

该矩阵V是通过调用eig上面的几行生成的，我从这个答案scipy.linalg.eig中替换了.

问题是 C0 明确定义为NxN矩阵（ndarray），但在我的代码V中是N**2xN**2矩阵。我无法知道V原始代码中的形状。

作为参考，下面是维基百科代码，下面是我尝试用 python 重写它（使用 scipy 和 numpy），然后是归因于上述行的错误。

来自维基百科的 Matlab 代码

N = 20; % The number of pixels along a dimension of the image
A = zeros(N, N); % The image
Adj = zeros(N * N, N * N); % The adjacency matrix

% Use 8 neighbors, and fill in the adjacency matrix
dx = [- 1, 0, 1, - 1, 1, - 1, 0, 1];
dy = [- 1, - 1, - 1, 0, 0, 1, 1, 1];
for x = 1:N
    for y = 1:N
        index = (x - 1) * N + y;
        for ne = 1:length(dx)
            newx = x + dx(ne);
            newy = y + dy(ne);
            if newx > 0 && newx <= N && newy > 0 && newy <= N
                index2 = (newx - 1) * N + newy;
                Adj(index, index2) = 1;
            end
        end
    end
end

% BELOW IS THE KEY CODE THAT COMPUTES THE SOLUTION TO THE DIFFERENTIAL EQUATION
Deg = diag(sum(Adj, 2)); % Compute the degree matrix
L = Deg - Adj; % Compute the laplacian matrix in terms of the degree and adjacency matrices
[V, D] = eig(L); % Compute the eigenvalues/vectors of the laplacian matrix
D = diag(D);

% Initial condition (place a few large positive values around and
% make everything else zero)
C0 = zeros(N, N);
C0(2:5, 2:5) = 5;
C0(10:15, 10:15) = 10;
C0(2:5, 8:13) = 7;
C0 = C0(:);

C0V = V'*C0; % Transform the initial condition into the coordinate system
% of the eigenvectors
for t = 0:0.05:5
    % Loop through times and decay each initial component
    Phi = C0V .* exp(- D * t); % Exponential decay for each component
    Phi = V * Phi; % Transform from eigenvector coordinate system to original coordinate system
    Phi = reshape(Phi, N, N);
    % Display the results and write to GIF file
    imagesc(Phi);
    caxis([0, 10]);
     title(sprintf('Diffusion t = %3f', t));
    frame = getframe(1);
    im = frame2im(frame);
    [imind, cm] = rgb2ind(im, 256);
    if t == 0
        imwrite(imind, cm, 'out.gif', 'gif', 'Loopcount', inf, 'DelayTime', 0.1);
    else
        imwrite(imind, cm, 'out.gif', 'gif', 'WriteMode', 'append', 'DelayTime', 0.1);
    end
end

我尝试的翻译

import numpy as np
import matplotlib.pyplot as plt

import scipy.linalg as la


N = 20 # The number of pixels along a dimension of the image
A = np.zeros((N, N)) # The image
Adj = np.zeros((N**2, N**2)) # The adjacency matrix

# Use 8 neighbors, and fill in the adjacency matrix
dx = [- 1, 0, 1, - 1, 1, - 1, 0, 1]
dy = [- 1, - 1, - 1, 0, 0, 1, 1, 1]
for x in range(N):
    for y in range(N):
        index = x * N + y
        for ne in range(len(dx)):
            newx = x + dx[ne]
            newy = y + dy[ne]
            if (newx >= 0 and newx < N 
                and newy >= 0 and newy < N):
                index2 = newx * N + newy;
                Adj[index, index2] = 1


# BELOW IS THE KEY CODE THAT COMPUTES THE SOLUTION TO THE DIFFERENTIAL EQUATION
Deg = np.diag(np.sum(Adj, 1)) # Compute the degree matrix
L = Deg - Adj # Compute the laplacian matrix in terms of the degree and adjacency matrices
D, V = la.eig(L) # Compute the eigenvalues/vectors of the laplacian matrix
D = np.diag(D)


# Initial condition (place a few large positive values around and
# make everything else zero)
C0 = np.zeros((N, N))
C0[1:4, 1:4] = 5
C0[9:14, 9:14] = 10
C0[1:5, 7:12] = 7

#C0 = C0(:) #This doesn't seem to do anything?


# matlab:C0V = V'*C0; % Transform the initial condition into the coordinate system
# of the eigenvectors

C0V = V.T * C0 # ???

错误

----------------------------------------------------------------------
ValueError                           Traceback (most recent call last)
<ipython-input-10-c3014dc90c9f> in <module>
      2 # of the eigenvectors
      3 
----> 4 C0V = V.T * C0

ValueError: operands could not be broadcast together with shapes (400,400) (20,20) 

编辑

看来@hpaulj 已经确定了我的问题。我省略的行是一个ravel操作C0 = C0(:)，即将20, 20矩阵转换为400, 1向量。所以我需要

C0 = C0.ravel()

标签： pythonmatlabnumpy