Index in position 1 is invalid. Array indices must be positive integers or logical values.

% calculation of Emax
T_epsilon = zeros(Nx,Ny);
for i =2:Nx-1
    for j = 2:Ny-1
        T_epsilon(i,j) = T(i,j) - A((dy*i)/dy1,(dx*j)/dx1);
        E_max = max(abs(T_epsilon));
        disp(E_max);
    end
end

I used this code to calculate a difference between two matrices T of order 101-by-101 and A of order 1001-by-1001, for the corresponding value points of T matrices such that the resultant matrix will be of order 101-by-101. But I'm getting an error:

Index in position 1 is invalid. Array indices must be positive integers or logical values.

in line:T_epsilon(i,j) = T(i,j) - A((dy*i)/dy1,(dx*j)/dx1);

Where have I gone wrong?

[Torsten]

If I should make a guess: dy*i/dy1 and (dx*j)/dx1 are not integer values as required to index into a matrix.

[Kushagra Saurabh]

Yes, you are correct, so how to calculate the integer value of such variable, how to give this in code?

[Torsten]

What are the numerical values of dy/dy1 and dx/dx1 ?

[Kushagra Saurabh]

actually i want to subtract to matrices which are not of the same order though. my first matrix T is of order 101 by 101 and the other one A is 1001 by 1001 and i want the resultant matrix E_epsilon to be of order 101 by 101. I know matrix subtaction is done when the matrices is of same order but I have to calculate this E_epsilon matrix. so what i thought, is to calculate the corresponding elements similar to 101 by 101 matrix into A matrix using this A((dy*i)/dy1,(dx*j)/dx1) , (which i think is not correct though) and then subtract from T asT_epsilon(i,j) = T(i,j) - A((dy*i)/dy1,(dx*j)/dx1);

where i,j is row and elements of 101 by 101 matrix and dy=1/Ny-1 and dx=1/Nx-1 where Nx =101 = Ny and dx1 = 1/nx-1 and dy1 = 1/ny-1 where nx=ny =1001

[Torsten]

First condense A to a 100x100 matrix (e.g. by taking the mean of subsequent 10x10 matrices). Then apply your loop.

[Kushagra Saurabh]

sorry for trouble again,I tried applying 'mat2cell' function but i think i'am doing it wrong ,could you please tell how to condense A as you said?

thanks

[Torsten]

A = rand(1000);
B = zeros(100);
for i = 1:10
    for j = 1:10
        B(i,j) = mean(mean(A((i-1)*10+1:i*10,(j-1)*10+1:j*10)));
    end
end
B

B = 100×100

0.4607 0.5500 0.5175 0.4988 0.5034 0.4868 0.5293 0.4834 0.4965 0.4872 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5426 0.5242 0.4628 0.4923 0.4933 0.5141 0.4904 0.5009 0.4594 0.4817 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.4722 0.4809 0.5043 0.4865 0.4926 0.5197 0.4563 0.5304 0.4503 0.5154 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5139 0.5235 0.5194 0.5526 0.5013 0.4902 0.4870 0.4790 0.4820 0.4932 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5023 0.5050 0.4929 0.5181 0.4926 0.4986 0.5080 0.5237 0.4926 0.4807 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.4678 0.5428 0.5189 0.4717 0.5028 0.5037 0.5250 0.4999 0.5043 0.5100 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.4570 0.4781 0.5157 0.5091 0.5191 0.4759 0.4955 0.5278 0.4514 0.5040 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.4203 0.4700 0.5152 0.5164 0.4807 0.4879 0.5169 0.4948 0.5427 0.4985 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.5152 0.5349 0.4994 0.5229 0.4923 0.5498 0.5090 0.4632 0.4846 0.5030 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0.4830 0.4640 0.5185 0.5146 0.5197 0.5435 0.5107 0.4744 0.4842 0.4807 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

[Kushagra Saurabh]

thank you so much sir, i used your suggestion into my code but it's still running and had been more than 10 min..

clear all

close all

clc

% % geometry of domain

Nx = 101;

Ny = 101;

nx = 1000;

ny = 1000;

dx = 1/(Nx-1);

dy = 1/(Ny-1);

dx1 = 1/(nx-1);

dy1 = 1/(ny-1);

X = 0:dx:1;

Y = 0:dy:1;

x = 0:dx1:1;

y = 0:dy1:1;

% % initial condition

T = zeros(Nx,Ny);

T(50,20) = 2.5;

T(25,25) = -0.5;

T(75,10) = -2.5;

A = zeros(nx,ny);

A(500,200) = 2.5;

A(250,250) = -0.5;

A(750,100) = -2.5;

% % boundary condition

TL = 1;

TR = cos(6*(3*pi*Y)/2)+1;

TR1 = cos(6*(3*pi*y)/2)+1;

