1

I am trying to integrate on a surface the following equation: intensity=(2*J1(z)/z)^2 with z=A*sqrt((x-mu1)^2+(y-mu2)^2), A(L) a constant of x and y and J1 the first order bessel function. To do so I use the dblquad function as below:

resultinf = dblquad(lambda r,phi:intensity(mu1,mu2,L,r,phi),0,inf,lambda phi:0,lambda phi:2*pi)

The only important parameters here are r and phi in polar coordinates (the others depends of other parameters unimportant here), with x=rcos(phi) and y=rsin(phi) But when I try to integrate the function I get this message:

C:\pyzo2013b\Lib\pyzo-packages\scipy\integrate\quadpack.py:289: UserWarning: The maximum number of subdivisions (50) has been achieved. If increasing the limit yields no improvement it is advised to analyze the integrand in order to determine the difficulties. If the position of a local difficulty can be determined (singularity, discontinuity) one will probably gain from splitting up the interval and calling the integrator on the subranges. Perhaps a special-purpose integrator should be used.
warnings.warn(msg)

And then a completely unacurrate result followed by:

C:\pyzo2013b\Lib\pyzo-packages\scipy\integrate\quadpack.py:289: UserWarning: The integral is probably divergent, or slowly convergent. warnings.warn(msg)

I do understand the meaning of this messages but I got two questions:

  • Is there any mean for me to avoid this subdivision error other than just dividing my integrated intervalls in smaller segment (I'd like to check my other results by comparing them to the norm over an infinite domain and I won't be able to do so if I can't integrate over an infinite domain properly)? Maybe with a special purpose integrator? But I don't know what it is or how to use them.
  • Why do I get a warning about a divergent integral or a singularity knowing that J1(z)/(z) converge to 1 when z converge to zero (just like a sinus cardinal)?

Does anybody has an answer?

Here is the complete useful lines of codes (all the others parameters are defined otherwise):

def intensity(mu1,mu2,L,r,phi):# distribution area for a diffracted beam
   x=r*cos(phi)
   y=r*sin(phi)
   X=x-mu1
   Y=y-mu2
   R=sqrt(X**2+Y**2)
   scaled_R = R*Dt *pi/(lambd*L)
   return (4*(special.jv(1,scaled_R)**2/scaled_R**2)


resultinf = dblquad(lambda r,phi:intensity(mu1,mu2,L,r,phi),0,inf,lambda phi:0,lambda phi:2*pi)
print(resultinf)

(I have modified it on the advise of gboffi for the sake of a better understanding of the function.)

Vincent
  • 26
  • 5
  • Am I wrong or your results do not depend on `mu1` or `mu2`? if you name the (scaled) distance from `P = (mu1, mu2)` with `R` the integrand is, except for powers, `jv(1, R)/R` and does not depend on `phi`, hence you can avoid the double integral. – gboffi Jul 24 '17 at 16:21
  • @gboffi mu1 and mu2 are supposed to be here to move the center of my curve and are generated randomly according to a normal distribution thanks to the function random.gauss which itself depends of the time with sigma(t). In fact at the beginning I was integrating on X and Y but I realised that I would be advise to use polar coordinates to reduce the area of integration. But now I'm not sure that mu 1 and 2 only interact with the beam (curve) center, I may have forgotten to transform them with polar coordinates. – Vincent Jul 24 '17 at 18:42
  • The idea here in one dimension is to control the beam center (red curve) [link](https://drive.google.com/file/d/0B20eJey-YBadNmJmRHpsa3hENEE/view?usp=sharing) with mu the mean value. – Vincent Jul 24 '17 at 18:46
  • Let's try to simplify your expression: `x, y = r cos phi, r sin phi` ፨ `X, Y = x-mu1, y-mu2` ፨ `R = sqrt(X**2+Y**2` ፨ `scaled_R = R*Dt *π/(λ*L)` ፨ `integrand = 4*jv(1, scaled_R)**2 / scaled_R**2` ፨ the integrand depends only on the distance from what you've called the centre of the beam. Is it possible that what you have written is not what you meant to write? – gboffi Jul 24 '17 at 20:41
  • Well the idea was to take the 1D Franhofer diffraction equation for a circular aperture and to use it in 2 dimensions. The original equation is given by `integrand = 4*jv(1, scale_R)**2 / scaled_R**2` with `scaled_R=d*Dt *π/(λ*L)`. Now instead of d, I use R as you described it. I could just then integrate on x and y and get this [link](https://drive.google.com/file/d/0B20eJey-YBadU1V6VnlLRWlEb3c/view?usp=sharing) result for the function shape (here mu1=mu2=0). But I choose to integrate on the polar coordinates. The function value does not depends of phi but I need to integrate on a surface – Vincent Jul 25 '17 at 09:51
  • Plus the beam center moves with the time (random walk) so when mu1=mu2=0 this is true but as soon as the center displace the equation does depends of phi (I integrate on a circle centered in (0,0)). – Vincent Jul 25 '17 at 09:55
  • `∬ f(r, ϕ) dA = ∬ f(r, ϕ) dr r dϕ = ∫ r dr ∫ f(r, ϕ) dϕ` ፨ but `f(r, ϕ)` doesn't depend on `ϕ`, so it's a constant w/r to the inner integral. ፨ The extremes of integration of the inner integral are `ϕ = 0` and `ϕ = 2 π`, so we can perform the inner integration and write: `∫ r dr ∫ f(r, ϕ) dϕ= ∫ 2 π r f(r) dr` were in our last result the extremes of integration are, of course, `r=0` and `r=infty` – gboffi Jul 26 '17 at 09:03
  • I don't understand why you say that the function doesn't depend on phi. For me R=sqrt(X^2+Y^2)=sqrt(r^2+µ1^2+µ2^2-2r(µ1cos(phi)+µ2sin(phi))) and you can't get rid of the phi dependence unless µ1=µ2=0 (I may have made a mistake somewhere, but I don't see it for the moment) . For any other values the center is not in (0,0) and thus when I integrate on a circle center on (0,0) there is a dependance with the angle. – Vincent Jul 26 '17 at 14:11

0 Answers0