Applying composite midpoint rule in two dimensions
1 view (last 30 days)
Show older comments
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
Answers (1)
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
See Also
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!