TT = 1;

TB = 1+X;

TB1 = 1+x;

T(:,1) = TL;

T(Ny,:) = TT;

T(:,Nx) = TR;

T(1,:) = TB;

A(:,1) = TL;

A(ny,:) = TT;

A(:,nx) = TR1;

A(1,:) = TB1;

T_new(Nx,Ny) = 0;

A_new(nx,ny) = 0;

TL = 1;

TR = cos(6*(3*pi*Y)/2)+1;

TR1 = cos(6*(3*pi*y)/2)+1;

TT = 1;

TB = 1+X;

TB1 = 1+x;

T_new(:,1) = TL;

T_new(Ny,:) = TT;

T_new(:,Nx) = TR;

T_new(1,:) = TB;

A_new(:,1) = TL;

A_new(ny,:) = TT;

A_new(:,nx) = TR1;

A_new(1,:) = TB1;

error_mag = 5;

error_req = 1e-04;

error_mag1 = 5;

error_req1 = 1e-04;

iteration = 0;

% calculation for T_new

while error_mag > error_req

for i = 2:Nx-1

for j=2:Ny-1

T_new(i,j) = (T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))/4;

T_new(50,20) = 2.5;

T_new(25,25) = -0.5;

T_new(75,10) = -2.5;

TL = 1;

TR = cos(6*(3*pi*Y)/2)+1;

TT = 1;

TB = 1+X;

T_new(:,1) = TL;

T_new(Ny,:) = TT;

T_new(:,Nx) = TR;

T_new(1,:) = TB;

iteration = iteration +1;

end

end

% calculation of error magnitude

for i= 2:Nx-1

for j = 2:Ny-1

error_mag = abs(T(i,j)-T_new(i,j));

end

end

%assigning new to old

T = T_new;

end

% Calculation for A_new

while error_mag1 > error_req1

for i1 = 2:nx-1

for j1 = 2:ny-1

A_new(i1,j1) = (A(i1+1,j1)+A(i1-1,j1)+A(i1,j1+1)+A(i1,j1-1))/4;

A_new(500,200) = 2.5;

A_new(250,250) = -0.5;

A_new(750,100) = -2.5;

TL = 1;

TR1 = cos(6*(3*pi*y)/2)+1;

TT = 1;

TB1 = 1+x;

A_new(:,1) = TL;

A_new(ny,:) = TT;

A_new(:,nx) = TR1;

A_new(1,:) = TB1;

iteration = iteration +1;

end

end

% calculation of error magnitude

for i1 = 2:nx-1

for j1 =2:ny-1

error_mag1 = abs(A(i1,j1)-A_new(i1,j1));

end

end

%assigning new to old

A = A_new;

end

% calculation of Emax

B = zeros(101,101);

for i =1:10

for j = 1:10

B(i,j) = mean(mean(A((i-1)*10+1:i*10,(j-1)*10+1:j*10)));

end

end

T_epsilon = zeros(Nx,Ny);

for i =2:Nx-1

for j = 2:Ny-1

T_epsilon(i,j) = T(i,j) - B(i,j);

E_max = max(abs(T_epsilon));

disp(E_max);

end

end

[Torsten]

What is the purpose of your code ? Do you solve Laplace's equation on two different meshes ?

[Kushagra Saurabh]

Yes, and I have to find the maximum value of truncation error...

Answer 1

[Image Analyst]

You swapped x and y. Remember x is columns which is the second index, not the first index of a matrix. Matrixes are references (row, column) or (y, x), not (x,y) which you are trying to do.

Not sure of your formula so look at it carefully but the code should look something like this:

[Ny, Nx] = size(A);  % Nx IS NOT THE FIRST RETURN ARGUMENT!!!
% calculation of Emax
T_epsilon = zeros(Nx,Ny);
for col = 2 : Nx-1
    for row = 2 : Ny-1
        T_epsilon(row, col) = T(row, col) - WHATEVER....
        E_max = max(abs(T_epsilon));
        fprintf('E_max = %f\n', E_max);
    end
end

See the FAQ for more info:

https://matlab.fandom.com/wiki/FAQ#%22Subscript_indices_must_either_be_real_positive_integers_or_logicals.%22

[Kushagra Saurabh]

I know sir, I have done previous calculations using i= Nx-1 ,j=Ny-1 and it worked well there

[Image Analyst]

I guess you know now, but you didn't when you posted the code because your indexes to T_epsilon were swapped so it would only for for a square matrix. Can you post the code after you made the corrections I mentioned? Also give code for computing dx, etc., and a .mat file with A and T in it.

I'm really confused by this statement you made "my first matrix T is of order 101 by 101 and the other one A is 1001 by 1001" and " dy=1/Ny-1 and dx=1/Nx-1 where Nx =101 = Ny and dx1 = 1/nx-1 and dy1 = 1/ny-1 where nx=ny =1001"

