MATLAB Answers

K_Dot
0

Applying composite midpoint rule in two dimensions

Asked by K_Dot
on 17 Jul 2017
Latest activity Answered by Saurabh Gupta on 20 Jul 2017
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?

  0 Comments

Sign in to comment.

1 Answer

Answer by 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;

  0 Comments

Sign in to comment.