I am trying to understand why the dimensions of this calculation are wrong?

Question

I am solving the poisson equation and want to plot the error of the exact solution vs. number of grid points. my code is:

function [Ntot,err] = poisson(N)


nx = N;                % Number of steps in space(x)
ny = N;                % Number of steps in space(y)       
Ntot = nx*ny;
niter = 1000;           % Number of iterations 
dx = 2/(nx-1);          % Width of space step(x)
dy = 2/(ny-1);          % Width of space step(y)
x = -1:dx:1;             % Range of x(-1,1) 
y = -1:dy:1;             % Range of y(-1,1) 
b = zeros(nx,ny);       
dn = zeros(nx,ny);      


% Initial Conditions
d = zeros(nx,ny);                  
u = zeros(nx,ny);

% Boundary conditions
d(:,1) = 0;
d(:,ny) = 0;
d(1,:) = 0;                  
d(nx,:) = 0;


% Source term
b(round(ny/4),round(nx/4)) = 3000;
b(round(ny*3/4),round(nx*3/4)) = -3000;


i = 2:nx-1;
j = 2:ny-1;

% 5-point difference (Explicit)
for it = 1:niter
    dn = d;
    d(i,j) = ((dy^2*(dn(i + 1,j) + dn(i - 1,j))) + (dx^2*(dn(i,j + 1) + dn(i,j - 1))) - (b(i,j)*dx^2*dy*2))/(2*(dx^2 + dy^2));
    u(i,j) = 2*pi*pi*sin(pi*i).*sin(pi*j);

    % Boundary conditions 
    d(:,1) = 0;
    d(:,ny) = 0;
    d(1,:) = 0;                  
    d(nx,:) = 0;
end


%     
% 
% err = abs(u - d);

the error I get is:

Subscripted assignment dimension mismatch.

Error in poisson (line 39)

u(i,j) = 2*pi*pi*sin(pi*i).*sin(pi*j);

I am not sure why it is not calculating u at every grid point. I tried taking it out of the for loop but that did not help. Any ideas would be appreciated.

Answer 1

This is because i and j are both 1-by-(N-2) vectors, so u(i, j) is an (N-2)-by-(N-2) matrix. However, the expression 2*pi*pi*sin(pi*i).*sin(pi*j) is a 1-by-(N-2) vector.
The dimensions obviously don't match, hence the error.

I'm not sure, but I'm guessing that you meant to do the following:

u(i,j) = 2 * pi * pi * bsxfun(@times, sin(pi * i), sin(pi * j)');

Alternatively, you can use basic matrix multiplication to produce an (N-2)-by-(N-2) like so:

u(i, j) = 2 * pi * pi * sin(pi * i') * sin(pi * j); %// Note the transpose

P.S: it is recommended not to use "i" and "j" as names for variables.