Matlab中integral2运行时vertcat维度不匹配错误求助

问题背景 #

编写二维有限元（FEM）程序，逻辑为：先将计算域划分为N1*N2单元，构建单元与节点信息矩阵，再通过嵌套循环组装刚度矩阵。运行时出现维度不匹配错误，测试单个点fun(0,1)时函数运行正常，但integral2调用失败，需排查解决。

程序代码 #

主程序：FEM_2D_TRI_QUA_1.m #

N1=2; 
N2=2; 
N=2*N1*N2;
top=1;
bottom=0;
left=0;
right=1;
h1=(right-left)/N1;
h2=(top-bottom)/N2;
T=zeros(3,2*N1*N2);
N10=N1+1;
N20=N2+1;

for i=1:N2
    for j=1:2*N1
        if mod(j,2)==1
            T(:,2*N1*(i-1)+j)=[(i-1)*N10+ceil(j/2);i*N10+ceil(j/2);(i-1)*N10+ceil(j/2)+1];
        else
            T(:,2*N1*(i-1)+j)=[(i-1)*N10+j/2+1;i*N10+j/2;i*N10+j/2+1];
        end
    end
end
P=zeros(2,(N1+1)*(N2+1));
for i=1:N1+1
    for j=1:N2+1
        P(:,(i-1)*(N1+1)+j)=[left+(i-1)*h1;bottom+(j-1)*h2];
    end
end
N1_1= 2*N1;
N2_1= 2*N1;
N_1= 2*N1_1*N2_1;

h1_1=(right-left)/N1_1;
h2_1=(top-bottom)/N2_1;
N10_1=N1_1+1;
N20_1=N2_1+1;
P_b=zeros(2,(N1_1+1)*(N2_1+1));
for i=1:N1_1+1
    for j=1:N2_1+1
        P_b(:,(i-1)*(N1_1+1)+j)=[left+(i-1)*h1_1;bottom+(j-1)*h2_1];
    end
end
boundarynodes=zeros(2,N2_1+1+N1_1+N2_1+N1_1-1);
for j=1:2*N2_1+2*N1_1
    boundarynodes(1,j)=-1;
end

for i=1:N1_1+1
    boundarynodes(2,i)=(i-1)*(N1_1+1)+1;
end
for i=N1_1+2:N1_1+1+N2_1
    boundarynodes(2,i)=(N1_1+1-1)*(N1_1+1)+i-(N1_1);
end
for i=N1_1+1+N2_1+N1_1:-1:N1_1+N2_1+2
    boundarynodes(2,i)=((-i+2*N1_1+N2_1+2)-1)*(N1_1+1)+N1_1+1;
end
for i=N1_1+N2_1+N2_1+N1_1:-1:N1_1+2+N2_1+N1_1
    boundarynodes(2,i)=-i+2*N1_1+2*N2_1+2;
end

T_b=zeros(6,2*N1*N2);
c2i=@(i,j) (i-1)*(N1_1+1)+j;
for i=1:N1
    for j=1:2*N2
        if mod(j,2)==1
            i_0=2*i-1;
            j_0=2*ceil(j/2)-1;
            T_b(:,2*N1*(i-1)+j)=[c2i(i_0,j_0);c2i(i_0+2,j_0);c2i(i_0,j_0+2);c2i(i_0+1,j_0);c2i(i_0+1,j_0+1);c2i(i_0,j_0+1)];
        else
            i_0=2*i-1;
            j_0=2*(j/2+1)-1;
            T_b(:,2*N1*(i-1)+j)=[c2i(i_0,j_0);c2i(i_0+2,j_0-2);c2i(i_0+2,j_0);c2i(i_0+1,j_0-1);c2i(i_0+2,j_0-1);c2i(i_0+1,j_0)];
        end
    end
