Problem with quad2d integration

Max
2012 1 25
1 个回答
5 次查看（30 天）

0 个投票

Hi!
I have got a problem with the quad2d function, integrating a double-integral. My function is a function in two scalar variables. I get the error message regarding a dimension missmatch. The function is:
f=@(xq,zq)exp(1i*2*pi/lambda*(sqrt((xp0-xq)^2+yp0^2+(zp0-zq)^2)+sqrt((xp-xq)^2+yp^2+(zp-zq)^2)))*(yp0/(sqrt((xp-xq)^2+yp^2+(zp-zq)^2))-yp/sqrt((xp0-xq)^2+yp0^2+(zp0-zq)^2))/(sqrt((xp0-xq)^2+yp0^2+(zp0-zq)^2)*sqrt((xp-xq)^2+yp^2+(zp-zq)^2));
where xp0,yp0,zp0 and xp,yp,zp are scalar constants.pi and lambda are constant as well. the function is basically and exponential function multiplied with a dotproduct of 2 vectors. I know that my error has something to do with the function not being able to handle vector inputs. so i tried using elemetwise power .^ and elemtentwise .* and ./ but then i get an error:
Warning: Reached the maximum number of function evaluations (2000). The result fails
the global error test.
> In quad2d at 241
In test2 at 22
where test2() is just the function that declares the f function and doing the integration quad2d(f,xq,zq,0,1,0,2). I actually put the elementwise . operator almost randomly because i dont really know how to vectorize my scalar function.
I'd be gratetful for any help! Max

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

回答（1 个）

Andrew Newell
Andrew Newell 2012-1-25

1 个投票

When in doubt, just put a dot in front of all ^'s, *'s and /'s. Also, next time you ask a question, it would help if you provided code with numerical values for all parameters. This code gives an answer:
lambda=rand;
xp0=rand;
xp=rand;
yp0=rand;
yp=rand;
zp0=rand;
zp=rand;
f = @(xq,zq) exp(1i.*2.*pi./lambda.* ...
(sqrt((xp0-xq).^2+yp0.^2+(zp0-zq).^2)+sqrt((xp-xq).^2+yp.^2+(zp-zq).^2))) ...
.*(yp0./(sqrt((xp-xq).^2+yp.^2+(zp-zq).^2))-yp./sqrt((xp0-xq).^2+yp0.^2+(zp0-zq).^2)) ...
./(sqrt((xp0-xq).^2+yp0.^2+(zp0-zq).^2).*sqrt((xp-xq).^2+yp.^2+(zp-zq).^2));
quad2d(f,0,1,0,2)

6 个评论
显示 4更早的评论 隐藏 4更早的评论

Max
Max 2012-1-27
hi,
thanks for the quick response! i got it to work now.
but still, sometimes i get the global error test fail warning mentioned above:
<code>
Warning: Reached the maximum number of function evaluations (2000). The result fails
the global error test.
> In quad2d at 241
In test2 at 22
</code>
i cant figure out what it is relatet to. i guess it changes with the integration area. maybe this is due to the rapdily oscillating exp-function.
any ideas how i can avoid this?
max
Andrew Newell
Andrew Newell 2012-1-27
How about next time it happens you post the parameters so I can try to reproduce it?
Max
Max 2012-1-28
hi,
for example try:
[code]
function test3()
a=1;
P0=[0,-1,0];
xp0=P0(1);
yp0=P0(2);
zp0=P0(3);
lambda=4*10^(-7);
P=[0,0,0];
xp=P(1);
yp=P(2);
zp=P(3);
f=@(xq,zq)exp(1i.*2.*pi./lambda.*(sqrt((xp0-xq).^2+yp0.^2+...
(zp0-zq).^2)+sqrt((xp-xq).^2+yp.^2+(zp-zq).^2)));
U=(quad2d(f,-10*10^(-6),10*10^(-3),-10*10^(-3),10*10^(-6)));
end
[/code]
in the MATLAB help it says, that there might be a singularity in the integration area, but dont see a connection between the exp fnction and a singularity since the error message doesnt seem to show up at a certian bound.
Andrew Newell
Andrew Newell 2012-1-28
It's always a good idea to plot the function you're going to integrate. Try adding this code inside your function:
xlb = -1e-5; xub = 1e-2; ylb = -1e-2; yub=1e-5;
x = linspace(xlb,xub,100); y = linspace(ylb,yub,100);
[X,Y] = meshgrid(x,y);
Z = f(X,Y);
surf(X,Y,real(Z))
You'll see a surface that looks like a bed of nails. Why? because the real and imaginary parts of f are a cosine and sine of a function that is multiplied by 1/lambda = 0.25e7. It oscillates so fast that you can't possibly integrate it. Try using lambda close to 1, or even as low as 4e-4, and you'll have no problem.
Max
Max 2012-1-28
ok thanks! i thought it had something to do with the oszillation. but scaling the lambda to 1 means that i have to adjust all other paramters to 10^7 since lambda stands for the wavelenght of an optic source. when i scale all other paramters up, including the integration area, doenst this reproduce the same problem? and additional to that, the time needed for the calculation is going to be a lot longer? since my program is already really slow, this would be a big problem. but anyways, i will try adjusting the parameters tomorow.
Andrew Newell
Andrew Newell 2012-1-28
It won't help to adjust all the parameters - you'll still have the same rapidly oscillating function. You have a function that oscillates rapidly but probably averages to near zero, which means that the roundoff errors will kill any numerical approximation. You'll need to do some careful thinking.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

提问:

Max
Max
2012-1-25

Community Treasure Hunt

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

Start Hunting!

Translated by