Applying composite midpoint rule in two dimensions

9 次查看（过去 30 天）
K_Dot
K_Dot 2017-7-17
回答： Saurabh Gupta 2017-7-20
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?

请先登录，再回答此问题。

回答（1 个）

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

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by