Высоко колебательный 2D-продукт Bessel Function Integrand - PullRequest
Новая
       109

Высоко колебательный 2D-продукт Bessel Function Integrand

0 голосов
/ 20 марта 2020

Моя задача - численно интегрировать следующую функцию, я бы выложил уравнение, но мне не хватает «репутации». Это можно найти в приведенном ниже коде.

Согласно scipy.integrate.nquad подынтегральная функция сильно колеблется в обеих переменных. Я хочу вычислить эту функцию и построить ее относительно z. Проблема в том, что этот интеграл не сходится к нулю, как ожидалось. График, который я имею, показывает, что он сходится асимптотически, однако мы считаем, что он должен сходиться гораздо быстрее, но при больших z подынтегральное выражение колеблется быстрее, что приводит к нашей проблеме. Я не могу опубликовать сюжет, потому что мне не хватает «репутации».

Я пытался использовать mpmath.quados c, но в документации не показано, как можно сделать двойной интеграл с квадро c , Я пытался установить предельные точки очень высоко (5000) в scipy.integrate.nquad, но это не помогло. Мой код следующий:

Мой вопрос: как лучше всего с этим справиться? т.е. какие методы я мог бы использовать, чтобы получить правильное поведение этого интеграла? Я видел некоторые другие сообщения о стекловом потоке об этой проблеме, но тот факт, что я вычисляю двойной интеграл, я думаю, усложняет ситуацию. 

def IntegrandJp(z, dl):
return integrate.nquad(
    lambda y,x: - x * scipy.special.jn(0, x * z)*np.exp((-x**2/2))*scipy.special.jn(0, x * y)*np.exp((-( y ** 2) / 2) - dl * y),[[0, np.inf],[0, np.inf]],opts={'limit':5000})[0]

def parallel_runs_Jp(z):
for i, dl in enumerate(DL):
    p = multiprocessing.Pool(8)
    IJp_dl = partial(IntegrandJp, dl)
    Jp_dict[dl] = p.map(IJp_dl, z)
return Jp_dict

dl - это просто реальный параметр, который исходит из От 0 до 1,5 с шагом 0,1. Он отвечает за семейство кривых на участке.

...