Matlab数组尺寸不兼容错误排查与修复请求（拉普拉斯有限差分法）

运行基于拉普拉斯有限差分法的Matlab代码时，持续收到“Arrays have incompatible sizes for this operation”错误提示。工作区中所有beta矩阵与h矩阵均显示为201x203维度，理论上尺寸兼容，但怀疑错误源于循环后半段添加不可渗透边界条件前的短暂维度变化。

原代码如下：

% Defining beta values
i=2:ny1-1; 
j=2:nx1-1;

oldbeta1(i,j) = (deltay*(k(i,j)+k(i-1,j))/(2*deltax)); %units: mm^2
beta1 = padarray(oldbeta1,[1,1],0,'post');

oldbeta2(i,j) = (deltay*(k(i+1,j)+k(i,j))/(2*deltax)); %units: mm^2
beta2 = padarray(oldbeta2,[1,1],0,'post');

oldbeta3(i,j) = (deltax*(k(i,j-1)+k(i,j))/(2*deltay)); %units: mm^2
beta3 = padarray(oldbeta3,[1,1],0,'post');

oldbeta4(i,j) = (deltax*(k(i,j)+k(i,j+1))/(2*deltay)); %units: mm^2
beta4 = padarray(oldbeta4,[1,1],0,'post');

% Initial Pressure Conditions (Zero Matrix) and pressure boundary conditions
h = zeros([201 203]); %h is pressure head
h(1,:) = 0; %pressure at exit (atmospheric, mm)
h(ny1,:) = InletPressure; %pressure at entrance (mm)

% Calculating Pressure Head 
for it = 1:niter
    hn = h;
    h(i,j) = ((beta1.*hn(i-1,j))+(beta2.*hn(i+1,j))+(beta3.*hn(i,j-1))+(beta4.*hn(i,j+1))) / (beta1+beta2+beta3+beta4);
    % Impermeable Edge Boundary Conditions
    h(:,1) = h(:,2); 
    h(:,nx1) = h(:,nx1-1); 
    h(1,:) = 0; 
    h(ny1,:) = InletPressure;
end

错误原因分析 #

核心问题出在循环内的压力头计算行：你使用了完整的beta矩阵（201x203）与hn的子矩阵（尺寸为(ny1-2)x(nx1-2)）进行逐元素乘法，两者维度不匹配。虽然工作区显示beta和h整体尺寸一致，但在赋值h(i,j)时，你需要的是beta矩阵中对应内部网格点的子矩阵，而非整个矩阵。

修复后的代码 #

% Defining beta values
i=2:ny1-1; 
j=2:nx1-1;

% 预先初始化oldbeta矩阵，避免自动扩展带来的尺寸偏差
oldbeta1 = zeros(ny1, nx1);
oldbeta1(i,j) = (deltay*(k(i,j)+k(i-1,j))/(2*deltax)); %units: mm^2
beta1 = padarray(oldbeta1,[1,1],0,'post');

oldbeta2 = zeros(ny1, nx1);
oldbeta2(i,j) = (deltay*(k(i+1,j)+k(i,j))/(2*deltax)); %units: mm^2
beta2 = padarray(oldbeta2,[1,1],0,'post');

oldbeta3 = zeros(ny1, nx1);
oldbeta3(i,j) = (deltax*(k(i,j-1)+k(i,j))/(2*deltay)); %units: mm^2
beta3 = padarray(oldbeta3,[1,1],0,'post');

oldbeta4 = zeros(ny1, nx1);
oldbeta4(i,j) = (deltax*(k(i,j)+k(i,j+1))/(2*deltay)); %units: mm^2
beta4 = padarray(oldbeta4,[1,1],0,'post');

% Initial Pressure Conditions (Zero Matrix) and pressure boundary conditions
h = zeros([201 203]); %h is pressure head
h(1,:) = 0; %pressure at exit (atmospheric, mm)
h(ny1,:) = InletPressure; %pressure at entrance (mm)

% Calculating Pressure Head 
for it = 1:niter
    hn = h;
    % 关键修改：使用beta矩阵的对应子矩阵beta1(i,j)，而非整个beta1
    h(i,j) = ((beta1(i,j).*hn(i-1,j))+(beta2(i,j).*hn(i+1,j))+(beta3(i,j).*hn(i,j-1))+(beta4(i,j).*hn(i,j+1))) / (beta1(i,j)+beta2(i,j)+beta3(i,j)+beta4(i,j));
    % Impermeable Edge Boundary Conditions
    h(:,1) = h(:,2); 
    h(:,nx1) = h(:,nx1-1); 
    h(1,:) = 0; 
    h(ny1,:) = InletPressure;
end

额外优化建议 #

• 预先初始化oldbeta矩阵：避免Matlab自动扩展矩阵带来的潜在尺寸偏差，确保oldbeta初始尺寸与ny1、nx1一致。
• 确认ny1和nx1的值：确保ny1=201、nx1=203，与h矩阵的尺寸对应，避免因变量不匹配导致的隐形错误。

内容的提问来源于stack exchange，提问作者Brittney