I got the error: Subscript indices must either be real positive integers or logicals. Could I have your help please?

Hi supporters.

I have the error mentioned above in line 33 of the following script

function f=femRHS2d(x,y,u)
element_order = 3;
xmin = 0;  xmax=1; ymin = 0; ymax = 1;
nx=3;
ny=3;
node_num = nx * ny;
element_num = 2 * ( nx - 1 ) * ( ny - 1 );
node_xy = zeros(2 , node_num);
element_node = zeros(element_order , element_num);
v=zeros(node_num , 1);
u_old=ones(node_num , 1);
%
k = 0; 
for j = 1 : ny
    for i = 1 : nx
        
        k = k + 1;
        
        node_xy(1,k) = ( ( nx - i ) * xmin + ( i - 1 ) * xmax ) / ( nx - 1 );
        
        node_xy(2,k) = ( ( ny - j ) * ymin + ( j - 1 ) * ymax ) / ( ny - 1 );
        
    end
end
%
for e = 1 : element_num
    
    i1 = element_node(1,e);
    i2 = element_node(2,e);
    i3 = element_node(3,e);
    
    32-    area = 0.5 * ...
        33-      ( node_xy(1,i1) * ( node_xy(2,i2) - node_xy(2,i3) ) ...
        34-     + node_xy(1,i2) * ( node_xy(2,i3) - node_xy(2,i1) ) ...
        35-      + node_xy(1,i3) * ( node_xy(2,i1) - node_xy(2,i2) ) );
    
    %
    %  Consider each quadrature point.
    %  Here, we use the midside nodes as quadrature points.
    %
    for q1 = 1 : 3
        
        q2 = mod ( q1, 3 ) + 1;
        
        nq1 = element_node(q1,e);
        nq2 = element_node(q2,e);
        
        xq = 0.5 * ( node_xy(1,nq1) + node_xy(1,nq2) );
        yq = 0.5 * ( node_xy(2,nq1) + node_xy(2,nq2) );
        wq = 1.0 / 3.0;
        %
        %  Consider each test function in the element.
        %
        for ti1 = 1 : element_order
            
            ti2 = mod ( ti1,     3 ) + 1;
            ti3 = mod ( ti1 + 1, 3 ) + 1;
            
            nti1 = element_node(ti1,e);
            nti2 = element_node(ti2,e);
            nti3 = element_node(ti3,e);
            
            qi = 0.5 * ( ...
                (node_xy(1,nti3)-node_xy(1,nti2))*(yq-node_xy(2,nti2)) ...
                -(node_xy(2,nti3)-node_xy(2,nti2))*(xq-node_xy(1,nti2))) ...
                / area;
            dqidx = - 0.5 * ( node_xy(2,nti3) - node_xy(2,nti2) ) / area;
            dqidy =   0.5 * ( node_xy(1,nti3) - node_xy(1,nti2) ) / area;
            
            v  = rhs(u_old);
            
            f(nti1) = f(nti1) + area * wq * rhs * qi;
            
        end
        
    end
    
end

Answer 1

[Torsten]

You initialized "element_node" as a matrix of zeros. That's why you try to access node_xy(1,0),node_xy(2,0) in the area calculation.

But MATLAB indexing starts at 1, not 0, which results in the error message.

[Qais Al Faraei]

Thank you Mr. Torsten

I got your idea. But do you mean I have to change the element_node matrix from zeroes to something else?

[Torsten]

If I knew what the matrix "element_node" is supposed to do, maybe I could answer.

[Qais Al Faraei]

If you could help me, here is the full original code where I took it from

function fem2d_poisson_rectangle_linear ( nx, ny )

%*****************************************************************************80

%

%% fem2d_poisson_rectangle_linear() solves the Poisson equation in a rectangle.

%

% Discussion:

%

% This program solves

%

% - d2U(X,Y)/dx2 - d2U(X,Y)/dy2 = F(X,Y)

%

% in a rectangular region in the plane.

%

% Along the boundary of the region, Dirichlet conditions

% are imposed:

%

% U(X,Y) = G(X,Y)

%

% The code uses continuous piecewise linear basis functions on

% triangles determined by a uniform grid of NX by NY points.

%

% u = sin ( pi * x ) * sin ( pi * y ) + x

%

% dudx = pi * cos ( pi * x ) * sin ( pi * y ) + 1

% dudy = pi * sin ( pi * x ) * cos ( pi * y )

%

