MATLAB projects for Computer Science/Engineering using these numerical method techniques?

i am reciving one error and could not solve it, could you please check it. many thanks

the error is

Index exceeds array bounds.

Error in probblem51 (line 53)

Sv(j)=-1*(phiold(j,i+1)+phi(j,i-1));

code :

%A two dimensional inviscid, incompressible fluid is flowing steadily
%through a chamber between inlet and the outlet. Obtain the solution
% using
%both PSOR and the ADI methods.%

close all, clc

%Defining the constants
L=6;                % length
H=4;                % Height
JN=21;              % Maximum number of grid points along y
IM=31;              % Maximum number of grid points along x
phi=zeros(JN,IM);
deltax=0.2;
deltay=0.2;
B=deltax/deltay;

%Initial Guess
for i=1:IM
for j=1:JN
phi(j,i)=50;
end
end

Errormax=0.01;
er=1;
count=0;
while er
    phiold(j,i)=phi(j,i);
    count=count+1;
    Beta=B^2;
    alpha=-2*(1+Beta);
    e=ones(IM,1);
    f=ones(JN,1);
    
    for i=2:IM-2
        for j=2:JN-1
            Su(i)=-Beta*(phiold(j+1,i)+phi(j-1,i));        
                if i==IM-1
                    Su(i)=-Beta*(phiold(j+1,i)+phi(j-1,i))-phiold(j,IM);
                end
                if i==2
                    Su(i)=-Beta*(phiold(j+1,i)+phi(j-1,i))-phiold(j,1);
                end
        end
            A=spdiags([e alpha*e e], -1:1,IM-2,IM-2);
            phiold=Su(i)*inv(A);

    end

    for j=2:JN-2
        for i=2:IM-1
                Sv(j)=-1*(phiold(j,i+1)+phi(j,i-1));
                if j==JN-1
                    Sv(j)=-1*(phiold(j,i+1)-phi(j,i-1))-phiold(JN,i);
                end
                if j==2
                    Sv(j)=-1*(phiold(j,i+1)-phi(j,i-1))-phiold(1,i);
                end
        end

        B=spdiags([f*Beta alpha*f Beta*f], -1:1, JN-2,JN-2);
        phi=Sv(j)*inv(B);
    end

 
end

[Walter Roberson]

    for i=2:IM-2
        for j=2:JN-1
            Su(i)=-Beta*(phiold(j+1,i)+phi(j-1,i));        
                if i==IM-1
                    Su(i)=-Beta*(phiold(j+1,i)+phi(j-1,i))-phiold(j,IM);
                end
                if i==2
                    Su(i)=-Beta*(phiold(j+1,i)+phi(j-1,i))-phiold(j,1);
                end
        end
            A=spdiags([e alpha*e e], -1:1,IM-2,IM-2);
            phiold=Su(i)*inv(A);
    end

You are overwriting all of phiold for each i iteration. Are you sure that is appropriate?

You are overwriting Su(i) for each different j iteration. Are you sure that is appropriate?

[mehmet salihi]

to be honest i am not sure, but there is not problem with Su, the problem is with Sv. do you have any idea about this

[Walter Roberson]

The error message is on a line that indexes into phiold and phi. The first double loop repeatedly overwrites all of phiold, so it is appropriate to first prove that the way it is overwritten is correct, that the size of phiold is kept consistent. Only then is it appropriate to study the structure of the second double loop, which overwrites all of phi, to see whether the size of phi is being kept consistent.

Hint: NO, the size of phiold is not being kept consistent.

[mehmet salihi]

i see , thank you for the explanation, i will recheck the structure of the topic. thanks again

Answer 1

[Walter Roberson]

phiold is initialized as 21 x 31, by way of the statement

    phiold(j,i)=phi(j,i);

which, the first time through, relies upon the values of j and i that happen to be "left over" after the loops

for i=1:IM
for j=1:JN
phi(j,i)=50;
end
end

MATLAB does define that after a for loop that does not re-assign the value of a loop index, that afterwards the loop variable will be left as the last value it was assigned... so after the double loop, i will be IM and j will be JN as those are the last things those were assigned. But it is not good practice to count on that; it makes it difficult to read and debug. If you are wanting to initialize phiold as IM by JN, then initialize it that size, instead of counting on left-over j and i values.

Later you

A=spdiags([e alpha*e e], -1:1,IM-2,IM-2);

which defines A as being an IM-2 by IM-2 matrix -- so 29 x 29.

phiold=Su(i)*inv(A);

That overwrites the 21 x 31 phiold with a 29 x 29 phiold, since Su(i) is a scalar and inv(A) is the same size as A which is 29 x 29.

    for j=2:JN-2
        for i=2:IM-1
                Sv(j)=-1*(phiold(j,i+1)+phi(j,i-1));
 

There you access up to i+1 as the second index of phiold, and i can be at most IM-1, so you will be indexing up to IM as the second index. But phiold is not IM wide: it is IM-2 wide, because A is IM-2 wide and constant times inv(A) is going to be IM-2 x IM-2.

And you are expecting the first index of phiold to be at most JN-2, so at most 21-2 = 19, even though phiold is now 29 x 29.

    for j=2:JN-2
        for i=2:IM-1
                Sv(j)=-1*(phiold(j,i+1)+phi(j,i-1));
                if j==JN-1
                    Sv(j)=-1*(phiold(j,i+1)-phi(j,i-1))-phiold(JN,i);
                end

The code in that if does access up to phiold(JN,i) which would be phiold(21,i) . However, it is conditional upon j==JN-1 in a circumstances that for j=2:JN-2 so j is at most JN-2 and so j==JN-1 will never be true.So that test will always be false and accessing at JN does not apply.

The way you access phiold in the loops dealing with Sv strongly suggest that the code thinks phiold is 21 x 31, rather than the 29 x 29 that you calculate by way of the spdiags and inv(A) .