Why am I not able to use decimal numbers in my code? (Index in position 1 is invalid. Array indices must be positive integers or logical values)

I'm getting the error "Index in position 1 is invalid. Array indices must be positive integers or logical values", meaning that my numbers with decimals won't work (which I need), as in the rho, and some others I need to change, and I'm kind of unsure of what to change.

clc
clear all;
close all;

loop=1;
Nx=50;
Ny=50;
M=4000; 
%--------------------------------------------------------------------------
%--Initialize V-matrix and Rho-matrix--------------------------------------
%---V(i,j)=potential inside the 2-D box------------------------------------
%---rho(i,j)=charge density inside the 2-D box-----------------------------
%---i is along X-axis and j is along Y-axis--------------------------------
%--------------------------------------------------------------------------
V(1:Nx,1:Ny)=0.0;
rho(1:Nx,1:Ny)=0.0;
rho(0.05,0.025)=0; % charge at the center of the box
%--------------------------------------------------------------------------
% Boundary conditions
V(1,:)=0; 
V(Nx,:)=1; 
V(:,1)=0; 
V(:,Ny)=0; 

 X=1:Nx;
 Y=1:Ny;
 contour(X,Y,V,0:1:50,'linewidth',2)
 h=gca; 
get(h,'FontSize') 
set(h,'FontSize',12)
colorbar('location','eastoutside','fontsize',12);
axis([1 Nx 1 Ny])
xlabel('X','fontSize',12);
ylabel('Y','fontSize',12);
title('Electric Potential Distribution, V(X,Y)','fontsize',12);
fh = figure(1);
set(fh, 'color', 'white'); 

[Walter Roberson]

Are you sure you want to go as far out as rho(x=50, y=50) = 0.0 when you are also working at the scale of rho(x=0.05, y=0.025) = 0 ? In order to place everything at grid coordinates, you would have to use indexes rho(x,y) -> Rho(row = 40*x, column = 40*y) which would take you out to row 2000 and column 2000, and rho(x=0.05, y=0.025) would then map to Rho(row=2, column=1) .

But your comment puts the charge (0) "at the center of the box", implying that the box would have to extend from about x=-50 to x=+50 and y=-50 to y=+50, which would require using grid coordinates such as Rho(row=80*x+4001, column=80*y+4001) for a total of 8001 rows and columns.... are you sure that is what you want?

[Nicholas Gillam]

I haven't changed that part yet since I'm stuck trying to get things to non integer values, but I need the values of the boundaries to be x=y=0.1.

Answer 1

[Star Strider]

This is the assignment that throws the error:

rho(0.05,0.025)=0; % charge at the center of the box

Note that 0.05 and 0.025 are not integers.

Since ‘rho’ is a (50x50) double array, perhaps you intend:

rho(25, 25)=0; % charge at the center of the box

With that change, your code runs without error, although the result is likely less than what you want.

[Nicholas Gillam]

