粒子-体素模拟中循环逻辑异常问题求助

我正在运行一个模拟，统计由M×M×M体素构成的L×L×L立方体中每个体素内的粒子数量，粒子总数为N。

原代码 #

function CountParticles(L,M,N)
    coordlist = {};
for counter = 1:N
    coordlist{counter} = L*rand(1,3);
end
step = L/M;
positlist={};
posx = 1;
posy = 1;
posz = 1;
mult = 1;
for counter = 1:N
    cellarray = coordlist{counter};
    checkx = cellarray(1,1);
    checky = cellarray(1,2);
    checkz = cellarray(1,3);
    if checkx > step*mult
        mult = mult+1;
    else 
        posx = mult;
        mult = 1;
    end
    if checky > step*mult
        mult = mult+1;
    else
        posy = mult;
        mult = 1;
    end
    if checkz > step*mult
        mult = mult+1;
    else
        posz = mult;
        mult = 1;
    end
    checkempty = isempty(positlist{0,1,1});
    disp(checkempty)
    if checkempty == 1
        positlist{posx,posy,posz} = 1;
    else
        positlist{posx,posy,posz} = positlist{posx,posy,posz}+1;
    end
  end

初始问题 #

原代码逻辑是：生成N个粒子的3维坐标，计算体素步长step=L/M，通过比较粒子坐标与step*mult确定所属体素索引，再统计每个体素的粒子数。但遇到两个问题：

• 体素索引posx/posy/posz未按预期更新，数值随机混乱
• 报错：Index in position 1 exceeds array bounds. Error in CountParticles (line 42) checkempty = isempty(positlist{1,1,1});

修改后代码（已解决报错，但索引仍异常） #

function CountParticles(L,M,N)
coordlist = cell(1,N);
for counter = 1:N
coordlist{counter} = L*rand(1,3);
end
step = L/M;
positlist=cell(M,M,M);
posx = 1;
posy = 1;
posz = 1;
mult = 1;
for counter = 1:N
cellarray = coordlist{counter};
checkx = cellarray(1,1);
checky = cellarray(1,2);
checkz = cellarray(1,3);
if checkx > step*mult
    mult = mult+1;
else 
    posx = mult;
    mult = 1;
end
if checky > step*mult
    mult = mult+1;
else
    posy = mult;
    mult = 1;
end
if checkz > step*mult
    mult = mult+1;
else
    posz = mult;
    mult = 1;
end
checkempty = isempty(positlist{posx,posy,posz});
if checkempty == 1
    positlist{posx,posy,posz} = 1;
else
    positlist{posx,posy,posz} = positlist{posx,posy,posz}+1;
end
end
disp(positlist)
end

索引异常的原因 #

核心问题在于体素索引计算逻辑错误：

• 原代码用单个if判断坐标与step*mult的大小，仅能让mult自增一次，无法找到正确的体素索引。例如粒子x坐标为3*step时，代码只会判断3*step>step*1，将mult设为2后就结束判断，最终posx=2，但实际应属于第3个体素。
• 处理x/y/z坐标时复用了mult变量，虽然每次处理后重置为1，但单个if的逻辑本身就无法遍历到正确的索引值。

解决建议 #

替换错误的索引计算逻辑，直接通过坐标与步长的关系推导体素索引：

方法1：使用ceil函数 #

Matlab中ceil(x)返回大于等于x的最小整数，结合步长计算：

posx = ceil(checkx / step);
posy = ceil(checky / step);
posz = ceil(checkz / step);

注：若坐标恰好等于k*step（k为整数），ceil会返回k，对应第k个体素，符合逻辑。

方法2：使用floor函数（避免边界浮点问题） #

若担心浮点精度导致的边界误判，可添加eps修正：

posx = floor((checkx - eps) / step) + 1;
posy = floor((checky - eps) / step) + 1;
posz = floor((checkz - eps) / step) + 1;

优化后的完整代码 #

function CountParticles(L,M,N)
    % 生成粒子坐标
    coordlist = cell(1,N);
    for counter = 1:N
        coordlist{counter} = L*rand(1,3);
    end
    
    step = L/M;
    % 初始化M×M×M的计数数组，默认值设为0更高效
    positlist = zeros(M,M,M);
    
    for counter = 1:N
        checkx = coordlist{counter}(1,1);
        checky = coordlist{counter}(1,2);
        checkz = coordlist{counter}(1,3);
        
        % 计算正确的体素索引
        posx = ceil(checkx / step);
        posy = ceil(checky / step);
        posz = ceil(checkz / step);
        
        % 直接累加计数，无需判断空值
        positlist(posx, posy, posz) = positlist(posx, posy, posz) + 1;
    end
    
    disp(positlist)
end

注：将positlist改为数值数组zeros(M,M,M)，比cell数组更高效，且无需判断空值，直接累加即可。

内容的提问来源于stack exchange，提问作者Nicolò Salvadori