% d2udx2 = - pi * pi * sin ( pi * x ) * sin ( pi * y )

% d2udy2 = - pi * pi * sin ( pi * x ) * sin ( pi * y )

%

% rhs = 2 * pi * pi * sin ( pi * x ) * sin ( pi * y )

%

% THINGS YOU CAN EASILY CHANGE:

%

% 1) Change NX or NY, the number of nodes in the X and Y directions.

% 2) Change XL, XR, YB, YT, the left, right, bottom and top limits

% of the rectangle.

% 3) Change the exact solution in the EXACT routine, but make sure you also

% modify the formula for RHS in the assembly portion of the program.

%

% HARDER TO CHANGE:

%

% 4) Change from "linear" to "quadratic" triangles;

% 5) Change the region from a rectangle to a general triangulated region;

% 6) Handle Neumann boundary conditions.

%

% Licensing:

%

% This code is distributed under the GNU LGPL license.

%

% Modified:

%

% 01 November 2010

%

% Author:

%

% John Burkardt

%

% Input:

%

% integer NX, NY, the number of nodes in the X and Y directions.

%

% Local:

%

% sparse real A(NODE_NUM,NODE_NUM), the finite element system matrix.

%

% real B(NODE_NUM), the finite element right hand side.

%

% real C(NODE_NUM), the finite element coefficient vector.

%

% integer ELEMENT_NODE(3,ELEMENT_NUM), the indices of the nodes

% that form each element.

%

% integer ELEMENT_NUM, the number of elements.

%

% integer NODE_NUM, the number of nodes.

%

% real NODE_XY(2,NODE_NUM), the X and Y coordinates of each node.

%

% real XL, the X coordinate of the left boundary.

%

% real XR, the X coordinate of the right boundary.

%

% real YB, the Y coordindate of the bottom boundary.

%

% real YT, the Y coordinate of the top boundary.

%

timestamp ( );

element_order = 3;

xl = 0.0;

xr = 1.0;

yb = 0.0;

yt = 1.0;

fprintf ( 1, '\n' );

fprintf ( 1, 'FEM2D_POISSON_RECTANGLE_LINEAR\n' );

fprintf ( 1, ' MATLAB/Octave version %s\n', version ( ) );

fprintf ( 1, '\n' );

fprintf ( 1, ' Solution of the Poisson equation:\n' );

fprintf ( 1, '\n' );

fprintf ( 1, ' - Uxx - Uyy = F(x,y) inside the region,\n' );

fprintf ( 1, ' U(x,y) = G(x,y) on the boundary of the region.\n' );

fprintf ( 1, '\n' );

fprintf ( 1, ' The region is a rectangle, defined by:\n' );

fprintf ( 1, '\n' );

fprintf ( 1, ' %g = XL<= X <= XR = %g\n', xl, xr );

fprintf ( 1, ' %g = YB<= Y <= YT = %g\n', yb, yt );

fprintf ( 1, '\n' );

fprintf ( 1, ' The finite element method is used, with piecewise\n' );

fprintf ( 1, ' linear basis functions on 3-node triangular\n' );

fprintf ( 1, ' elements.\n' );

fprintf ( 1, '\n' );

fprintf ( 1, ' The corner nodes of the triangles are generated by an\n' );

fprintf ( 1, ' underlying grid whose dimensions are\n' );

fprintf ( 1, '\n' );

fprintf ( 1, ' NX = %d\n', nx );

fprintf ( 1, ' NY = %d\n', ny );

%

% NODE COORDINATES

%

node_num = nx * ny;

fprintf ( 1, ' Number of nodes = %d\n', node_num );

node_xy = zeros(2,node_num);

k = 0;

for j = 1 : ny

for i = 1 : nx

k = k + 1;

node_xy(1,k) = ( ( nx - i ) * xl ...

+ ( i - 1 ) * xr ) ...

/ ( nx - 1 );

node_xy(2,k) = ( ( ny - j ) * yb ...

+ ( j - 1 ) * yt ) ...

/ ( ny - 1 );

end

end

%

% ELEMENT array

%

% Organize the nodes into a grid of 3-node triangles.

%

element_num = 2 * ( nx - 1 ) * ( ny - 1 );

fprintf ( 1, ' Number of elements = %d\n', element_num );

element_node = zeros ( element_order, element_num );

k = 0;

for j = 1 : ny - 1

for i = 1 : nx - 1

%

% (I,J+1)-

% | \

