filling a matrix using for loops

 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

Answer 1

[Walter Roberson]

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.

[joshua payne]

im sorry, but my goal is to not name a scalar, i want calculate the elments of psi for j and k then allocate them to a matrix. correpsonding to the j,k element

[Walter Roberson]

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]

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