numeric double integral take too long to compute

Hi, Im having a problem computing a double integral numerically when i determine the integration borders to be very big. the function needs a huge number of iterations to compute, taking couple of hours for each time. my code has it looped over, as i need it to compute a couple hundred times (the goal it to compute a near-field approximation of interference of a wave function on a razor edge at a distance D. thats why i need the integration borders to be very big in relation to lambda).

hers the code: x,y are the variables, and xp(i), yp(j) are numbers determined for each loop cycle. D, lambda are constants.

r = @(x,y) sqrt((x-xp(i)).^2+(y-yp(j)).^2+D^2);
fun = @(x,y) amp * exp(w .* r(x,y)) ./ sqrt(r(x,y));
A = quad2d(fun, -10^5*lambda, 10^5*lambda, 0, 10^5*lambda,'MaxFunEvals',10^7);

what can i do?

thanks

[John D'Errico]

Without looking carefully at your kernel, I'd wonder if a Gauss-Hermite integration would be a better choice than an adaptive quadrature.

[David Goodmanson]

Hi amit, could you provide some information on what the variable w looks like?

[John D'Errico]

Note that amp is completely irrelevant. However, at the very least, removing it from the integration saves a multiply.

[amit nestor]

of course, i forgot about those two variables. but they dont make a big difference - w is a complex given constant, and amp=1, not relevant in this case. doesnt make a big difference if i take amp out or not.

this integral is enclosed in two for loops:

for i = 1:100
  for j = 1:100
      r = @(x,y) sqrt((x-xp(i)).^2+(y-yp(j)).^2+D^2);
      fun = @(x,y) amp * exp(w .* r(x,y)) ./ sqrt(r(x,y));
      A = quad2d(fun, -10^4*lambda, 10^4*lambda, 0,10^3*lambda,'MaxFunEvals',10^7);

thats the entire relevant code piece - each integral is calculated for a point (xp,yp) between 1 to 100.

[David Goodmanson]

Hi amit, This does not look like the standard Fresnel diffraction calculation. Do you have a reference you could supply for this calculation?

Assuming that the limits have no significance other than having to be very large (which is tied to the value of w), there appears to be no reason why this integral will depend on xp(i).

[amit nestor]

sure, we based our model on the article: http://www.people.fas.harvard.edu/~djmorin/waves/interference.pdf on the near field approximation (page 26, formula 48). this is for a plane wave in one dimension. r(y) representing the distance from the center of the screen to any point y on the integration area, which is the slit. we modeled it for any point (p) on the screen, so we had to calculate the vector r(y) for each point p. we represented it by yp(i) for each point in the desired area, and basic geometry led us to that expression of r(y):

r = @(y) sqrt((y-yp(i)).^2 + D^2)

now we wanted to model it for 2 dimensional slit and screen. now instead of using just y and yp we had to factor in x and xp. this led us to the expression above for r(x,y), and the integration simply takes account also x-dimension, given the boundaries.

in each calculation xp(i) and yp(j) are given points to which the integral computes and returns a value, that is the field in that point on the screen. the integration is only done on x and y as variables. w is also a constant, depending on the wave length alone (its not omega, just a letter we used for comfort).

hope it helps, thanks

[David Goodmanson]

I took a look at the reference, which is about two dimensional situations such as fig. 26 or (in your case) fig. 35, i.e. with propagation to the right and the aperture determined by the vertical coordinate y. The slits are infinitely long in the third dimension z, perpendicular to the page. It's a 2d problem, nothing varies in z, and you only need a single integral for the solution as in eqn. 48.

Fresnel diffraction is treated a bit differently than this in most textbooks so it might be a good idea to take a look at some of them. The approximation they arrive at is simpler than eqn. 48.

If you include the third dimension, that would mean you have, say, a rectangular aperture in a screen, with the height of the aperture as in Fig. 26 and the width determined by z. Your limits of integration appear to be approximations to infinity rather than the boundaries of a very, very, very, very, large aperture. If you do have an infinitely wide aperture in the third dimension you don't need anything more than the single integral.

If you have a reasonably sized rectangular aperture in a plate then you do need the third dimension and you have to do a double integral but it's not in the same form that you borrowed from eqn. 48. For a 2d problem the denominator ~ sqrt(r), but for a 3d problem the denominator ~ r. I think it would make sense not to rely entirely on Morin but to get the perspective of a more standard treatment from commonly used textbooks. Or even the 'Fresnel Diffraction' discussion in wikipedia.