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

3 visualizzazioni (ultimi 30 giorni)
岳
岳 il 21 Ott 2022
Modificato: Torsten il 21 Ott 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

Accedi per rispondere a questa domanda.

Risposte (1)

Torsten
Torsten il 21 Ott 2022
Modificato: Torsten il 21 Ott 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

Categorie

Scopri di più su Get Started with MATLAB in Help Center e File Exchange

Tag

Prodotti


Release

R2019a

Community Treasure Hunt

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

Start Hunting!

Translated by