I actually need non integer values (hence the 0.05, 0.025). I know I have been able to do it before on similar code (which I so convientiently don't have anymore), and can't remember how to do it.

[Walter Roberson]

The only way it could have worked is if you had defined a custom storage class that permitted non-integer indices.

Answer 2

[Walter Roberson]

factor = 40;
loop=1;
Nx=50;
Ny=50;
M=4000; 
%--------------------------------------------------------------------------
%--Initialize V-matrix and Rho-matrix--------------------------------------
%---V(i,j)=potential inside the 2-D box------------------------------------
%---rho(i,j)=charge density inside the 2-D box-----------------------------
%---i is along X-axis and j is along Y-axis--------------------------------
%--------------------------------------------------------------------------
V(1:Nx*factor,1:Ny*factor)=0.0;
rho(1:Nx*factor,1:Ny*factor)=0.0;
rho(0.05*factor,0.025*factor)=0; % charge at the center of the box
%--------------------------------------------------------------------------
% Boundary conditions
V(1,:)=0; 
V(Nx*factor,:)=1; 
V(:,1)=0; 
V(:,Ny*factor)=0; 
 X=1:Nx*factor;
 Y=1:Ny*factor;
 contour(X,Y,V,0:1:50,'linewidth',2)
 h=gca; 
get(h,'FontSize') 
set(h,'FontSize',12)
colorbar('location','eastoutside','fontsize',12);
axis([1 Nx 1 Ny])
xlabel('X','fontSize',12);
ylabel('Y','fontSize',12);
title('Electric Potential Distribution, V(X,Y)','fontsize',12);
fh = figure(1);
set(fh, 'color', 'white'); 

The contour will be empty because almost everything is 0 and the exceptions that are 1 are only directly against the edges. But at least it will not give invalid index.

[Nicholas Gillam]

I have seen what I am supposed to get, and there will be contours along the region gradually going down from 1 to 0 along the y-axis.

[Walter Roberson]

You assign 0 to all locations in V except for V(Nx*factor,:)=1; which is to say the bottom row. Nothing else is non-zero. There is no gradient to be found, just a large matrix of all 0 that suddenly jumps to 1. You ask for contours at 0 and contours at 1. The contour for 0 is going to occupy everything except the bottom margin. The contour at 1 is going to be a thin line across the bottom.

[Nicholas Gillam]

I'm not sure I understand your wording right. The graph should look almost identical to this (with more contours for a more accurate graph):

I am doing a different approach to the same problem (via FDM). I know that there will be a gradient as well, as the electric field is found via a gradient (and I also know before hand that it exists).

[Walter Roberson]

Review your code:

V(1:Nx*factor,1:Ny*factor)=0.0;
%above makes V Nx*factor talll and Ny*factor wide, and initializes it all to 0

V(1,:)=0; 
%above sets top edge of V to 0, but it was already 0

V(Nx*factor,:)=1; 
%above sets bottom edge of V to 1, but only one pixel thick

V(:,1)=0;
%above sets left edge of V to 0, but it was already 0
%except for V(Nx*factor,1) which was set to 1 on previous line

V(:,Ny*factor)=0;
%above sets right edge of V to 0, but it was already 0
%except for V(Nx*factor,Ny*factor) which was set to 1 two lines ago

 X=1:Nx*factor;
 %above does not change V
 
 Y=1:Ny*factor;
 %above does not change V
 
 contour(X,Y,V,0:1:50,'linewidth',2)
 %above does a contour based on V

So, what part of V is non-zero? The answer is that only the bottom row is non-zero (it is 1), and even then the left and right edges are 0.

Now, you might make a claim that contour should see the 0 beside the bottom row and the 1 in the bottom row, and make a gradient between those two. But if so then

1. The entire gradient would have to occur in the space between the second last row and the last row as all of the rest of V has been set to definite 0; and
2. Any such gradient would be a values of V between 0 and 1, but your contour asks to plot contour values only at values 0, 1, 2, ... 50.

The 0:1:50 you passed to contour() is the contour levels to plot. You ask to plot up to 50, but the largest value you store in the array is 1.

There is simply no way at all that you should expect a contour like that given your code.

What could be the case is that the values might fill out if you were to run a number of steps of some finite element mesh code. For example, hypothetically the bottom edge could be programmed as a heat source, in which case after a number of steps of heat distribution, the material could have heated up a number of degrees in a way that matched the pattern you are expecting.

[Nicholas Gillam]

You know it just now also occured to me (enough that I feel so incredibly stupid now) that my actual code for the FDM is not in here. Got to fix that first.

[Nicholas Gillam]

Added the FDM method code, though I am still stuck with my initial source of my post, and added some comments

I also switched the j and i terms so I don't have to invert matrices.

clc
clear all;
close all;
%-------------------------------------------------------------------------
%--Define dimension of the 2-d box and -----------------------------------
%-------------------------------------------------------------------------
loop=1;
Nx=0.1;
Ny=0.05;
M=1000; % maximum iteration value
%--------------------------------------------------------------------------
%--Initialize V-matrix and Rho-matrix--------------------------------------
%---V(i,j)=potential inside the 2-D box------------------------------------
%---rho(i,j)=charge density inside the 2-D box-----------------------------
%---i is along X-axis and j is along Y-axis--------------------------------
%--------------------------------------------------------------------------
V(1:Nx,1:Ny)=0.0;
rho(1:Nx,1:Ny)=0.0;
rho(0.05,0.025)=0.0; % charge at the center of the box
%--------------------------------------------------------------------------
% Boundary conditions
V(1,:)=0; 
V(Nx,:)=1; 
V(:,1)=0; 
V(:,Ny)=0; 
%-------------------------------------------------------------------------
w=cos(pi/Nx)+cos(pi/Ny); % Converging Term
Ncount=0;
loop=1;
while loop==1;
   Rmin=0; 
for j=2:Nx-1
for i=2:Ny-1
   Residue=w.*(0.25.*(V(j-1,i)+V(j+1,i)+V(j,i-1)+V(j,i+1)+rho(j,i))-V(j,i));
   Rmin=Rmin+abs(Residue);
   V(j,i)=V(j,i)+Residue;
end
end
end
Rmin=Rmin/(Nx*Ny); % Average Residue per grid point

%------------------------------------------------------------------------
% Plot the result
%------------------------------------------------------------------------
 X=1:Nx;
 Y=1:Ny;
 contour(X,Y,V,0:1:50,'linewidth',2)
 h=gca; 
get(h,'FontSize') 
set(h,'FontSize',12)
colorbar('location','eastoutside','fontsize',12);
axis([1 Nx 1 Ny])
xlabel('X','fontSize',12);
ylabel('Y','fontSize',12);
title('Electric Potential Distribution, V(X,Y)','fontsize',12);
fh = figure(1);
set(fh, 'color', 'white'); 

[Walter Roberson]

Your "while loop==1;" is an infinite loop.

[Walter Roberson]

I had to take some guesses because your coordinates were confused. Some places you were using 1 to mean first index, but other places you were using factional values between 0 and 1 to indicate coordinates, so I had to guess about which coordinates were intended to be absolute and which were intended to be scaled.

The resulting plot is not what you were hoping for, but at least it gets closer.

Over 1000 iterations, the maximum value V gets to is 1.0, but V also goes negative. If you went to 10000 iterations then somewhere in there V would creep up to 1.008 . It would take a long long time before it reached 50 (if it can at all.)

clearvars

%-------------------------------------------------------------------------
%--Define dimension of the 2-d box and -----------------------------------
%-------------------------------------------------------------------------
factor = 4000;
coord2ind = @(c) round(c*factor);
nx = 0.1;
ny = 0.05;
Nx = coord2ind(nx);
Ny = coord2ind(ny);
M=1000; % maximum iteration value
%--------------------------------------------------------------------------
%--Initialize V-matrix and Rho-matrix--------------------------------------
%---V(i,j)=potential inside the 2-D box------------------------------------
%---rho(i,j)=charge density inside the 2-D box-----------------------------
%---i is along X-axis and j is along Y-axis--------------------------------
%--------------------------------------------------------------------------
V(1:Nx, 1:Ny)=0.0;
rho(1:Nx, 1:Ny)=0.0;
rho(coord2ind(0.05), coord2ind(0.025))=0.0; % charge at the center of the box
%--------------------------------------------------------------------------
% Boundary conditions
V((1),:)=0;
V((Nx),:)=1;
V(:,(1))=0;
V(:,(Ny))=0;
%-------------------------------------------------------------------------
w=cos(pi/nx)+cos(pi/ny); % Converging Term
Rmins = zeros(1,M);
for loop=1:M
    Rmin=0;
    for j=2:Nx-1
        for i=2:Ny-1
            Residue=w.*(0.25.*(V(j-1,i)+V(j+1,i)+V(j,i-1)+V(j,i+1)+rho(j,i))-V(j,i));
            Rmin=Rmin+abs(Residue);
            V(j, i) = V(j, i)+Residue;
        end
    end
    Rmins(loop) = Rmin;
end
Rmins = Rmins/(Nx*Ny); % Average Residue per grid point, recorded by iteration
%------------------------------------------------------------------------
% Plot the result
%------------------------------------------------------------------------
X = (1:Nx);
Y = (1:Ny);

%contour(X, Y, V.', 0:1:50, 'linewidth', 2)
contour(Y, X, V, 10, 'linewidth', 2)

h=gca;
set(h,'FontSize',12)
colorbar('location','eastoutside','fontsize',12);
axis(h, [1 Nx 1 Ny])
xlabel(h, 'X','fontSize',12);
ylabel(h, 'Y','fontSize',12);
title(h, 'Electric Potential Distribution, V(X,Y)','fontsize',12);
fh = figure(1);
set(fh, 'color', 'white');

[Nicholas Gillam]

nx = 0.1 (x-axis maximum) and ny = 0.05 (y-axis maximum) are the boundaries in which everything happens, the scaling of the code needs to occur such that these conditions are met. While I'm at least getting something somewhat closer to what I'm looking for, I have changed a bit of the code you provided in attempt to make it look like what it should. However, the "factor" in line 5 must be greater than or equal to 100, or otherwise I'm getting the same error as before. Setting it to 100 has gotten closer to what I need for the boundaries, however as I just noted, it cannot go lower than 100. Attempts at making nx and ny different values also gives me the same error as well.

[Walter Roberson]

You have a grid of discrete points, each of fixed width and height. Somewhere in that grid is a point that corresponds to x=0.05, y=0.025, but that does not mean that the point should be indexed with those coordinates.

As your maximum is x=0.1 and y=0.05, then your statement

V(0.05,0.025)=0.0;

should instead become

V(floor(end/2), floor(end/2)) = 0.0;

and X should not be 1:Nx and Y should not be 1:Ny : instead you should

Nx = 100;
Ny = 100;
maxx = 0.1;
maxy = 0.05;
X = linspace(0, maxx, Nx);
Y = linspace(0, maxy, Ny);

and now maxx/Nx is the x resolution and maxy/Ny is the y resolution. If you want square "pixels" then you should define Nx with twice the value of Ny.

[Nicholas Gillam]

It is now giving me an error that "undefined function or variable 'V'', but I'm defining that there am I not?

[Walter Roberson]

Attached.

But are you sure that you want to put a charge of 0 at the center point? Why bother, considering that the array starts all 0 (except the edges) ?

[Nicholas Gillam]

I just put it in the middle cause I wasn't too sure about how the FDM would handle it. Solving the problem by doing normal methods for PDEs, I didn't really need to specify, so you're right by saying why bother. While the graph shown is showing more of what I need, it's still not it, so I'm going to see what I can do with it. Though I'm kind of stumped onto why I'm getting almost no data for x>0.05. Thanks for helping and being patient with me.