% | \

% | \ |

% (I,J)---(I+1,J)

%

k = k + 1;

element_node(1,k) = i + ( j - 1 ) * nx;

element_node(2,k) = i + 1 + ( j - 1 ) * nx;

element_node(3,k) = i + j * nx;

%

% (I,J+1)--(I+1,J+1)

% | \ |

% \ |

% \ |

% -(I+1,J)

%

k = k + 1;

element_node(1,k) = i + 1 + j * nx;

element_node(2,k) = i + j * nx;

element_node(3,k) = i + 1 + ( j - 1 ) * nx;

end

end

%

% ASSEMBLE THE SYSTEM

%

% Assemble the coefficient matrix A and the right-hand side B of the

% finite element equations, ignoring boundary conditions.

%

b = zeros(node_num,1);

a = sparse ( [], [], [], node_num, node_num );

for e = 1 : element_num

i1 = element_node(1,e);

i2 = element_node(2,e);

i3 = element_node(3,e);

area = 0.5 * ...

( node_xy(1,i1) * ( node_xy(2,i2) - node_xy(2,i3) ) ...

+ node_xy(1,i2) * ( node_xy(2,i3) - node_xy(2,i1) ) ...

+ node_xy(1,i3) * ( node_xy(2,i1) - node_xy(2,i2) ) );

%

% Consider each quadrature point.

% Here, we use the midside nodes as quadrature points.

%

for q1 = 1 : 3

q2 = mod ( q1, 3 ) + 1;

nq1 = element_node(q1,e);

nq2 = element_node(q2,e);

xq = 0.5 * ( node_xy(1,nq1) + node_xy(1,nq2) );

yq = 0.5 * ( node_xy(2,nq1) + node_xy(2,nq2) );

wq = 1.0 / 3.0;

%

% Consider each test function in the element.

%

for ti1 = 1 : element_order

ti2 = mod ( ti1, 3 ) + 1;

ti3 = mod ( ti1 + 1, 3 ) + 1;

nti1 = element_node(ti1,e);

nti2 = element_node(ti2,e);

nti3 = element_node(ti3,e);

qi = 0.5 * ( ...

( node_xy(1,nti3) - node_xy(1,nti2) ) * ( yq - node_xy(2,nti2) ) ...

- ( node_xy(2,nti3) - node_xy(2,nti2) ) * ( xq - node_xy(1,nti2) ) ) ...

/ area;

dqidx = - 0.5 * ( node_xy(2,nti3) - node_xy(2,nti2) ) / area;

dqidy = 0.5 * ( node_xy(1,nti3) - node_xy(1,nti2) ) / area;

rhs = 2.0 * pi * pi * sin ( pi * xq ) * sin ( pi * yq );

b(nti1) = b(nti1) + area * wq * rhs * qi;

%

% Consider each basis function in the element.

%

for tj1 = 1 : element_order

tj2 = mod ( tj1, 3 ) + 1;

tj3 = mod ( tj1 + 1, 3 ) + 1;

ntj1 = element_node(tj1,e);

ntj2 = element_node(tj2,e);

ntj3 = element_node(tj3,e);

qj = 0.5 * ( ...

( node_xy(1,ntj3) - node_xy(1,ntj2) ) * ( yq - node_xy(2,ntj2) ) ...

- ( node_xy(2,ntj3) - node_xy(2,ntj2) ) * ( xq - node_xy(1,ntj2) ) ) ...

/ area;

dqjdx = - 0.5 * ( node_xy(2,ntj3) - node_xy(2,ntj2) ) / area;

dqjdy = 0.5 * ( node_xy(1,ntj3) - node_xy(1,ntj2) ) / area;

a(nti1,ntj1) = a(nti1,ntj1) ...

+ area * wq * ( dqidx * dqjdx + dqidy * dqjdy );

end

end

end

end

%

% BOUNDARY CONDITIONS

%

% If the K-th variable is at a boundary node, replace the K-th finite

% element equation by a boundary condition that sets the variable to U(K).

%

k = 0;

for j = 1 : ny

for i = 1 : nx

k = k + 1;

if ( i == 1 | i == nx | j == 1 | j == ny )

[ u, dudx, dudy ] = exact ( node_xy(1,k), node_xy(2,k) );

a(k,1:node_num) = 0.0;

a(k,k) = 1.0;

b(k) = u;

end

end

end

%

% SOLVE the linear system A * C = B.

%

c = a \ b;