Problem with quad2d integration

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

Answer 1

[Andrew Newell]

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)

[Max]

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]

How about next time it happens you post the parameters so I can try to reproduce it?

[Max]

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]

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]

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]

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.