How do i construct this matrix with else and elseif statements? Or should i be using cases?

I am trying to construct a matrix of coefficients using a storage method from Peric's book, "Computational Methods for Fluid Dynamics", using a storage location system index 'l', however I cannot get it to work! I was able to construct the matrix using a long method, and the result is below, but it only works for one matrix size... whereas the Peric method works for all...

nX = 6;
nYp = 6;
uaW = zeros(nX,nYp);
uaE = zeros(nX,nYp);
uaP = zeros(nX,nYp);
uaN = zeros(nX,nYp);
uaS = zeros(nX,nYp);
u = zeros(nX,nYp);

uaW(:,:) = 1;
uaE(:,:) = 2;
uaP(:,:) = 3;
uaN(:,:) = 4;
uaS(:,:) = 5;

%%Construct the coefficient matrix
  nj=size(u(3:nX-1,2:nYp-1),2); % number of columns
  ni=size(u(3:nX-1,2:nYp-1),1); % number of rows
  A=zeros(ni*nj);
  for i=1:ni
      for j=1:nj;
  l = (i-1)*nj + j; % storage locater
 A(l,l) =  uaP(i+2,j+1) ;  % always used

    if l-nj>=1 % if this happens, i want uaW to be used, otherwise not.

    A(l,l-nj) = -uaW(i+2,j+1) ;

    elseif l>=2 % if this happens i want uaS to also be used, otherwise not.
   A(l,l-1) = -uaS(i+2,j+1);

    elseif l+1<ni*nj % if this happens i want uaN to also be used, otherwise not.
    A(l,l+1) = -uaN(i+2,j+1);

    elseif l+nj<ni*nj % if this happens i want uaE to also be used, otherwise not.
     A(l,l+nj) = -uaE(i+2,j+1)  ;
     end
        end
    end

The matrix should look like this:

>> full(A)
ans =
     5    -3     0     0    -2     0     0     0     0     0     0     0
    -4     5    -3     0     0    -2     0     0     0     0     0     0
     0    -4     5    -3     0     0    -2     0     0     0     0     0
     0     0    -4     5    -3     0     0    -2     0     0     0     0
    -1     0     0    -4     5    -3     0     0    -2     0     0     0
     0    -1     0     0    -4     5    -3     0     0    -2     0     0
     0     0    -1     0     0    -4     5    -3     0     0    -2     0
     0     0     0    -1     0     0    -4     5    -3     0     0    -2
     0     0     0     0    -1     0     0    -4     5    -3     0     0
     0     0     0     0     0    -1     0     0    -4     5    -3     0
     0     0     0     0     0     0    -1     0     0    -4     5    -3
     0     0     0     0     0     0     0    -1     0     0    -4     5

I'd really appreciate any help.

Best Regards, Michael

Answer 1

[Sean de Wolski]

@Everyone else: wayyyy too much effort:

toeplitz([5 -4 0 0 -1 0 0 0 0 0 0 0],[5 -3 0 0 -2 0 0 0 0 0 0 0])

Answer 2

[Babak]

N=10;
A = zeros(N,N);

for i=1:N

for j=1:N

if i=j

A(i,j)=5;

end

if i+1=j

A(i,j) = 3

end

.
.
.

end

end

[Michael Moor]

Thank-you so much for the answer. However using this arrangement will not correctly account for the variation in position, i.e. for a my nX=20 and nYp =20, the matrix rlationship is: [-uaW 0 0 0 0 0 0 -uaS uaP -uaN 0 0 0 0 0 0 -uaE]

whereas in the nX = 6 nYp =6, the relationship is, [-uaW 0 0 -uaS uaP -uaN 0 0 -uaE]

For this reason, the storage locater method from Peric is used. The indexing scheme is correct, e.g. A(l,l-1), however during the loop, this will lead to negative indexing, as well as inserting values into the matrix that i do not want, e.g. the last entry in the matrix MUST be uaP, or "5" in this illustration.

Answer 3

[Matt Fig]

Here are a couple of suggestions. First the old-fashioned approach:

A = diag(5*ones(1,12)) + ...
diag(-3*ones(1,11),1) + ...
diag(-2*ones(1,8),4) + ...
diag(-4*ones(1,11),-1) + ...
diag(-1*ones(1,8),-4)

Next use SPDIAGS, since it looks like you want to end up sparse anyways.

R = repmat([-1 -4 5 -3 -2],12,1);
R([9:12 24 37 49:52]) = 0;
D = [-4 -1 0 1 4];
A = spdiags(R,D,12,12) % compare full(A) to previous A.

[Michael Moor]

Thank-you for your answer. Unfortunately though, if the matrix size changes, the spacing of the diagonals also changes. Also, I only used matrices of those numbers to highlight which matrix they came from, in my code, these numbers are all different values. This is why I have the storage locater "l". The indexing scheme is correct, e.g. A(l,l-1), however during the loop, this will lead to negative indexing, as well as inserting values into the matrix that i do not want, e.g. the last entry in the matrix MUST be uaP, or "5" in this illustration.

My issue is thus trying to program out these unwanted entries.

Best Regards, Michael

Answer 4

[Andrei Bobrov]

nX = 6;
nYp = 6;
uaW = 1;
uaE = 2;
uaP = 3;
uaN = 4;
uaS = 5;
nj=size(u(3:nX-1,2:nYp-1),2);
ni=size(u(3:nX-1,2:nYp-1),1);
A=zeros(ni*nj);
for i=1:ni
    for j=1:nj;
        l = (i-1)*nj + j;
        A(l,l) =  uaS ;
        if l >= 2
             A(l,l-1) = -uaN;
             A(l-1,l) = -uaP;
             if l-nj>=1
                 A(l,l-nj) = -uaW ;
                 A(l-nj,l) = -uaE ;
             end
        end
    end
end

OR like Matt

A = full(spdiags(repmat([-1 -4 5 -3 -2],12,1),[-4 -1:1 4],12,12));

[Michael Moor]

I think that this may be correct, I have afapted it such that the coefficients come from their respective matrix locations, and will respond as soon as i know!