so if

[Ny, Nx] = size(A)

then if A is 1001 by 1001 then Nx and Ny are both 1001. So why do you say they're 101 (the same size as T).

How do you compute dx1? Like this

nx = 1001;
dx1 = 1/nx-1

dx1 = -0.9990

or this

dx1 = 1/(nx-1)

dx1 = 1.0000e-03

And since T and A are different sizes, exactly how to you plan on taking a 10x10 block of A to fit it into a single element of T (since T is 10 times smaller)???

Are you just making careless typos? If so please be more careful because that just wastes everyone's time.

[Kushagra Saurabh]

I am so sorry for my mistakes, I deeply apologize, yes dx1= 1/(nx-1) . Yes, I am confused how to fit 10×10 block of A into single element of T.

[Image Analyst]

OK, but you continue to drag this out. Did you overlook the part of my comment where I said "Can you post the code after you made the corrections I mentioned? Also give code for computing dx, etc., and a .mat file with A and T in it."?

I'll check back later for the code, with your corrected/improved terminilogy, and for the attached MAT file.

[Kushagra Saurabh]

clear all
close all
clc
% % geometry of domain
Nx = 101;
Ny = 101;
nx = 1000;
ny = 1000;
dx = 1/(Nx-1);
dy = 1/(Ny-1);
dx1 = 1/(nx-1);
dy1 = 1/(ny-1);
X = 0:dx:1;
Y = 0:dy:1;
x = 0:dx1:1;
y = 0:dy1:1;
% % initial condition
T = zeros(Nx,Ny);
T(50,20) = 2.5;
T(25,25) = -0.5;
T(75,10) = -2.5;
A = zeros(nx,ny);
A(500,200) = 2.5;
A(250,250) = -0.5;
A(750,100) = -2.5;

% % boundary condition
TL = 1;
TR = cos(6*(3*pi*Y)/2)+1;
TR1 = cos(6*(3*pi*y)/2)+1;
TT = 1;
TB = 1+X;
TB1 = 1+x;
T(:,1) = TL;
T(Ny,:) = TT;
T(:,Nx) = TR;
T(1,:) = TB;
A(:,1) = TL;
A(ny,:) = TT;
A(:,nx) = TR1;
A(1,:) = TB1;

T_new(Nx,Ny) = 0;
A_new(nx,ny) = 0;
TL = 1;
TR = cos(6*(3*pi*Y)/2)+1;
TR1 = cos(6*(3*pi*y)/2)+1;
TT = 1;
TB = 1+X;
TB1 = 1+x;
T_new(:,1) = TL;
T_new(Ny,:) = TT;
T_new(:,Nx) = TR;
T_new(1,:) = TB;
A_new(:,1) = TL;
A_new(ny,:) = TT;
A_new(:,nx) = TR1;
A_new(1,:) = TB1;

error_mag = 5;
error_req = 1e-04;
error_mag1 = 5;
error_req1 = 1e-04;
iteration = 0;

% calculation for T_new
while error_mag > error_req
    for i = 2:Nx-1
        for j=2:Ny-1
            T_new(i,j) = (T(i+1,j)+T(i-1,j)+T(i,j+1)+T(i,j-1))/4;
            T_new(50,20) = 2.5;
            T_new(25,25) = -0.5;
            T_new(75,10) = -2.5;
            TL = 1;
            TR = cos(6*(3*pi*Y)/2)+1;
            TT = 1;
            TB = 1+X;
            
            T_new(:,1) = TL;
            T_new(Ny,:) = TT;
            T_new(:,Nx) = TR;
            T_new(1,:) = TB;
            iteration = iteration +1;
        end
    end
    % calculation of error magnitude
    
    for i= 2:Nx-1
        for j = 2:Ny-1
            error_mag =  abs(T(i,j)-T_new(i,j));
        end
    end
    %assigning new to old
    T = T_new;
end



% Calculation for A_new
while error_mag1 > error_req1
    for i1 = 2:nx-1
        for j1 = 2:ny-1
            A_new(i1,j1) = (A(i1+1,j1)+A(i1-1,j1)+A(i1,j1+1)+A(i1,j1-1))/4;
            A_new(500,200) = 2.5;
            A_new(250,250) = -0.5;
            A_new(750,100) = -2.5;
            TL = 1;
            TR1 = cos(6*(3*pi*y)/2)+1;
            TT = 1;
            TB1 = 1+x;
            A_new(:,1) = TL;
            A_new(ny,:) = TT;
            A_new(:,nx) = TR1;
            A_new(1,:) = TB1;
            iteration = iteration +1;
        end
    end
    
    % calculation of error magnitude
    
    for i1 = 2:nx-1
        for j1 =2:ny-1
            error_mag1 =  abs(A(i1,j1)-A_new(i1,j1));
        end
    end
    %assigning new to old
    A = A_new;