end
A=sparse((N1_1+1)*(N2_1+1),(N1_1+1)*(N2_1+1));
for n=1:N
    y_n2=P_b(2,T_b(2,n));
    y_n3=P_b(2,T_b(3,n));
    y_n1=P_b(2,T_b(1,n));
    y_n4=P_b(2,T_b(4,n));
    y_n5=P_b(2,T_b(5,n));
    y_n6=P_b(2,T_b(6,n));
    x_n3=P_b(1,T_b(3,n));
    x_n2=P_b(1,T_b(2,n));
    x_n1=P_b(1,T_b(1,n));
    x_n4=P_b(1,T_b(4,n));
    x_n5=P_b(1,T_b(5,n));
    x_n6=P_b(1,T_b(6,n));
    Y_n2=P(2,T(2,n));
    Y_n3=P(2,T(3,n));
    Y_n1=P(2,T(1,n));
    X_n3=P(1,T(3,n));
    X_n2=P(1,T(2,n));
    X_n1=P(1,T(1,n));
    xmin=X_n1;
    xmax=xmin+h1;
    J_n=(X_n2-X_n1)*(Y_n3-Y_n1)-(X_n3-X_n1)*(Y_n2-Y_n1);
    x_hat=@(x,y) ((Y_n3-Y_n1)*(x-X_n1)-(X_n3-X_n1)*(y-Y_n1))/J_n;
    y_hat=@(x,y) ((X_n2-X_n1)*(y-Y_n1)-(Y_n2-Y_n1)*(x-X_n1))/J_n;
    if mod(n,2)==1 
       ymin= Y_n1;
       ymax=@(x) Y_n3+((Y_n3-Y_n2)/(X_n3-X_n2))*(x-X_n3);
    else
       ymin=@(x) Y_n2+((Y_n2-Y_n1)/(X_n2-X_n1))*(x-x_n2);
       ymax= Y_n1;
    end
    for i=1:6
        for j=1:6
            fun=@(x,y)...
                    (p_i_x(x_hat(x,y),y_hat(x,y),i))*((Y_n3-Y_n1)/J_n)+...
                    (p_i_y(x_hat(x,y),y_hat(x,y),i))*((Y_n1-Y_n2)/J_n)*...
                    (p_i_x(x_hat(x,y),y_hat(x,y),j)*((Y_n3-Y_n1)/J_n)+...
                    (p_i_y(x_hat(x,y),y_hat(x,y),j))*((Y_n1-Y_n2)/J_n))+...
                    ((p_i_x(x_hat(x,y),y_hat(x,y),i))*((X_n1-X_n3)/J_n)+...
                    (p_i_y(x_hat(x,y),y_hat(x,y),i))*((X_n2-X_n1)/J_n))*...
                    ((p_i_x(x_hat(x,y),y_hat(x,y),j))*((X_n1-X_n3)/J_n)+...
                    (p_i_y(x_hat(x,y),y_hat(x,y),j))*(X_n2-X_n1)/J_n);
            r=integral2(fun,xmin,xmax,ymin,ymax);
            A(T_b(j,n),T_b(i,n))=A(T_b(j,n),T_b(i,n))+r;
        end
    end
end
A

依赖函数 #

p_i_x.m

function r=p_i_x(x,y,i)
    i_x=[4*x+4*y-3;4*x-1;0;-8*x-4*y+4;4*y;-4*y];
    r=i_x(i);
end

p_i_y.m

function r=p_i_y(x,y,i)
    i_y=[4*x+4*y-3;0;4*y-1;-4*x;4*x;-8*y-4*x+4];
    r=i_y(i);
end

错误信息 #

FEM_2D_TRI_QUA_1
Error using vertcat
Dimensions of arrays being concatenated are not consistent.

Error in p_i_x (line 3)
    i_x=[4*x+4*y-3;4*x-1;0;-8*x-4*y+4;4*y;-4*x];

Error in
FEM_2D_TRI_QUA_1>@(x,y)(p_i_x(x_hat(x,y),y_hat(x,y),i))*((Y_n3-Y_n1)/J_n)+(p_i_y(x_hat(x,y),y_hat(x,y),i))*((Y_n1-Y_n2)/J_n)*(p_i_x(x_hat(x,y),y_hat(x,y),j)*((Y_n3-Y_n1)/J_n)+(p_i_y(x_hat(x,y),y_hat(x,y),j))*((Y_n1-Y_n2)/J_n))+((p_i_x(x_hat(x,y),y_hat(x,y),i))*((X_n1-X_n3)/J_n)+(p_i_y(x_hat(x,y),y_hat(x,y),i))*((X_n2-X_n1)/J_n))*((p_i_x(x_hat(x,y),y_hat(x,y),j))*((X_n1-X_n3)/J_n)+(p_i_y(x_hat(x,y),y_hat(x,y),j))*(X_n2-X_n1)/J_n)

Error in integral2Calc>integral2t/tensor (line 228)
        Z = FUN(X,Y);  NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
    [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 106)
    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Error in FEM_2D_TRI_QUA_1 (line 120)
            r=integral2(fun,xmin,xmax,ymin,ymax);

排查与解决方案 #

错误原因 #

integral2会向传入的函数传递数组型的x和y参数（而非单个标量），但p_i_x和p_i_y中构造数组时，使用了标量0，导致垂直拼接时维度不匹配：比如4*x+4*y-3是M×N数组，而0是标量，vertcat要求所有行的维度完全一致，因此触发报错。测试单个点时x和y是标量，0和其他标量可以正常拼接，所以不会报错。

修复代码 #

将两个依赖函数中的标量0替换为与输入x/y同维度的零数组，确保所有行维度匹配：

修改后p_i_x.m

function r=p_i_x(x,y,i)
    i_x=[4*x+4*y-3;4*x-1;zeros(size(x));-8*x-4*y+4;4*y;-4*y];
    r=i_x(i);
end

修改后p_i_y.m

function r=p_i_y(x,y,i)
    i_y=[4*x+4*y-3;zeros(size(x));4*y-1;-4*x;4*x;-8*y-4*x+4];
    r=i_y(i);
end

验证逻辑 #

修改后，当x和y为任意维度数组时，i_x和i_y的每一行都与输入维度一致，垂直拼接不会报错，integral2可以正常传递数组参数进行逐元素运算，满足函数要求。

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