matlab 使用矩阵作为正方形点阵，并以定义的概率随时间激活每个相邻元素

egdjgwm8 于 2022-11-15 发布在 Matlab

关注(0)|答案(1)|浏览(263)

**注意：我提供的代码是用MatLab编写的，但是，如果您对R-Studio比较熟练，我非常乐意为您提供指针/代码。我只是碰巧稍微更好地掌握了带循环的MATLAB的语法。

我之前有一个问题是关于我在R-STUDIO中正在做的类似事情：Bond Percolation on the square lattice in R-studio
但由于我更加流利地使用MatLab，我想我应该试着从头开始编写我想要的代码。
因此，对于我的代码，我基本上想做这样的事情：

长度为N*N的状态向量，其中0=易感，1=传染性，2=已恢复，其中不能被再次感染；
循环，以便在每个时间步，代码检查列表中第一个受感染节点(状态向量)的相邻邻居(非对角线)，如果有任何受影响的节点，则应用概率p，如果发生传输，则更新状态向量；
检查所有邻居后，将该受感染节点的状态向量从1更新为2；转到下一个受感染节点；
每一代的感染人数就是每个时间步结束时状态向量中的1的数量。

我现在的代码如下所示：

%% Setting inital conditions
;
N=6;                % Set initial matrix size
n=3;                % n and m are initial activation node
m=4
t=10;               % time frame of infection
p = 2/3             % probability of transmission
X = zeros(N)        % Set initial matrix using size N

%% For loop iterating over time t
for i = 1:t
    
    X(n,m) = 1;
    
    if X(n,m) == 1
        % if loops to check if a node is susceptable (equals 0) and if it
        % is assign 0 or 1 based on a predefined probability. (Using
        % binornd(). 
        if n+i > size(X,1)
            break;
        elseif X(n+i,m)==0
            X(n+i,m)= binornd(1,p);
        end
        %
        if m+i > size(X,2);
            break;
        elseif X(n,m+i)==0
            X(n,m+i)= binornd(1,p);
        end
        %
        if n-i < 1
            break;
        elseif X(n-i,m)==0
            X(n-i,m)= binornd(1,p);
        end
        %
        if m-i < 1
            break;
        elseif X(n,m-i)==0
            X(n,m-i)= binornd(1,p);
        end
        
        
    end
    X(n,m) = 2
end

所以，我的代码只检查来自起始节点的十字形的值，我想要做的是在for循环的每个时间步，我想检查矩阵中的所有值，以寻找‘受感染的’节点‘，并在每个受感染的节点上迭代if条件。

即在时间步骤1，我们有1个感染节点，每个节点有2/3的机会感染其周围的4个节点，在此示例中，我们在该步骤结束时感染3个节点，初始感染节点恢复(值2)。

在时间步骤2，三个节点感染其相邻节点，最初的三个感染节点恢复(值2)。
重复此操作，直到For循环结束。

一些额外的好处是将这些节点画在正方形格子图上，甚至开发出一张图，显示随着感染概率的增加，感染人数如何接近1。
谢谢你的阅读。

matlab

来源：https://stackoverflow.com/questions/74022176/using-a-matrix-as-a-square-lattice-and-activating-each-adjacent-element-over-t

1条答案

按热度按时间

8i9zcol21#

要运行单个模拟，这里有一种矢量化方法，它更新尚未感染和当前感染的每个时间步的索引。它只在时间步长上循环。它返回模拟结束时的晶格状态。

function X = SimInfect(N, n, m, t, p)
% set up constants
    N2 = N + 2;
    idx = 2:N + 1;
    offset = [-1, 1, -N2, N2];
    % the probabilities of becoming infected given i - 1 infected
    % neighbors, with i being the index of P
    P = [0, 1 - (1 - p).^(1:4)];
% set up variables
    I = N2*m + n + 1; % the indices of the currently infected
    % Set up the lattice. Previously infected = NaN. Pad with NaN.
    X = NaN(N2, N2);
    X(idx, idx) = 0;
    X(I) = 1;
    U = find(X == 0); % the indices of the not yet infected
% run the simulation
    for i = 1:t
        % for each not-yet-infected, do they become infected this round?
        bln = rand(1, length(U)) < P(sum(X(U + offset), 2, 'omitnan') + 1);
        X(U(bln)) = 1; % update the state of the newly infected
        X(I) = NaN; % update the state of the previously infected
        I = U(bln); % update the indices of the currently infected
        U = U(~bln); % update the indices of the not-yet-infected
    end

    X = X(idx, idx); % remove the padding
    X(isnan(X)) = 2; % re-code previously infected as 2
end

为了生成统计数据，我们可以进一步向量化所需的复制次数：

function nInfected = SimInfectMult(N, n, m, t, p, reps)
    N2 = N + 2;
    P = [0, 1 - (1 - p).^(1:4)];
    NN = N2^2;
    offset = [-1, 1, -N2, N2];
    I = N2*m + n + 1:NN:NN*reps;
    X = NaN(N2, N2, reps);
    X(2:N + 1, 2:N + 1, :) = 0;
    X(I) = 1;
    U = find(X == 0);

    for i = 1:t
        bln = rand(1, length(U)) < P(sum(X(U + offset), 2, 'omitnan') + 1);
        X(U(bln)) = 1;
        X(I) = NaN;
        I = U(bln);
        U = U(~bln);
    end

    nInfected = N^2 - histcounts(U, 1:NN:NN*(reps + 1));
end

10万次模拟只花了我的机器一秒钟多一点的时间：

>> timeit(@() SimInfectMult(6, 3, 4, 10, 2/3, 1e5))

ans =

    1.2663

还请注意，模拟可以用多个受感染的对象进行初始化，例如：

X = SimInfect(6, [1, 3], [6, 4], 10, 2/3);

赞(0）回复(0）举报 2022-11-15

我来回答

matlab 使用矩阵作为正方形点阵，并以定义的概率随时间激活每个相邻元素

1条答案

相关问题

热门标签

最新问答