Error using horzcat. Dimensions of matrices being concatenated are not consistent while using integral2

3 views (last 30 days)
岳
岳 on 21 Oct 2022
Edited: Torsten on 21 Oct 2022
I want to code a code to calculate the numerical solution of PDE,here is my code
N1=2; %横向分解成2个元素
N2=2; %纵向分解成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);%生成初始T矩阵
N10=N1+1;
N20=N2+1;
%先把每个element的行列转化为每个element对应的起始node的坐标，然和把坐标转化为index
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
%建立坐标到index函数
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
This code depend on two function
function 1
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
function2
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
then error in the title occurs when I execute r=integral2(fun,xmin,xmax,ymin,ymax) in the main programme

Sign in to answer this question.

Answers (1)

Torsten
Torsten on 21 Oct 2022
Edited: Torsten on 21 Oct 2022
Ran in:
fun=@(x,y) ...
The x and y for which you must evaluate your function usually are matrices of the same size, not simple scalar values.
Thus something like
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
does not work in this case.
Define your own function "fun" and evaluate "fun" for the matrices x and y elementwise.
This code should work:
N1=2; %横向分解成2个元素
N2=2; %纵向分解成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);%生成初始T矩阵
N10=N1+1;
N20=N2+1;
%先把每个element的行列转化为每个element对应的起始node的坐标，然和把坐标转化为index
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
%建立坐标到index函数
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)...
arrayfun(@(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),X,Y);
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
A =
(1,1) 0.4167 (2,1) -0.7500 (3,1) 0.0833 (6,1) -0.0833 (7,1) -0.0833 (11,1) -0.0833 (1,2) -1.0000 (2,2) 1.0000 (3,2) -1.0000 (6,2) -0.3333 (7,2) -0.3333 (11,2) -0.3333 (1,3) 0.1667 (2,3) -0.6667 (3,3) 1.4167 (4,3) -0.7500 (5,3) 0.0833 (6,3) -0.0000 (7,3) 0.0000 (8,3) -0.7500 (9,3) -0.0833 (12,3) -0.0000 (13,3) 0.0833 (3,4) -1.0000 (4,4) 1.0000 (5,4) -1.0000 (8,4) -0.3333 (9,4) -0.3333 (13,4) -0.3333 (3,5) 0.1667 (4,5) -0.6667 (5,5) 1.0000 (8,5) -0.0000 (9,5) 0.0000 (10,5) -0.6667 (14,5) -0.0000 (15,5) 0.1667 (1,6) -0.0000 (2,6) -0.0000 (3,6) 0.0000 (6,6) 1.3333 (7,6) -1.3333 (11,6) -0.0000 (1,7) 0.3333 (2,7) 0.3333 (3,7) 0.3333 (6,7) -1.0000 (7,7) 4.3333 (8,7) -1.3333 (11,7) 0.3333 (12,7) -1.3333 (13,7) -0.0000 (3,8) -0.6667 (4,8) -0.0000 (5,8) 0.0000 (7,8) -1.3333 (8,8) 4.0000 (9,8) -1.3333 (11,8) -0.0000 (12,8) -0.0000 (13,8) -0.6667 (3,9) 0.3333 (4,9) 0.3333 (5,9) 0.3333 (8,9) -1.0000 (9,9) 4.3333 (10,9) -1.3333 (13,9) 0.3333 (14,9) -1.3333 (15,9) -0.0000 (5,10) -0.6667 (9,10) -1.3333 (10,10) 2.6667 (13,10) -0.0000 (14,10) -0.0000 (15,10) -0.6667 (1,11) 0.0833 (2,11) 0.0833 (3,11) 0.0833 (6,11) 0.0833 (7,11) 0.0833 (8,11) -0.0000 (11,11) 1.0000 (12,11) -1.4167 (13,11) 0.2500 (16,11) -0.0833 (17,11) -0.0833 (21,11) -0.0833 (3,12) -0.0000 (7,12) -1.3333 (8,12) -0.0000 (11,12) -1.6667 (12,12) 3.6667 (13,12) -1.6667 (16,12) -0.3333 (17,12) -0.3333 (21,12) -0.3333 (3,13) 0.2500 (4,13) 0.0833 (5,13) 0.0833 (7,13) -0.0000 (8,13) -0.5833 (9,13) 0.0833 (10,13) -0.0000 (11,13) 0.3333 (12,13) -1.3333 (13,13) 3.0000 (14,13) -1.4167 (15,13) 0.2500 (16,13) -0.0000 (17,13) 0.0000 (18,13) -0.7500 (19,13) -0.0833 (22,13) -0.0000 (23,13) 0.0833 (5,14) -0.0000 (9,14) -1.3333 (10,14) -0.0000 (13,14) -1.6667 (14,14) 3.6667 (15,14) -1.6667 (18,14) -0.3333 (19,14) -0.3333 (23,14) -0.3333 (5,15) 0.1667 (9,15) -0.0000 (10,15) -0.6667 (13,15) 0.3333 (14,15) -1.3333 (15,15) 2.0000 (18,15) -0.0000 (19,15) 0.0000 (20,15) -0.6667 (24,15) -0.0000 (25,15) 0.1667 (11,16) -0.0000 (12,16) -0.0000 (13,16) 0.0000 (16,16) 1.3333 (17,16) -1.3333 (21,16) -0.0000 (11,17) 0.3333 (12,17) 0.3333 (13,17) 0.3333 (16,17) -1.0000 (17,17) 4.3333 (18,17) -1.3333 (21,17) 0.3333 (22,17) -1.3333 (23,17) -0.0000 (13,18) -0.6667 (14,18) -0.0000 (15,18) 0.0000 (17,18) -1.3333 (18,18) 4.0000 (19,18) -1.3333 (21,18) -0.0000 (22,18) -0.0000 (23,18) -0.6667 (13,19) 0.3333 (14,19) 0.3333 (15,19) 0.3333 (18,19) -1.0000 (19,19) 4.3333 (20,19) -1.3333 (23,19) 0.3333 (24,19) -1.3333 (25,19) -0.0000 (15,20) -0.6667 (19,20) -1.3333 (20,20) 2.6667 (23,20) -0.0000 (24,20) -0.0000 (25,20) -0.6667 (11,21) 0.0833 (12,21) 0.0833 (13,21) 0.0833 (16,21) 0.0833 (17,21) 0.0833 (18,21) -0.0000 (21,21) 0.5833 (22,21) -0.6667 (23,21) 0.1667 (13,22) -0.0000 (17,22) -1.3333 (18,22) -0.0000 (21,22) -0.6667 (22,22) 2.6667 (23,22) -0.6667 (13,23) 0.2500 (14,23) 0.0833 (15,23) 0.0833 (17,23) -0.0000 (18,23) -0.5833 (19,23) 0.0833 (20,23) -0.0000 (21,23) 0.1667 (22,23) -0.6667 (23,23) 1.5833 (24,23) -0.6667 (25,23) 0.1667 (15,24) -0.0000 (19,24) -1.3333 (20,24) -0.0000 (23,24) -0.6667 (24,24) 2.6667 (25,24) -0.6667 (15,25) 0.1667 (19,25) -0.0000 (20,25) -0.6667 (23,25) 0.1667 (24,25) -0.6667 (25,25) 1.0000
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
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

Categories

Find more on Get Started with MATLAB in Help Center and File Exchange

Tags

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!