filling a matrix using for loops

Question

joshua payne 2023-2-24

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1918485-filling-a-matrix-using-for-loops

编辑： joshua payne 2023-2-24

 clc
clear all
Rmax=10;
Rcyl=.5;
Uinf=1;
Nr=50;
Ntheta=50;
theta=linspace(0,2*pi,Ntheta)
%discretize
dr=(Rmax-Rcyl)/(Nr-1);
dtheta= 2*pi/(Ntheta-1);
rj=zeros(Nr,1); %allocating rj vector
rtheta=zeros(Ntheta,1); %allocating rtheta vector
for jj=1:Nr
    rj(jj)=Rcyl+(jj-1)*dr;
end
for kk=1:Ntheta
    rtheta(kk)=(kk-1)*dtheta;
end
w=2/(1+sqrt(1-cos(pi/(Nr+1))*cos(pi/(Ntheta+1))));
term=((2./dr.^2)+(2./(rj.^2).*dtheta.^2)); %the term with psi on left side in eqn8
psi=zeros(Ntheta,Nr); %old psi
psi(1:Ntheta,1)=0; %BC's on wall side
psi(1:Ntheta,Nr)=Uinf*Rmax*sin(theta); %BC's on far field
psi(Ntheta,1:Nr)=0 ;%BC's on top
psi(1,1:Nr)= 0 ;%BC's on bottom
%SOR formula for new psi
psinew=zeros(Ntheta,Nr); %allocating new psi values
bracket=psi
for j=2:Nr-1
    for k=2:Ntheta-1
        psi(j,k)=((1./rj(j)).*((psi(j+1,k)-psi(j-1,k))./(2.*dr))+((psi(j+1,k)-psi(j-1,k))./(dr.^2))+(1./(rj(j).^2)).*((psi(j,k+1)+psi(j,k+1)))./(dtheta.^2)).*term;
        
        deltapsi=bracket/((2/(dr^2))+(2./((rj.^2).*dtheta^2)))-psi(j,k)
        psinew(j,k)=psi(j,k)+w.*deltapsi
    end
end

im trying build a psinew matrix, based off psi(j,k) where psi(j,k) becvomes the bracket term. i get the error that "Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 50-by-1." psi matrix should be a 50 by 50

any help is appreciated

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Walter Roberson 2023-2-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1918485-filling-a-matrix-using-for-loops#answer_1179385

在 MATLAB Online 中打开

Rmax=10;
Rcyl=.5;
Uinf=1;
Nr=50;
Ntheta=50;
theta=linspace(0,2*pi,Ntheta)
theta = 1×50
         0    0.1282    0.2565    0.3847    0.5129    0.6411    0.7694    0.8976    1.0258    1.1541    1.2823    1.4105    1.5387    1.6670    1.7952    1.9234    2.0517    2.1799    2.3081    2.4363    2.5646    2.6928    2.8210    2.9493    3.0775    3.2057    3.3339    3.4622    3.5904    3.7186
%discretize
dr=(Rmax-Rcyl)/(Nr-1);
dtheta= 2*pi/(Ntheta-1);
rj=zeros(Nr,1); %allocating rj vector

That is a vector

rtheta=zeros(Ntheta,1); %allocating rtheta vector
for jj=1:Nr
    rj(jj)=Rcyl+(jj-1)*dr;
end
for kk=1:Ntheta
    rtheta(kk)=(kk-1)*dtheta;
end
w=2/(1+sqrt(1-cos(pi/(Nr+1))*cos(pi/(Ntheta+1))));
term=((2./dr.^2)+(2./(rj.^2).*dtheta.^2)); %the term with psi on left side in eqn8

That uses all of rj, and rj is a vector, so term is non-scalar

psi=zeros(Ntheta,Nr); %old psi
psi(1:Ntheta,1)=0; %BC's on wall side
psi(1:Ntheta,Nr)=Uinf*Rmax*sin(theta); %BC's on far field
psi(Ntheta,1:Nr)=0 ;%BC's on top
psi(1,1:Nr)= 0 ;%BC's on bottom
%SOR formula for new psi
psinew=zeros(Ntheta,Nr); %allocating new psi values
bracket=psi
bracket = 50×50
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0     0
whos
  Name          Size            Bytes  Class     Attributes

  Nr            1x1                 8  double              
  Ntheta        1x1                 8  double              
  Rcyl          1x1                 8  double              
  Rmax          1x1                 8  double              
  Uinf          1x1                 8  double              
  bracket      50x50            20000  double              
  cmdout        1x33               66  char                
  dr            1x1                 8  double              
  dtheta        1x1                 8  double              
  jj            1x1                 8  double              
  kk            1x1                 8  double              
  psi          50x50            20000  double              
  psinew       50x50            20000  double              
  rj           50x1               400  double              
  rtheta       50x1               400  double              
  term         50x1               400  double              
  theta         1x50              400  double              
  w             1x1                 8  double              
for j=2:Nr-1
    for k=2:Ntheta-1
        psi(j,k)=((1./rj(j)).*((psi(j+1,k)-psi(j-1,k))./(2.*dr))+((psi(j+1,k)-psi(j-1,k))./(dr.^2))+(1./(rj(j).^2)).*((psi(j,k+1)+psi(j,k+1)))./(dtheta.^2)).*term;
 

That uses all of term on the right-hand side, and term is a vector, so the right hand side is a vector, but the left-hand side names a scalar destination.

       
        deltapsi=bracket/((2/(dr^2))+(2./((rj.^2).*dtheta^2)))-psi(j,k)
        psinew(j,k)=psi(j,k)+w.*deltapsi
    end
end
Unable to perform assignment because the size of the left side is 1-by-1 and the size of the right side is 50-by-1.

3 个评论
显示 1更早的评论隐藏 1更早的评论

Walter Roberson 2023-2-24

在 MATLAB Online 中打开

psi(j,k)=((1./rj(j)).*((psi(j+1,k)-psi(j-1,k))./(2.*dr))+((psi(j+1,k)-psi(j-1,k))./(dr.^2))+(1./(rj(j).^2)).*((psi(j,k+1)+psi(j,k+1)))./(dtheta.^2)).*term;

j and k are scalars

rj(j) is a scalar so 1./rj(j) is a scalar.

psi(j+1,k) is a scalar, and so is psi(j-1,k) so the subtraction is a scalar.

dr is a scalar.

So ((1./rj(j)).*((psi(j+1,k)-psi(j-1,k))./(2.*dr)) is operations only on scalars, giving a scalar result.

psi(j+1,k) is a scalar, as is psi(j-1,k) so the subtraction is a scalar. dr is a scalar. rj(j) is a scalar. So everything in ((psi(j+1,k)-psi(j-1,k))./(dr.^2))+(1./(rj(j).^2)) is operations on scalars, so that part is a scalar.

Likewise everything in (1./(rj(j).^2)).*((psi(j,k+1)+psi(j,k+1)))./(dtheta.^2)) is operations on scalars.

So every calculation on the line is a scalar. Until you get to the .*term part. term is a 50x1 vector. So you are multiplying some calculated scalar by a vector, which is going to give you a vector result. And you are trying to store that vector result into a scalar output location.

Perhaps you should be indexing term by something in the calculation, so that you get a scalar for that part?

joshua payne 2023-2-24

编辑：joshua payne 2023-2-24

oh that makes sense.

i know my goal is is use the rj thats a vector then multiply its elements by the psi matrix elements and fill it using the for loops. i do not want anything to be a scalar. im not sure how to go about it then but maybe i require more knowledge

请先登录，再进行评论。