Applying composite midpoint rule in two dimensions

12 views (last 30 days)
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?

Answers (1)

Saurabh Gupta
Saurabh Gupta on 20 Jul 2017
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;

Categories

Find more on Creating and Concatenating Matrices in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!