end


% calculation of Emax
T_epsilon = zeros(Nx,Ny);
for i =2:Nx-1
    for j = 2:Ny-1
        T_epsilon(i,j) = T(i,j) - A((dy*j)/dy1,(dx*i)/dx1);
        E_max = max(abs(T_epsilon));
        disp(E_max);
    end
end

This is my entire code... sorry,I missed it earlier.

[Kushagra Saurabh]

here's .m file

[Image Analyst]

Still not exactly sure what you want to do (I didn't delve into it in detail, sorry) but if you want to process a matrix by processing it in separate blocks -- 10 vertical and 10 horizontal so you have 100 tiles -- you can use blockproc. See attached demos where I do several things with each block, like take it's mean or standard deviation, etc.

[Kushagra Saurabh]

Thank you sir

Answer 2

[Walter Roberson]

A((dy*i)/dy1,(dx*j)/dx1)

having a division in a calculation of coordinates is very often a problem. You can sometimes get away with it if you are working with pure integers and there is reason why the divisor is known to divide in exactly, such as if the divisor is 1, but if there are any non-integer quantities involved in the calculation then at some point you will encounter problems with round-off error giving you an index that is not exactly an integer.

It looks to me as if you are doing a matrix rescaling operation. If so then typically you would need to calculate each output as a linear combination of adjacent matrix elements. For example if you had coordinates (12.9, 15.4) then you would not round those to (13,15) and would instead calculate A(12,15)*0.1*0.6 + A(13,15)*0.9*0.6 + A(12,16)*0.1*0.4 + A(13,16)*0.9*0.4

[Kushagra Saurabh]

yes you're correct,but i 'am still facing a problem to put this into code. what my intension is, I want to subtract two matrices which are not of the same order though. my first matrix T is of order 101 by 101 and the other one A is 1001 by 1001 and i want the resultant matrix E_epsilon to be of order 101 by 101. I know matrix subtaction is done when the matrices is of same order but I have to calculate this E_epsilon matrix. so what i thought, is to calculate the corresponding elements similar to 101 by 101 matrix into A matrix using this A((dy*i)/dy1,(dx*j)/dx1) , (which i think is not correct though) and then subtract from T asT_epsilon(i,j) = T(i,j) - A((dy*i)/dy1,(dx*j)/dx1);

where i,j is row and elements of 101 by 101 matrix and dy=1/Ny-1 and dx=1/Nx-1 where Nx =101 = Ny and dx1 = 1/nx-1 and dy1 = 1/ny-1 where nx=ny =1001.

could you please help me here.

[Walter Roberson]

imresize() or resample()

[Walter Roberson]

A((dy*j)/dy1,(dx*i)/dx1)

That is the wrong way to process arrays; the indices are very likely not going to be integers. You could round() or floor() or ceil() to be sure they are integers, but you have to ask yourself whether that is really what you want to do.

Suppose you had an image that is all black except that every 10'th pixel horizontally and vertically is white, and you are downsizing it 10 to 1. So suppose it just happens that (dy*j)/dy1 and (dx*i)/dx1 are exact integers so that in this hypothetical situation you do not need to worry about rounding coordinates or interpolating between adjacent pixels.

Now in this hypothetical situation, would downsizing by taking every 10th pixel be the right thing to do?

You have two possibilities in the hyptothetical situation: either the integer coordinates you pick out happen to index the white pixels every time, or else the integer coordinates happen to index black pixels every time. The first would result in an all-white output: is an all-white output a fair representation of an image that is 99% black? The second would result in an all-black output: is an all-black output a fair representation of an image that regularly has some white in it?

Now suppose that you have the same input image that is white at locations where the coordinates are multiples of 10, but your downsizing happens to be picking out every 11'th pixel. (11,11) is black, (11,22) is black, (11,33) is black... eventually you get to (11,110) which might potentially be white because the 110 is a multiple of 10, but because the 11 is not a multiple of 10 it is still not white. So the output would be all black except at what was originall (110,110), (110,220), (110,330) and other locations that are multiples of 110 on both axes. Is that a fair representation?

Some methods of downsizing images do not index a single pixel as the source, and instead take the mean2() of a block of input pixels. So in the case of 10:1 reduction, with 99 of the 100 pixels in the 10 x 10 block being black and 1 being white, the output would be (W/100) where W is the intensity associated with a white pixel. If you are using uint8 where white pixels are 255, that would be (255/100) which would round to 3 (out of 255) . Almost black but slightly brighter than pure black. You would get out the average intensity of the block, which in some applications might be appropriate.

Some resizing techniques involve taking a fourier of wavelet representation of the image array and reconstructing at a subset of locations.