Интеграция бесселевой функции в python: проблема подразделения

Question

Я пытаюсь интегрировать на поверхности следующее уравнение: интенсивность =(2*J1(z)/z)^2 с z=A*sqrt((x-mu1)^2+(y-mu2)^2), A(L) постоянная x и y и J1 функция бесселева первого порядка. Для этого я использую функцию dblquad, как показано ниже:

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

Единственными важными параметрами здесь являются r и phi в полярных координатах (остальные зависят от других неважных здесь параметров), причем x = rcos (phi) и y = r sin (phi). Но когда я пытаюсь интегрировать функцию, я получаю это сообщение:

C: \ pyzo2013b \ Lib \ pyzo-packages \ scipy \ integrate \ quadpack.py: 289: Предупреждение пользователя: достигнуто максимальное количество подразделений (50). Если увеличение лимита не приводит к улучшению, рекомендуется проанализировать подынтегральное выражение, чтобы определить трудности. Если можно определить положение локальной сложности (сингулярность, разрыв), то, вероятно, выиграет от разбиения интервала и вызова интегратора на поддиапазонах. Возможно, следует использовать специальный интегратор.
warnings.warn (MSG)

И тогда совершенно неточный результат, за которым следуют:

C: \ pyzo2013b \ Lib \ pyzo-packages \ scipy \ integrate \ quadpack.py: 289: Предупреждение пользователя: Интеграл, вероятно, расходится или медленно сходится. warnings.warn(MSG)

Я понимаю значение этих сообщений, но у меня есть два вопроса:

• Есть ли какой-то смысл для меня избежать этой ошибки подразделения, кроме простого деления моих интегрированных интервалов на меньший сегмент (я хотел бы проверить другие мои результаты, сравнив их с нормой по бесконечной области, и я не смогу сделать так что, если я не могу должным образом интегрироваться в бесконечную область)? Может быть, со специальным интегратором? Но я не знаю, что это такое или как их использовать.
• Почему я получаю предупреждение о расходящемся интеграле или сингулярности, зная, что J1(z)/(z) сходятся к 1, когда z сходится к нулю (как синусовый кардинал)?

У кого-нибудь есть ответ?

Вот полные полезные строки кодов (все остальные параметры определены иначе):

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)

(Я изменил его по совету gboffi для лучшего понимания функции.)