如何对Scilab中的所有矩阵元素执行操作？

Question

我试图模拟随着时间的推移无限板上的热量分布。为此，我写了一个Scilab脚本。现在，它的关键点是计算所有板点的温度，并且必须在我想要观察的每个时刻进行：

for j=2:S-1
    for i=2:S-1
        heat(i, j) = tcoeff*10000*(plate(i-1,j) + plate(i+1,j) - 4*plate(i,j) + plate(i, j-1) + plate(i, j+1)) + plate(i,j);          
    end;
end 

问题是，如果我想为一个100x100点的盘子做，那就意味着，这里（它只适用于内部，没有边界条件），我将不得不循环98x98 = 9604次，每一次计算给定i,j点的热量。如果我想观察那个，比如100个secons，1步，我必须重复100次，共计960,400次迭代。这花了很长时间，我想避免它。高达50x50的盘子，这一切都发生在合理的4-5秒时间范围内。

现在我的问题是 - 是否有必要使用for循环完成所有这些？在Scilab中是否有任何内置的聚合函数，这将让我对矩阵的所有元素执行此操作？我还没有找到方法的原因是，每个点的结果都取决于其他矩阵点的值，这使得我使用嵌套循环。关于如何让它更快被欣赏的任何想法。

Answer 1

在我看来，你想要计算你的热场和某种扩散模式的二维相互关系。这种模式可以被认为是“过滤器”内核，这是使用线性过滤器矩阵修改图像的常用方法。你的“过滤器”是：

F=[0,1,0;1,-4,1;0,1,0];

如果您安装图像处理工具箱（IPD），您将有一个MaskFilter函数来执行此2D相互关系。

S=500;
plate=rand(S,S);
tcoeff=1;

//your solution with nested for loops
t0=getdate();
for j=2:S-1
  for i=2:S-1
    heat(i, j) = tcoeff*10000*(plate(i-1,j)+plate(i+1,j)-..
    4*plate(i,j)+plate(i,j-1)+plate(i, j+1))+plate(i,j);          
  end
end
t1=getdate();
T0=etime(t1,t0);
mprintf("\nNested for loops: %f s (100 %%)",T0);

//optimised nested for loop
F=[0,1,0;1,-4,1;0,1,0];   //"filter" matrix
F=tcoeff*10000*F;
heat2=zeros(plate);
t0=getdate();
for j=2:S-1
  for i=2:S-1
    heat2(i,j)=sum(F.*plate(i-1:i+1,j-1:j+1));
  end
end
heat2=heat2+plate;
t1=getdate();
T2=etime(t1,t0);
mprintf("\nNested for loops optimised: %f s (%.2f %%)",T2,T2/T0*100);

//MaskFilter from IPD toolbox
t0=getdate();
heat3=MaskFilter(plate,F);
heat3=heat3+plate;
t1=getdate();
T3=etime(t1,t0);
mprintf("\nWith MaskFilter: %f s (%.2f %%)",T3,T3/T0*100);

disp(heat3(1:10,1:10)-heat(1:10,1:10),"Difference of the results (heat3-heat):");

请注意，MaskFilter在应用滤镜之前填充图像（原始矩阵），据我所知，它使用跨越边界的“镜像”阵列。您应该检查此行为是否适合您。

速度增加约为* 320（执行时间为原始代码的0.32％）。这够快吗？

从理论上讲，它可以通过两个二维傅立叶变换（可能有Scilab内置mfft）完成，但它可能不会比这更快。见这里：http://mailinglists.scilab.org/Image-processing-filter-td2618144.html#a2618168

Answer 2

请注意，矢量化操作和并行计算之间存在很大差异，正如我已经解释过here。尽管矢量化可能会稍微改善性能，但这与您通过GPU计算实现的效果无法比较（例如OpenCL）。我将尝试解释代码的矢量化形式，而不必过多详细介绍。考虑这些给出：

S = ...;
tcoeff = ...;

function Plate = plate(i, j)
    ...;
endfunction

function Heat = heat(i, j)
    ...;
endfunction

现在你可以定义一个meshgrid：

x = 2 : S - 1;
y = 2 : S - 1;
[M, N] = meshgrid(x,y);
Result = feval(M, N, heat);

feval是关键，它将在feval和M矩阵上播放N函数。

Answer 3

您的方案是应用于矩形网格的拉普拉斯算子的有限差分格式。如果您选择自由度的行或列式编号（这里是板（i，j））以便将它们视为向量，那么应用“离散”拉普拉斯算子可以通过乘以稀疏矩阵来完成在左边（非常快）这在以下文件中特别好解释：

https://www.math.uci.edu/~chenlong/226/FDMcode.pdf。

该实现在Matlab中描述，但很容易在Scilab中翻译。