Applying composite midpoint rule in two dimensions

I am trying to use the composite midpoint rule to evaluate a two dimensional integral as follows:

lambda = 2*pi
h = lambda/20; % step size
xx = 0:h:20; % x axis with step size h
yy = 0:h:20; % y axis with step size h
[X,Y] = meshgrid(xx,yy); % create meshgrid
rho = [X,Y]; % position vector
rhos = [lambda/2 10*lambda]; % location of a point
phi = @(x)(x>0).*exp(-1./x.^2);
K = @(X,Y) phi(X).*phi(1-X).*phi(Y).*phi(1-Y); 
Chi = @(X,Y) (X>0).*((K(X,Y)/kb).^2-1);
M = 20; 
[M1, M2] = meshgrid(M,M); % create mesh grid
h = [lambda lambda]./[M1 M2];

us = zeros(M1,M2); % preallocate memory
c = zeros(M1,M2);
for n = 1:M1
  for m = 1:M2
      c(n,m) = ([n m]-(1/2)).*h;
      us = -(kb^2/16)*h.*besselh(0,2,kb*norm(rho(n,m)-ck(n,m))).*Chi(rho).*...
          besselh(0,2,kb*norm(rho(n,m)-rhos));
  end
end

where in the composite midpoint rule we have h = (b-a)/m, but here my a = 0 and b = (lambda,lambda) and m = (M1, M2) is two dimensional.

I get the error

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

I was wondering if someone could help me overcome this problem and compute via the composite midpoint rule?

Answer 1

[Saurabh Gupta]

The error occurs at the following line

c(n,m) = ([n m]-(1/2)).*h;

The reason is that lhs is a scalar (single element of a matrix), whereas rhs is a 1x2 vector. You need to match the dimensions of lhs and rhs, though the actual solution would depend on what you expect of this operation.

If you want rhs to return the result as a scalar, you need to update rhs expression to return a scalar value.

On the other hand, if you want to store the 1x2 vector in the variable on lhs, you need to supply a 1x2 vector depending upon your data structure design. It could be something like

c(n,m:m+1) = ([n m]-(1/2)).*h;

or you may need a 3D data structure like this

c = zeros(M1,M2,2);
...
c(n,m,:) = ([n m]-(1/2)).*h;